ベクトルの内積、行列ベクトル積、行列積の区間演算はとても重要。すなわち $x, y\in\mathbb{IF}^n$, $A, B\in\mathbb{IF}^{n\times n}$ に対して、それぞれ$x^Ty \subset \alpha$, $Ax\subset z$, $AB\subset C$となる $$ \alpha\in\mathbb{IF},\quad z\in\mathbb{IF}^n,\quad C\in\mathbb{IF}^{n\times n} $$ を求める方法を考える。
最も素朴な方法は各演算を区間演算に変更した方式である。
$$ \alpha = \sum_{k = 1}^n x_ky_k $$$$ z_i = \sum_{k = 1}^n A_{i,k}x_k\quad (1\le i\le n) $$$$ C_{ij} = \sum_{k+1}^n A_{i,k}B_{k,j}\quad (1\le i,j\le n) $$特にJuliaのIntervalArithmetic.jlパッケージでは、この方法で区間演算が行なわれている。しかし、この方法だとプログラムの最適化が難しく、計算機の性能を十分に引き出した実装が難しく、計算時間がかかってしまう。
using IntervalArithmetic, BenchmarkTools, StaticArrays
n = 1000;
x = randn(n);
y = randn(n);
x_int = map(Interval, x);
y_int = map(Interval, y);
@btime alpha = dot(x_int,y_int) # Level 1
# x_int = IntervalBox(map(Interval, x));
# y_int = IntervalBox(map(Interval, y));
# @btime alpha = dot(x_int,y_int) # Level 1
# x_int = SVector{n}(map(Interval, x));
# y_int = SVector{n}(map(Interval, y));
# @btime alpha = dot(x_int,y_int) # Level 1
# setrounding(Interval, :fast)
using LinearAlgebra
n = 1000; A = randn(n,n); x = randn(n);
A_int = map(Interval, A); x_int = map(Interval, x); y = similar(x_int);
@time z = A_int*x_int; # Level 2
@time mul!(y,A_int,x_int) # Level 2
# x_int = IntervalBox(x_int);
# @time z = A_int*x_int; # Level 2
# @time mul!(y,A_int,x_int) # Level 2
println(y ⊆ z)
println(z ⊆ y)
# @btime z = $A_int*$x_int;
n = 100;
A = randn(n,n);
B = randn(n,n);
A_int = map(Interval, A);
B_int = map(Interval, B);
@time C = A_int*B_int; # Level 3
@time C = A_int*B_int;
@time C = A_int*B_int;
# A_int = SMatrix{n,n}(A_int);
# B_int = SMatrix{n,n}(B_int);
# @time C = A_int*B_int; # Level 3
# @time C = A_int*B_int;
# @time C = A_int*B_int;
このくらいの計算スピードならば、体感として、そこまで遅い印象はないが、そこそこ大きなサイズの行列積（Level 3 BLAS）は計算時間がかかってしまう。例えば、
n = 1000; A = randn(n,n); B = randn(n,n);
A_int = map(Interval, A); B_int = map(Interval, B);
@btime C = $A_int*$B_int; # Level 3
注意 BenchmarkToolsの@btimeや@benckmarkマクロを使うときは、
@btime A*B;
とせずに
@btime $A*$B;
と変数名に$をつける。これは計測時間に変数の呼び出しを含めるか指定するもので、前者が「変数の呼び出し時間を含めた計測」、後者は「行列積の計算時間の計測」を意味している。
さて、区間行列積(Level 3)を高速化したい。考えられる方法は
@simd命令）Tullioパッケージを使うと行列積を自動的に並列化して計算してくれる。小さいサイズの行列積はそこまで高速化されないが、ある程度行列サイズが大きくなると区間行列積の自動並列によって、計算時間が短くなることが観測される。ただし、JULIA_NUM_THREADSという環境変数で使用するスレッド数をあらかじめ定義しておく必要がある。例えばmacOSでzshというシェルを使っているなら、~/.zshrcというファイル内に
export JULIA_NUM_THREADS=8
と計算スレッド数を指定しておく必要がある。CPUのコア数分だけ指定したら良い？
using LoopVectorization, BenchmarkTools, IntervalArithmetic, Tullio
tullio_gemm(A, B) = @tullio C[i,j] := A[i,k] * B[k,j]
n = 100;
A = map(Interval,randn(n,n));
B = map(Interval,randn(n,n));
@btime C = $A*$B;
@btime C1 = tullio_gemm($A, $B);
# Threads.nthreads()
n = 1000;
A = map(Interval,randn(n,n)); B = map(Interval,randn(n,n));
# @btime C = A*B; # ≈ 40 sec
# setrounding(Interval, :tight)
@btime C1 = tullio_gemm(A,B); # ≈ 6 sec
次、行列、ベクトルの配列に
StaticArraysは小さいサイズでは多少計算時間が変わるが、要素数が多くなるほどコンパイルに時間がかかり、全体の経過時間は長くなってしまった。要素数が100以下ならStaticArraysが使えるとWebサイトに書いてありました。
残るは
まずJuliaの中で、BLASを使うにはLinearAlgebraパッケージが必要。
using LinearAlgebra, BenchmarkTools
# versioninfo()
# BLAS.vendor()
n = 5000;
A, B, C = randn(n,n), randn(n,n), zeros(n,n);
# @btime C = $A * $B;
@btime mul!(C, A, B);
このようにBLASを使うと、物凄いスピードで数値計算ができる。ただし、行列の内部の型がFloat64に限られるのと、丸め方向を変更した計算ができないため、区間演算の結果が粗くなる（区間幅が増加する）。以下、丸め方向を変更しない区間行列積を実装する。
定義 $\mathbf{u}=2^{-53}$を倍精度の単位相対丸めとする。$\mathbf{S}_{\min}=2^{-1074}$を倍精度浮動小数点数の正の最小数とする。$\mathbf{F}_{\min}=2^{-1022}$を倍精度浮動小数点数の正規化数の正の最小数とする。
注意 単位相対丸めは$1$と$1$よりも小さい最大の浮動小数点数との差を表す。つまり
$$ 1-\mathbf{u}<a<1 $$となる$a\in\mathbb{F}$は存在しない。
注意 単位相対丸めとは別に計算機イプシロン（unit roundoff, machine epsilon）という単位もあり、こちらは$1$と$1$よりも大きい最小の浮動小数点数との差を表し、64bit浮動小数点数では$2^{-52}$でこちらを$\mathbf{u}$と書く流儀もある（こっちが主流？）ので、注意が必要。Juliaではeps(Float64)とすると$2^{-52}$の数値が得られる。
u = 2.0^-53;
s_min = 2.0^-1074;
f_min = 2.0^-1022;
println(bitstring(1.0))
println(bitstring(1.0-u))
println(bitstring(s_min))
println(bitstring(f_min))
println(eps(Float64))
println(2*u)
定義 実数$a\in\mathbb{R}$に対するUnit in the First Place (ufp)は次のように定める。
$$ \mbox{ufp}(a):= 2^{\lfloor\log_2|a|\rfloor} (a\neq 0),\quad \mbox{ufp}(0)=0. $$関数 $\mbox{ufp}$ は実数$a$を2進数表記した先頭ビットを返す。4回の浮動小数点演算で求めることができる。関数 $\mbox{ufp}$ を使うと $\mbox{ufp}(a)\le a < 2\mbox{ufp}$ が成立する。
function ufp(P)
u = 2.0^(-53);
ϕ = 2.0^52 + 1;
q = ϕ * P;
T = (1 - u)*q;
return abs(q - T);
end
p = 1.0 - u;
println(bitstring(p))
println(bitstring(ufp(p)))
println(bitstring(2*ufp(p)))
定義 $c\in\mathbb{R}$として、
これらを使うと $a$, $b\in\mathbb{F}$, $\circ\in\{+,-,\times,\div\}$ で
$$ \mbox{pred}(\mbox{fl}(a\circ b))<a\circ b<\mbox{succ}(\mbox{fl}(a\circ b)) $$が成り立つ。ここで $\mbox{fl}:\mathbb{R}\to\mathbb{F}$ は入力の実数を最も近い浮動小数点数に写す写像で、$\mbox{fl}$ 内の演算は全て浮動小数点演算で計算されているとする。
また、pred関数, succ関数はベクトル・行列においても各要素に対するsucc, predを考える事で拡張することが出来る。以下では、文献1のRumpによる方法を実装する。これは複雑だが、なるべく非正規化数の演算を避けるように設計されている。
function succ(c)
s_min = 2.0^-1074;
u = 2.0^-53;
ϕ = u * (1.0 + 2.0 * u);
if abs(c) >= (1. / 2.) * u^(-2) * s_min # 2^(-969)(Float64)
e = ϕ * abs(c);
succ = c + e;
elseif abs(c) < (1. / 2.) * u^(-1) * s_min # 2^(-1022)(Float64)
succ = c + s_min;
else
C = u^(-1) * c;
e = ϕ * abs(C);
succ = (C + e) * u;
end
return succ
end
function pred(c)
s_min = 2.0^-1074;
u = 2.0^-53;
ϕ = u * (1.0 + 2.0 * u);
if abs(c) >= (1. / 2.) * u^(-2) * s_min # 2^(-969)(Float64)
e = ϕ * abs(c);
pred = c - e;
elseif abs(c) < (1. / 2.) * u^(-1) * s_min # 2^(-1022)(Float64)
pred = c - s_min;
else
C = u^(-1) * c;
e = ϕ * abs(C);
pred = (C - e) * u;
end
return pred
end
a = 0.1;
println(bitstring(succ(a)))
println(bitstring(a))
println(bitstring(pred(a)))
ベクトルの総和と内積のufpを用いた事前誤差解析を応用して、行列 $A\in\mathbb{F}^{m\times n}$, $B\in\mathbb{F}^{n\times p}$ の積 $AB$ を包含する中心半径型区間行列 $\langle C,R\rangle\in\mathbb{IF}^{m\times p}$ を計算する。
$$ C=\mbox{fl}(AB),\quad R=\mbox{fl}((n+2)\mathbf{u}\cdot \mbox{ufp}(|A||B|)) + \mathbf{F}_{\min}\cdot e_{m}e_{p}^T). $$ただし $e_n=(1,\dots,1)^T\in\mathbb{F}^{n}$ で $n$ は $2(n+\mathbf{u})<1$を満たすとする。このとき必要な行列積の計算は、$\mbox{fl}(AB)$, $\mbox{fl}(|A||B|)$ の2回である。これらをBLASを用いて計算することにより、高速な数値計算が実行可能である。そして、浮動小数点演算を行った結果 $C$ とその誤差半径を $R$ とし、行列積の真値を包含する。
function mm_ufp(A_mid, B_mid)
u = 2.0^(-53);
realmin = 2.0^(-1022);
n = size(A_mid,2);
if(2*(n+2)*u>=1)
error("mm_ufp is failed!(2(n+2)u>=1)")
end
C_mid = A_mid * B_mid;
C_rad = (n+2) * u * ufp.(abs.(A_mid)*abs.(B_mid)) .+ realmin;
return C_mid, C_rad;
end
using IntervalArithmetic, LinearAlgebra
n = 1000;
A = randn(n,n);
B = randn(n,n);
A_int = map(Interval, A);
B_int = map(Interval, B);
A_mid = mid.(A_int);
B_mid = mid.(B_int);
@btime C_mid, C_rad = mm_ufp(A_mid, B_mid);
C_mid, C_rad = mm_ufp(A_mid, B_mid);
C_int = C_mid .± C_rad;
注意 @tullioマクロで並列化した区間演算（8コア！）で約6秒かかる区間行列積が、BLASを使うと、約0.02秒と300倍程度高速化される。これがBLASの力だー！しかし、BLASの計算はFloat64のみサポートされるため、64ビットの浮動小数点数を用いた区間演算のとき以外は、このままでは高速化されない点を注意したい。例えば端点にBigFloatのような多倍長数値を用いたいときなどは、残念ながらBLASを用いたの高速な計算が出来ないのが現状である。また、ufpを用いる区間演算は過大評価になるため、区間幅の増大が避けられない事（概ね10倍大きな区間幅になる）にも注意が必要である。
C1 = tullio_gemm(A_int, B_int);
println(maximum(radius.(C1[:])))
println(maximum(C_rad[:]))
次に区間行列 $A=\langle A_m,A_r\rangle\in\mathbb{IF}^{m\times n}$, $B=\langle B_m,B_r\rangle\in\mathbb{IF}^{n\times p}$ の積 $AB$ を包含する区間行列 $\langle C,T\rangle\in\mathbb{IF}^{m\times p}$ を計算することを考える。はじめに区間行列積 $AB$ は
$$ AB\subseteq \langle A_mB_m,\tilde{T}\rangle,\quad \tilde{T}:=|A_m|B_r+A_r(|B_m|+B_r) $$と包含できる。丸め方向の変更ができる場合は、この中心と半径を丸めの向きを変えることで計算し厳密に包含することができる。しかし今回は丸め方向の変更ができないため、事前誤差評価によって中心と半径を厳密に包含する必要がある。いま $n$ は $2(n+\mathbf{u})<1$ を満たすとし、区間行列の中心同士の積 $A_mB_m$ は
$$ C=\mbox{fl}(A_mB_m),\quad R=\mbox{fl}((n+2)\mathbf{u}\cdot \mbox{ufp}(|A_m||B_m|)) + \mathbf{F}_{\min}\cdot e_{m}e_{p}^T). $$として $\langle C,R\rangle(=\mbox{mm_ufp}(A_m,B_m))$ で包含できる。次に $|A_m|B_r$ の上限を考える。
\begin{align*} &T_1:=\mbox{fl}(|A_m|B_r)\\ &||A_m|B_r-\mbox{fl}(|A_m|B_r)|\le \mbox{fl}((n+2)\mathbf{u}\cdot\mbox{ufp}(T_1) + \mathbf{F}_{\min}\cdot e_{m}e_{p}^T)=:T_2\\ &|A_m|B_r\le T_1+T_2 \end{align*}さらに $|B_m|+B_r$ は
$$ |B_m|+B_r\le\mbox{succ}(\mbox{fl}(|B_m|+B_r))=:T_3 $$で抑えられる。さらに $A_rT_3$ の上限は
\begin{align*} &T_4:=\mbox{fl}(A_rT_3)\\ &|A_rT_3-\mbox{fl}(A_rT_3)|\le \mbox{fl}((n+2)\mathbf{u}\cdot\mbox{ufp}(T_4) + \mathbf{F}_{\min}\cdot e_{m}e_{p}^T)=:T_5\\ &A_rT_3\le T_4+T_5. \end{align*}最後に誤差半径 $T$ は $T_1+T_2+T_3+T_4+T_5+R$の上限で得られ、
\begin{align*} &T_1+T_2+T_3+T_4+T_5+R\\ &\le \mbox{fl}(T_1+T_2+T_3+T_4+T_5+R) + \mbox{fl}(4\mathbf{u}(T_1+T_2+T_3+T_4+T_5+R))\\ &\le \mbox{succ}(\mbox{fl}(T_1+T_2+T_3+T_4+T_5+R) + \mbox{fl}(4\mathbf{u}(T_1+T_2+T_3+T_4+T_5+R)))=:T. \end{align*}したがって、$AB\subseteq \langle C,T\rangle$ となる $C\in\mathbb{F}^{m\times p}$ と $T\in\mathbb{F}^{m\times p}$ が得られる。
function imm_ufp(A_mid, A_rad, B_mid, B_rad)
u = 2.0^(-53);
realmin = 2.0^(-1022);
n = size(A_mid,2);
if(2*(n+2)*u>=1)
error("mm_ufp is failed!(2(n+2)u>=1)")
end
# C, R = mm_ufp(A_mid,B_mid);
C_mid = A_mid * B_mid;
R = (n+2) * u * ufp.(abs.(A_mid)*abs.(B_mid)) .+ realmin;
# T_1, T_2 = mm_ufp(abs.(A_mid), B_rad);
T1 = abs.(A_mid) * B_rad;
T2 = (n+2)*u*ufp.(T1) .+ realmin;
# T_3 = succ.(abs.(B_mid)+B_rad);
T3 = succ.(abs.(B_mid)+B_rad)
# T_4, T_5 = mm_ufp(A_r, T_3);
T4 = A_rad * T3;
T5 = (n+2)*u*ufp.(T4) .+ realmin;
rad_sum = R + T1 + T2 + T4 + T5;
C_rad = succ.(rad_sum + 4*u*ufp.(rad_sum));
return C_mid, C_rad;
end
using IntervalArithmetic, LinearAlgebra
# n = 1000;
# A = randn(n,n);
# B = randn(n,n);
# A_int = map(Interval, A);
# B_int = map(Interval, B);
A_mid = mid.(A_int);
A_rad = radius.(A_int);
B_mid = mid.(B_int);
B_rad = radius.(B_int);
@btime C1_mid, C1_rad = imm_ufp(A_mid, A_rad, B_mid, B_rad);
C1_mid, C1_rad = imm_ufp(A_mid,A_rad,B_mid,B_rad);
C1_int = C1_mid .± C1_rad;
# CC = A_mid*B_mid .∈ C1;
# sum(CC[:])
# println(size(C_int))
# println(size(C1))
# println(maximum(radius.(C1[:])))
println(maximum(C_rad[:]))
println(maximum(C1_rad[:]))
# println(size(C_int))
注意 mm_ufpとimm_ufpは入力の行列 $A$, $B$ に応じて使い分けることを計算速度、誤差半径の両面からおすすめする。
本資料も筆者が学生の頃に精度保証付き数値計算を教えて下さった柏木雅英先生の「数値解析特論」の講義資料を参考にしています。また、丸め方向を変更しないBLASを使った区間行列積について森倉悠介先生から助言をいただき、imm_ufpのMATLABコードをご提供していただきました。ここに感謝申し上げます。
さらに、以下のような文献・Web ページ等を参考にこの文章は書いています。