行列計算における確率的誤差解析～ 行列指数関数の計算を例として ～2022年12月23日
日本応用数理学会
三部会連携「応用数理セミナー」
電気通信大学 情報理工学研究科
山本有作

はじめに• 行列指数関数 exp(A)– 任意の行列 A に対して収束– 応用： 微分方程式の数値解法• 数値計算法（exp(A) 自体を計算するもの）– 打切り型テイラー展開– 有理関数近似– 複素周回積分– 固有値分解／シューア分解– Scaling & Squaring

・ Aのノルムが小さいときは効率的（ODEの数値解法では exp(AΔt) を計算）

・ 行列乗算のみで計算可能（HPC の観点から最適）2本発表の目的• 打切り型テイラー展開による exp(A) の計算法を例題として，近年 Higham らにより提案された新しい確率的誤差解析1)の実際について解説する．
• 併せて，以下の点についても触れる．
– 混合精度計算の利用– ブロッククリロフ部分空間法における残差ギャップの評価• 背景– 誤差制御型数値計算への要求（速度と精度のトレードオフ）– 超高速・低精度（半精度など）の計算ハードウェアの登場
1) N. J. Higham & T. Mary: A New Approach to Probabilistic Rounding Error Analysis, SIAM J. Sci. Comput.,Vol. 41, No. 5, pp. A2815-A2835 (2019).
3目次• はじめに• 行列指数関数の計算手法• 誤差解析
– 打切り誤差– 確定的丸め誤差解析– 混合精度演算の場合– 確率的丸め誤差解析• 数値実験
• 事後誤差解析への応用– ブロッククリロフ部分空間法における残差ギャップの評価• おわりに4

テイラー展開に基づく行列指数関数の計算• 計算手順– 次の手順により，exp(A) の近似値 fn(A) を計算すると仮定• 注意– 後に述べる混合精度演算では，ある k1（≧2）に対し，k > k1 のときのBk の計算を低精度演算で行う．
– 計算にはホーナー法を使うほうが一般的であるが，低精度で計算する部分を後回しにしてその影響を局限するため，上記の計算法を用いた．
5

誤差解析• 打切り誤差– 無限精度で計算したときの打切り誤差は，次のように評価される．

– ここで，ノルムは p ノルム（p≧1）またはフロベニウスノルムである．
6

（確定的）丸め誤差解析• 準備– 浮動小数点演算で計算した量を，fl(・) で表す．
– 丸め誤差の単位を u とする．
• fl(a ・○ b) = (1 + δ)(a ・○ b)， ただし，○・ = +, –, *, /， | δ | ≦ u．

– γn = nu / (1 – nu) とおく．
• (1 + γm)(1 + γn) ≦ 1 + γm+n などの関係式が成り立つ．

– A = (aij) に対し，|aij| を要素とする行列を |A| で表す．

– 行列 A, B に対し，A≦B は要素ごとに不等式が成り立つことを表す．

• 定理： 行列乗算の誤差上界– A∈Rm×n, B∈Rn×p のとき，行列乗算 AB の丸め誤差の上界は次式で与えられる．
| fl(AB) – AB | ≦ γn| A | | B |.
7

Ak の誤差^• Ak = fl(Ak) とし，ある正数 εk–1 > 0 についてが成り立つと仮定する．このとき，まず，次の式が成り立つ．

これと行列乗算の丸め誤差の公式を使うと，第 k 回の乗算の誤差 Ak–1 が持っていた誤差

これを εk と置ける
8

Ak の誤差（続き）• εk に関する漸化式– 前ページの結果より，– ε1 = 0 に注意してこれを解くと，εk に対する次の上界が得られる．

(1 + γm)(1 + γn) ≦ 1 + γm+n• Ak の誤差– 以上より，Ak の丸め誤差の上界が次のように得られる．

9

1/k! 倍の誤差^

• Bk = (1/k!)Ak, Bk = fl(Bk) とすると，|fl(a/b) – a/b| ≦ (1+u)(a/b) を利用

前ページの結果を利用| A^k | ≦ | Ak | + | A^k – Ak | として前ページの結果を利用(1 + γm)(1 + γn)≦ 1 + γm+n10

級数の和の誤差• n + 1 個の行列の加算に対する丸め誤差の上界は，次の公式により与えられる．

• ここで，Ck = B^ k = fl((1/k!)Ak) なので，^^

| Ck | ≦ | Bk | + | Bk – Bk | = (1/k!)| Ak | + | Bk – Bk |として前ページの結果を用いる．

11

全体の誤差• 全体の丸め誤差
級数の和の誤差1/k! 倍の誤差
Ak の計算の誤差• 全体の誤差12混合精度計算の場合• 計算法– exp(A) のテイラー展開では，高次の項ほど寄与が小さくなるため，この部分を低精度で計算しても，全体の精度が落ちないと期待される．
– A^k = fl(Ak) の計算において，k≦k1 のときは高精度で計算し，k > k1 のときは（高精度で計算した Ak1 に基づき）低精度で計算を行う．

• Ak の誤差– 低精度演算での丸め誤差の単位を us とし，γns = nus / (1 – nus) とおく．
– このとき，A^k （k > k1）の丸め誤差の上界は次のように与えられる．

13

混合精度計算の場合（続き）• 全体の丸め誤差2)
高精度計算部分の誤差

低精度計算部分の誤差• 切り替え点 k1 の決定– （高精度部分の最大誤差）≒（低精度部分の最大誤差）となるようにk1 を決めると，||A||F ≦ 1 の場合，,– すなわち，
2) Y. Yamamoto, S. Kudo and T. Hoshi: Error analysis of the truncated Taylor series expansion method forcomputing matrix exponential, JSIAM Lett, to appear.

14確率的丸め誤差解析• 確定的丸め誤差解析の問題点– 最悪ケースの解析• 計算過程での丸め誤差が，すべて最大値を取り，かつ，全体の誤差を増大させる方向の符号を持つと仮定．

– 長さ N のベクトルの内積の丸め誤差上界は O(Nu)．
• 実際の数値結果で多く見られる O(N1/2u) の振る舞いと乖離．

– N が非常に大きい場合，あるいは低精度計算の場合は，意味のある誤差上界が得られないことがある．

• 確率的丸め誤差解析– 丸め誤差を乱数として扱い，確率的な打ち消し合いを考慮．
– 長さ N のベクトルの内積の丸め誤差上界は O(N1/2u)．
– 上界は，ある（1に近い）確率で成り立つ確率的上界として得られる．
15

新しい確率的丸め誤差解析手法• 従来の確率的丸め誤差解析– 丸め誤差に一様分布や正規分布などの簡単な分布を仮定．
• 実際の丸め誤差の分布を必ずしも反映していない．

– O(u2) の項を無視して成り立つ解析– 中心極限定理に基づくため，N が大きい場合にのみ適用可能．

• 新しい確率的丸め誤差解析1)– 丸め誤差に特定の分布を仮定せず，期待値 0 と独立性のみを仮定．
– 高次の項についても厳密な解析– 集中不等式3)に基づくため，任意の N について適用可能．

1) N. J. Higham & T. Mary: A New Approach to Probabilistic Rounding Error Analysis, SIAM J. Sci. Comput.,Vol. 41, No. 5, pp. A2815-A2835 (2019).

3) 確率変数の和がその期待値からずれる確率の上界を，ずれの関数として与える不等式．
16

新しい確率的丸め誤差解析手法（1）• 定義（確率的誤差解析での γn）• 定理1（[1], Theorem 2.4）
– δ1, δ2, …, δn を平均 0 で絶対値が u 以下の独立な乱数とし，ρi = ±1（1≦ i ≦ n）とする．このとき，任意の λ > 0 に対して，次の式が次の P(λ) 以上の確率で成り立つ．

u
誤差解析で頻繁に現れる式.の確率的上界を与える．
17

新しい確率的丸め誤差解析手法（2）• 定義• 意味と利用法– 各々が確率 P(λ) で起こる独立な事象 n 個が同時に起こる確率の下界–の形の式が複数個所で現れる場合に，その全てが誤差上界以下となる確率を計算するのに使われる．

– 性質： 1 – {(1 – Q(λ, m)) + (1 – Q(λ, n))} = Q(λ, m+n)P(λ)

確率 1 – P(λ) 以下で起こる事象（誤差が上界を超える）

Q(λ, 3)の形の式が3カ所で現れる場合
18

新しい確率的丸め誤差解析手法（3）• 定理2（[1], Theorem 3.1 とその後の注）– a, b ∈ Rn とし，y = aTb とする．このとき，|fl(y) – y| ≦|a|T|b| が少なくとも確率 Q(λ, n) で成り立つ．
確定的誤差解析の γ をγ~n(λ) で置き換えた式n• 略証– 通常の誤差解析と同様，|εi| ≦ u, |δj| ≦ u に対し，次の式が成り立つ．

fl(y)– ここで，εi, δj を独立な乱数と見ると，定理1より，各項について確率P(λ) 以上で次の不等式が成り立つ．

– 項は n 個あるから，そのすべてが上界を満たす確率は Q(λ, n) 以上.
19

新しい確率的丸め誤差解析手法（4）• 定理3（[1], Theorem 3.4）
– A∈ Rm×n, A∈ Rn×p とし，C = AB とする．このとき，|fl(C) – C| ≦||A| |B| が少なくとも確率 Q(λ, mnp) で成り立つ．

• 証明– 行列乗算 C = AB は長さ n の内積 mp 個からなり，その全てで定理2の上界が成り立たなくてはならない．その確率は Q(λ, mnp) 以上．

20

exp(A) の計算の確率的誤差解析• 補題4（Ak の誤差）
– Q(λ, (k − 1)N3) 以上の確率で次の不等式が成り立つ．

≒• 補題4の証明– ある正数 ε~k–1 > 0 について，Q(λ, (k − 2)N3) 以上の確率で次の式が成り立つと仮定する．

– この式が成り立つとき，次の式も成り立つ．

21

補題4の証明 （続き）^

• さらに，定理3より，Ak−1 と A の乗算において，次の不等式が Q(λ, N3) 以上の確率で成り立つ．

• これらより，以上の確率で，次の不等式が成り立つ，

これを ~εk と置ける

22

補題4の証明 （続き）
• ε~k に関する漸化式– 前ページの結果より，．
– ε~1 = 0 に注意してこれを解くと，ε~ k に対する次の上界が得られる．
~γ の形にするために，ここを大きくした．

– すなわち，．
23

exp(A) の計算の確率的誤差解析• 1/k! 倍の誤差と級数の和の誤差– 1/k! 倍の誤差と級数の和の誤差についても考慮し，確率をアップデートすると，次のように全体の丸め誤差の確率的評価が得られる．

• 定理5–以上の確率で次の不等式が成り立つ．

注意
解析過程からわかるように，ある Ak が誤差上界を満たすという事象は，Al （k ≦ l）が誤差上界を満たすという事象に含まれるので， Ak の誤差の確率については，k = n のときの確率のみを考えればよい．
24

混合精度計算の場合• 計算法– A^k = fl(Ak) の計算において，k≦k1 のときは高精度で計算し，k > k1 のときは（高精度で計算した Ak1 に基づき）低精度で計算を行う．

• 定理6（結果のみ）–以上の確率で次の不等式が成り立つ．

高精度計算部分の誤差
低精度計算部分の誤差.
25確率的誤差解析のまとめ• ポイント(1) まず，確率を付与すべき事象（アルゴリズム中のある行列積の誤差が一定値以下となるなど）を定義し，その包含関係を整理する．

(2) 通常の誤差解析と同様にして，計算の各ステップで生じる誤差の上界を求める．ただし，γn を ~γn(λ) で置き換える．
(3) 誤差上界を新しく作るごとに，その成立確率についてもアップデートする．
(4) 複数の γn(λ) が現れる式を簡単化するには，次の公式を使う．
・通常の誤差解析におけるに相当．

(5) 誤差が求めた上界を超える確率が十分小さな値（10–10以下など）になるよう l の値を調節する．
26

数値実験
• 目的– テイラー展開法で exp(A) を計算し，実際の誤差と誤差上界を比較• 実験条件– テスト行列： [-0.5, 0.5] の一様乱数で生成し，スケーリング– 計算精度： 単精度（一部は半精度）– 実際の誤差： 倍精度の固有値分解法での結果を真値として算出– 確率的誤差上界については，誤差が上界を超える確率が 10–10 以下になるように，λ = 10 と設定• 計算機環境– Intel Xeon E2660 V2 processor （2.2MHz, 10コア）– Intel Fortran compiler– Intel Math Kernel Library （LAPACK/BLAS）
27

数値結果（1）• A のノルムを変えた場合– N = 1000, 2000 とし，|| A ||F を 10–3 から 1 まで変えて，実際の誤差，確定的誤差上界，確率的誤差上界を比較．
– テイラー展開の次数 n は，打切り誤差が 10–7 以下になるよう設定．
|| A ||F|| A ||F確定的誤差上界
誤差

確率的誤差上界
実際の誤差

n=1

n=2

n=4

N=1000n=9n=2

n=3

n=5

n=10

N=2000– 確率的誤差上界の方が実際の誤差とのギャップが小さい．
– 実際の誤差とのギャップは，ノルムが大きくなるにつれて増大

28

数値結果（2）• A のサイズを変えた場合– || A ||F = 1 と固定し，N を 200 から 3200 まで変えて，実際の誤差，確定的誤差上界，確率的誤差上界を比較．
N

確定的誤差上界

誤差

確率的誤差上界
実際の誤差– 確定的誤差上界は N に比例して増加– 確率的誤差上界は N1/2 に比例し，実際の誤差により近い結果
29数値結果（3）• 半精度演算による結果– N = 400, || A ||F = 1 とし，k > k1 のときの Ak の計算を半精度で実行– 半精度演算は，単精度演算の結果の仮数部を 11 ビットに切り捨ててエミュレート– 次数を k = k1 で打ち切った場合の結果も同時にプロット
N

単精度

誤差

k > k1 のとき半精度
k = k1 で打ち切り– 4次以上の項を半精度で計算しても，単精度とほぼ同じ精度を達成– しかし，完全な打ち切りと比べた精度向上はあまり大きくない
30

事後誤差解析への応用（1）• 事前誤差解析– 入力データから計算される量（ノルム，条件数など）のみを使って丸め誤差上界を導くタイプの誤差解析– 計算を行う前に丸め誤差を推定できる．
– これまで説明した誤差解析事例はすべて事前誤差解析• 事後誤差解析– 浮動小数点による計算結果を誤差上界に含むタイプの誤差解析– 計算後（or 計算途中）でないと誤差上界の値が求められない．

– 事前誤差解析よりもタイトな誤差上界が得られる場合が多い．
– 事前誤差解析が困難でも事後誤差解析は可能な場合がある．
^^• LU 分解の後退誤差 A + E = LU など
31

事後誤差解析への応用（2）• ブロッククリロフ部分空間法– 複数右辺を持つ方程式 AX = B （A∈Rn×n, B, X∈Rn×s）を解くためのクリロフ部分空間法– ブロッククリロフ部分空間を用いて解を求める．

• 残差ギャップの問題^^– ブロッククリロフ部分空間法では，残差ベクトル Rk を近似解 Xk から^^
Rk = B – AXk と計算するのではなく，漸化式により計算する．
^^– そのため，アルゴリズム中の残差 Rk と真の残差 B – AXk との間に食い違いが生じることがある．これを残差ギャップと呼ぶ．
– 残差ギャップは偽収束等の問題をもたらすため，その評価は重要4)
確率的な事後誤差解析を用いた残差ギャップの上界の導出
4) K. Aihara, A. Imakura and K. Morikuni: Cross-interactive residual smoothing for global and llock Lanczostype solvers for linear systems with multiple right-hand sides, SIMAX, Vol. 43, No. 3 (2022).

32

残差ギャップの上界の導出（1）• 解の更新式と残差の更新式– 列直交なブロック進行方向ベクトル Qk–1∈Rn×s と係数 ak–1∈Rs×s を用いてXk と Rk を次のように更新する．
Xk = Xk–1 + Qk–1ak–1
Rk = Rk–1 – (AQk–1)ak–1•ak–1 のノルムの評価– 漸化式より， Qk–1ak–1 = Xk – Xk–1.
– Qk–1が列直交行列であることより，事後評価
33

残差ギャップの上界の導出（2）^• Xk の確率的誤差評価– 次式で定義される誤差 E1, E2 を評価する．

行列乗算の誤差

加算の誤差– 定理3（行列乗算の誤差）より，確率 Q(λ, ns2) 以上で次が成り立つ．

∴
確率的事後誤差上界– 一方，確定的事後誤差上界– これらより，確率 Q(λ, ns2) 以上で，34

残差ギャップの上界の導出（3）^• Rk の確率的誤差評価– 次式で定義される誤差 Ek(3), Ek(4), Ek(5) を評価する．

– まず，行列 A の1行当たりの非ゼロ要素数の最大値を m とすると，（確率≧ Q(λ, nms)）∴– 次に，u2の項を無視すると，次の評価が成り立つ．

∴（確率≧ Q(λ, ns2)）
35

残差ギャップの上界の導出（4）^• Rk の確率的誤差評価（続き）– また，– 以上より，^^• Xk と Rk の表式（確率≧ Q(λ, n(m+s)s)）36

残差ギャップの上界の導出（5）• 残差ギャップの確率的上界（確率≧ Q(λ, n(m+2s)sk)）• 定理7– ブロッククリロフ部分空間法において，解 X^ k と残差 R^ k を列直交なブロック方向ベクトル Q^ k を用いて更新するとき，残差ギャップ (B – AX^ k)^– Rk は Q(λ, n(m+2s)sk) 以上の確率で上記の上界を満たす．
真の残差の推定や偽収束の検出への応用
37

おわりに

• 打切り型テイラー展開に基づく行列指数関数の計算法を例題として，
Higham らによって提案された確率的誤差解析の実際を解説した．
• 確率的誤差解析では，N1/2（N: 行列の次元数）に比例する上界が得られた．これは，実際の誤差の振る舞いと一致している．

• テイラー展開の高次の項を半精度で計算する混合精度計算，及びブロッククリロフ部分空間法における残差ギャップについても，誤差上界を導出した．

38