---
title: 第8回 配信講義　計算科学技術特論A（2025）
tags: 
author: [R-CCS　計算科学研究推進室](https://convert.docswell.com/user/2300203199)
site: [Docswell](https://www.docswell.com/)
thumbnail: https://bcdn.docswell.com/page/YE9XNNVP73.jpg?width=480
description: 第８回  6 月 5日  行列計算における高速アルゴリズム１ 本講義では、科学技術計算の基盤となる行列計算について、基本的なアルゴリズムの原理と高性能計算技術を紹介する。 第1回ではクリロフ部分空間に基づく連立1次方程式反復解法と固有値計算法の基礎について解説する。
published: May 26, 25
canonical: https://convert.docswell.com/s/2300203199/KQRWLN-2025-05-26-105928
---
# Page. 1

![Page Image](https://bcdn.docswell.com/page/YE9XNNVP73.jpg)

行列計算における高速アルゴリズム
- 大規模連立1次方程式の反復解法 2025年 6月 5日
計算科学技術特論A
電気通信大学 大学院情報理工学研究科
情報・ネットワーク工学専攻
山本 有作


# Page. 2

![Page Image](https://bcdn.docswell.com/page/GE82KKPXED.jpg)

本講義の構成
• 第1回（6月5日）
「大規模連立1次方程式の反復解法」
– 代表的な行列計算である連立1次方程式の解法について，現在
最も広く利用されているクリロフ部分空間法の基礎を紹介する．
• 第2回（6月12日）
「ポストペタ時代に向けた線形計算アルゴリズムの課題と
研究動向」
– 様々な行列計算について，ポストペタ計算機上での実装における
課題と，その解決のための取り組みを紹介する．
– また，行列計算における量子アニーリングマシンの活用の可能性
についても簡単に触れる．
（1/61）


# Page. 3

![Page Image](https://bcdn.docswell.com/page/LEL2ZZ1N7R.jpg)

目次
1. はじめに
2. クリロフ部分空間と射影
3. 代表的なクリロフ部分空間法
4. 前処理
5. 波束シミュレーション向けの解法
6. 終わりに
（2/61）


# Page. 4

![Page Image](https://bcdn.docswell.com/page/4JM8VVWQJW.jpg)

はじめに
（3/61）


# Page. 5

![Page Image](https://bcdn.docswell.com/page/PJR522YKE9.jpg)

大規模連立1次方程式
• 今回の講義で扱う問題
・ 連立1次方程式
Ax = b
・ ただし，A∈Rm×m， A は正則，
b∈Rm （右辺ベクトル）， x∈Rm （未知ベクトル）
・ 特に A が大規模な疎行列の場合を考える
• 応用
・ 差分法・有限要素法による偏微分方程式の解法
・ 連立常微分方程式の解法（回路シミュレーションなど）
・ 連立非線形方程式の解法
・ 他の線形計算の部品（shift-and invert Lanczos 法，櫻井-杉浦法）
（4/61）


# Page. 6

![Page Image](https://bcdn.docswell.com/page/PEXKDDG5JX.jpg)

係数行列の分類と数値解法
• 係数行列の分類
・ 密行列／疎行列／帯行列／構造を持つ行列
・ 対称行列／非対称行列
・ 正定値行列／その他の行列
• 行列の分類に応じた解法（例）
・ 幅の狭い帯行列ならば直接解法（LU分解法）が有利
・ 対称正定値行列ならばコレスキー分解が利用可能
・ 対称正定値行列ならば共役勾配法（CG法）が利用可能
参考文献：
・ 櫻井鉄也, 松尾宇泰, 片桐孝洋編： 「数値線形代数の数理とHPC」, 共立出版, 2018.
・ 張紹良編著： 「計算科学のための基本数理アルゴリズム」, 共立出版, 2019.
（5/61）


# Page. 7

![Page Image](https://bcdn.docswell.com/page/3EK5XX4RJD.jpg)

直接解法と反復解法
• 直接解法
・ 行列を A = LU とLU分解し，Ly = b，Ux = y を順に解く
・ A が疎行列の場合，ゼロでない要素のみに対して演算を行うことで，
演算量と記憶領域を削減可能（疎行列LU分解）
・ A が対称正定値行列の場合，コレスキー分解 A = LLT を用いることで，
記憶領域と演算量を半分にできる
• 特徴
・ 有限回の演算で解が求まる（所要時間がほぼ予測可能）
・ 悪条件の問題に強い
・ 行列の非零要素数に比べて、数倍～数十倍の記憶容量が必要なこと
が多い（フィルインの問題）
・ 反復解法に比べると，演算量大
・ 1回LU分解を行えば， b のみが異なる方程式を少ないコストで解ける
（6/61）


# Page. 8

![Page Image](https://bcdn.docswell.com/page/L73K55LG75.jpg)

直接解法と反復解法（続き）
• 反復解法
・ 真の解 x* に収束するベクトル列 x(0), x(1), x(2), … を順次生成していく
方法
• 特徴
・ 反復回数が事前にはわからない（ことが多い）
・ 収束しない場合がある
・ 行列 A を変形せず，行列ベクトル積 Ap の形でのみ用いるため，行列
の記憶領域は A の非ゼロ要素数の分のみ
# ただしベクトルの記憶領域は別に必要
・ うまく収束する場合，演算量は直接解法に比べて遥かに小さい
（7/61）


# Page. 9

![Page Image](https://bcdn.docswell.com/page/87D3DDNNJG.jpg)

定常反復法と非定常反復法
• 定常反復法
・ 線形の漸化式 x(n+1) = Cx(n) + d （C：定数行列，d：定数ベクトル）により
近似解を更新していく反復解法
・ 上記の反復が収束し，収束先が真の解 x* になるように C，d を選ぶ
・ 収束次数は1次。収束率は C の絶対値最大の固有値に依存
・ ガウス-ザイデル法，SOR法，ADI法，マルチグリッド法など
• 非定常反復法
・ x(n) → x(n+1) の漸化式が，x(n) に関して非線形・非定常な反復解法
・ 一般に，定常反復法よりも収束が速い
・ クリロフ部分空間法，チェビシェフ準反復法など
（8/61）


# Page. 10

![Page Image](https://bcdn.docswell.com/page/VJP4XX5N78.jpg)

行列の特異値と条件数
• 行列の特異値
・ A∈Rm×m に対し，ATA は対称非負定値行列であり，その m 個の固有
値はすべて 0 以上の実数となる
・ これらの平方根を大きい順に 1(A) ≧ 2(A) ≧ ・・・ ≧ m(A) ≧0と
書き，Aの特異値と呼ぶ
・ 行列 A が正則 ⇔ m(A) &gt; 0
• 行列の条件数
・ 正則行列 A∈Rm×m に対し，(A) = 1(A)/m(A) を A の条件数と呼ぶ
・ 単位行列の条件数は 1
・ 条件数の大きな連立1次方程式は，次の2つの意味で解きにくい
～
# 残差 r(n) = Ax(n) – b が小さくても，x(n) が真の解から遠い場合がある
# クリロフ部分空間法が収束しにくい
（9/61）


# Page. 11

![Page Image](https://bcdn.docswell.com/page/2EVXLLPYJQ.jpg)

クリロフ部分空間と射影
（10/61）


# Page. 12

![Page Image](https://bcdn.docswell.com/page/57GVXXPW7L.jpg)

クリロフ部分空間
• 定義
・ 行列 A∈Rm×m，ベクトル b∈Rm に対し，次のように定義される Rm の
部分空間 Kn(A; b) を A，b に関する n 次のクリロフ部分空間と呼ぶ
Kn(A; b) = span{b, Ab, A2b, …, An–1b}
• クリロフ部分空間法
・ x(0) = 0 から始めて，x(n)∈Kn(A; b) （n ≧1）となるように近似解の列
{x(n)} を構成していく反復解法を，クリロフ部分空間法と呼ぶ
# 0 以外の初期値から始める方法もある
・ K1(A; b) ⊆ K2(A; b) ⊆ ・・・ ⊆ Kn(A; b) なので，クリロフ部分空間法では，
探索する空間を広げながら近似解を更新していくことになる
・ Kn(A; b) の中で，どのような基準で近似解 x(n) を定めるかにより，様々
なクリロフ部分空間法がある
（11/61）


# Page. 13

![Page Image](https://bcdn.docswell.com/page/4EQ6882QEP.jpg)

クリロフ部分空間内での解の存在
• クリロフ部分空間の拡大
– Kk(A; b) ⊆ Kk+1(A; b) は明らか
– あるところで拡大が止まって Kk(A; b) = Kk+1(A; b) となるはず．
• 疑問
– 初めて Kk(A; b) = Kk+1(A; b) となるのはどんなときか？
– 拡大が止まることで，真の解 x* に到達できない可能性は？
KK44 = K5
K3
x3
K2
x2
K1
x1
x4
x*
x0 = 0
（12/61）


# Page. 14

![Page Image](https://bcdn.docswell.com/page/KJ4422LYJ1.jpg)

クリロフ部分空間内での解の存在（続き）
• 定理
– A を正則行列とする。このとき，Kk(A; b) = Kk+1(A; b) を満たす最小の
k に対し，Ax = b の解 x* は Kk(A; b) に含まれる．
• 証明
– 仮定より， Akb = c0b + c1Ab + ・・・ + ck–1Ak–1b を満たす c1, c2, …, ck–1
が存在する。ここで，c0 = 0 とすると，A の正則性より，k が最小である
ことと矛盾が生じるので，c0≠0。そこで，両辺を c0 で割って整理すると，
A((1/c0)Ak–1b – (ck–1/c0)Ak–2b – ・・・ – (c1/c0)b) = b。これは x*∈Kk(A; b)
を意味する．
部分空間の拡大が止まって真の解 x* に到達できないという心配は不要．
（13/61）


# Page. 15

![Page Image](https://bcdn.docswell.com/page/LE14MMLNEG.jpg)

残差多項式
• 残差の表式
・ Kn(A; b) における近似解を x(n) = d0b + d1Ab + ・・・ + dn–1An–1b と定め
たとする
・ このとき，残差 r(n) = Ax(n) – b は次のように書ける
r(n) = – b + d0Ab + d1A2b + ・・・ + dn–1Anb = n(A)b
ただし，
n(z) = – 1 + d0z + d1z2 + ・・・ + dn–1zn 定数項が –1の n 次多項式
これを残差多項式と呼ぶ
・ Kn(A; b) で x(n) を1つ定めるとは，残差多項式を1つ定めることと等価
・ 残差多項式は，クリロフ部分空間法の解析で重要な役割を果たす
（14/61）


# Page. 16

![Page Image](https://bcdn.docswell.com/page/GEWX33QM72.jpg)

クリロフ部分空間の正規直交基底の生成
• 正規直交基底の必要性
・ 近似解 x(n) を Kn(A; b) の基底ベクトルの線形結合として表す
・ 基底として正規直交基底を用いると，解法の理論的導出のし易さでも，
数値的安定性の面からも有利
• 正規直交基底の生成
・ 考え方
・ b を正規化して q1 = b /∥b∥2 を作成
・ v2 = Aq1
・ v2 を q1 に対して直交化し，正規化して q2 を作成
・ v3 = Aq2
・ v3 を q1, q2 に対して直交化し，正規化して q3 を作成
・ これは，逐次的に生成されるベクトル {vn} に対し，グラム・シュミット
法による直交化を行っていることに相当
・ これを Arnordi 過程と呼ぶ
（15/61）


# Page. 17

![Page Image](https://bcdn.docswell.com/page/47Z6VVWM73.jpg)

クリロフ部分空間の正規直交基底の生成（続き）
• Arnoldi 過程
新たなベクトルの生成（空間の拡張）
直交化
正規化
（16/61）


# Page. 18

![Page Image](https://bcdn.docswell.com/page/YJ6255G57V.jpg)

クリロフ部分空間の正規直交基底の生成（続き）
• Arnoldi 分解
・ Arnoldi 過程 より，Aqi は q1, q2, …, qi+1 を使って次のように書ける
Aqi = h1i q1 + h2i q2 + ・・・ + hi+1,i qi+1
・ そこで，
（m×n 列直交行列）
（(n+1)×n 行列）
と定義すると，次の式が成り立つ
～
AQn = Qn+1Hn
・ これを，A の n 次の Arnoldi 分解と呼ぶ
（17/61）


# Page. 19

![Page Image](https://bcdn.docswell.com/page/GJ52YYVGJ4.jpg)

クリロフ部分空間の正規直交基底の生成（続き）
• Arnoldi 分解（続き）
～
・ Arnoldi 分解の式 AQn = Qn+1Hn に左から QnT をかけ，Qn の各列が直
交することに注意すると，
（n×n ヘッセンベルグ行列）
・ これは，A を span{q1, q2, …, qn} = Kn(A; b) に射影した行列と見ること
ができる
# 固有値計算のための Arnoldi 法の原理
（18/61）


# Page. 20

![Page Image](https://bcdn.docswell.com/page/LE3K55DPE5.jpg)

固有値計算のためのArnoldi法
• クリロフ部分空間の中での固有ベクトルの近似
– Kn(A; b) の中で A の近似固有ベクトル x を求めることとし，x = Qy
（y∈Rn）とおく．また，対応する近似固有値を  とする
– 残差 Ax–x が Kn(A; x(0)) に直交するという条件（Ritz-Galerkin 条件）
を課すと，
QnT(Ax – x) = QnTAQny – QnTQn
= Hny – y = 0.
したがって，k×k 行列に対する小型の固有値問題
H n y = y
を解けばよい
これは，ヘッセンベルグ行列に対する固有値問題となる
– 近似固有値は ，近似固有ベクトルは x = Qny で与えられる
（19/61）


# Page. 21

![Page Image](https://bcdn.docswell.com/page/8ED3DDP37G.jpg)

A が対称行列の場合
• Lanczos 過程
・ QnTAQn = Hn は対称行列となるから，非ゼロ構造を考えると，
となり，Hn は対称3重対角行列 Tn となる
・ これは，Anoldi 過程において計算した内積 hji （i ≧ j+1）が自動的に 0
になることを意味する
・ すなわち，Aqi は qi と qi–1 に対してのみ直交化すればよい
・ これを利用して Arnoldi 過程の計算量を減らした方法を Lanczos 過
程と呼ぶ
（20/61）


# Page. 22

![Page Image](https://bcdn.docswell.com/page/V7P4XXNPE8.jpg)

A が対称行列の場合（続き）
• Lanczos 過程
新たなベクトルの生成（空間の拡張）
直交化
正規化
・ ただし，n = hnn, n = hn+1,n = hn,n+1 とおいた
• Lanczos 過程の特徴
・ qn と qn–1 のみから qn+1 が求まる（3項間漸化式）
・ 1反復あたりの演算量は n によらず一定
（21/61）


# Page. 23

![Page Image](https://bcdn.docswell.com/page/2JVXLLKV7Q.jpg)

固有値計算のためのLanczos法
• クリロフ部分空間の中での固有ベクトルの近似
– Kn(A; b) の中で A の近似固有ベクトル x を求めることとし，x = Qy
（y∈Rn）とおく．また，対応する近似固有値を  とする
– 残差 Ax–x が Kn(A; x(0)) に直交するという条件（Ritz-Galerkin 条件）
を課すと，
QnT(Ax – x) = QnTAQny – QnTQn
= Tny – y = 0.
したがって，k×k 行列に対する小型の固有値問題
Tn y =  y
を解けばよい
これは，3重対角行列に対する固有値問題となる
– 近似固有値は ，近似固有ベクトルは x = Qny で与えられる
（22/61）


# Page. 24

![Page Image](https://bcdn.docswell.com/page/5EGVXXQ1EL.jpg)

クリロフ部分空間内での近似解の決定
• 近似解の選び方
・ 近似解 x(n) としてクリロフ部分空間 Kn(A; b) 内のどんなベクトルを選
ぶかは，自由度がある
・ 主に，残差最小化， Ritz-Galerkin 法，Petrov-Galerkin 法の3種類
の原理が使われる
・ どの原理を使うかによって，クリロフ部分空間法の様々な変種が生ま
れる
（23/61）


# Page. 25

![Page Image](https://bcdn.docswell.com/page/4JQ688KNJP.jpg)

近似解の決定法 I： 残差最小化
• 方法
・ 残差 r(n) = Ax(n) – b が最小になる x(n) を Kn(A; b) の中で選ぶ
・ 最も自然なアプローチ
# 本当は誤差 e(n) = x(n) – x* を最小にしたいが，一般には困難
・ GMRES法，MINRES法などがこのタイプに属する
（24/61）


# Page. 26

![Page Image](https://bcdn.docswell.com/page/K74422P371.jpg)

近似解の決定法 II： Ritz-Galerkin法
• 方法
・ 残差 r(n) = Ax(n) – b が Kn(A; b) に直交するように x(n) を選ぶ
Ax(n) – b ⊥ Kn(A; b)
・ これは，残差のうち Kn(A; b) に正射影した成分が 0 であることを示す
• 残差の直交性
・ K1(A; b) ⊂ K2(A; b) ⊂ ・・・ ⊂ Kn(A; b) （真部分空間）とするとき，残差
r(0) = b, r(1), r(2), …, r(n–1) は Kn(A; b) の直交系をなす
# r(n–1)∈Kn(A; b) を用いて，n に関する帰納法により容易に示せる
（25/61）


# Page. 27

![Page Image](https://bcdn.docswell.com/page/LJ14MM9ZJG.jpg)

近似解の決定法 II： Ritz-Galerkin法（続き）
• A が対称正定値行列の場合の意味付け
・  (x) = (1/2) xTAx – xTb とおくと， Ax = b の解は  (x) の最小値
・ ここで， x(n)∈Kn(A; b) と制限すると，x(n) = Qny(n) （y(n)∈Rn）と基底の線
形結合で書けるので， (x(n)) = (1/2) (Qny(n))T AQny(n) – (Qny(n))Tb
・ これを最小にする y(n) は，微分して QnT(AQny(n) – b) = 0 より求まる
・ これは，ちょうど Ax(n) – b ⊥ Kn(A; b) と等価
・ すなわち，Ritz-Galerkin法はクリロフ部分空間 Kn(A; b)の中で  (x(n))
を最小化する x(n) を求めている
（26/61）


# Page. 28

![Page Image](https://bcdn.docswell.com/page/GJWX33Q6E2.jpg)

近似解の決定法 II： Ritz-Galerkin法（続き）
• A が対称正定値行列の場合の意味付け（続き）
・ いま，誤差 e(n) = x(n) – x* の A ノルムを考えると，
・ したがって，Ritz-Galerkin法は，誤差の A ノルムを最小化する x(n) を
求めていると考えることもできる
• Ritz-Galerkin 法に基づく解法
・ 対称正定値行列向け： CG法（共役勾配法）
・ 非対称行列向け： FOM法
# Ritz-Galerkin 法は関数空間での近似（有限要素法など）にも使われる
（27/61）


# Page. 29

![Page Image](https://bcdn.docswell.com/page/4EZ6VVWRE3.jpg)

近似解の決定法 III： Petrov-Galerkin法
• 方法
・ Rm の適当な部分空間の列 L1⊂L2⊂L3⊂ ・・・ （Ln は n 次元空間）を用
意し，残差 r(n) = Ax(n) – b が Ln に直交するように x(n) を選ぶ
Ax(n) – b ⊥ Ln
・ これは，残差のうち空間 Ln に正射影した成分が 0 であることを示す
・ Ln としては，b* を適当なベクトルとし，AT と b* に対するクリロフ部分
空間 Kn(AT; b*) を使うことが多い
• Petrov-Galerkin 法に基づく解法
・ Bi-CG法，QMR法
# Petrov-Galerkin 法は関数空間での近似にも使われる
・ Petrov-Galerkin 法から派生した解法： CGS法，Bi-CGSTAB法，
GPBi-CG法
（28/61）


# Page. 30

![Page Image](https://bcdn.docswell.com/page/Y76255G1EV.jpg)

残差多項式による特徴付け
• 残差最小化
・ 定数項が –1 の n 次多項式の集合を Pn とする
・ このとき，残差最小化は，各 n における残差多項式 n(z) を，次の最
小化問題の解になるように選ぶアプローチだと解釈できる
minn ∈Pn ∥n(A)b∥2
• Ritz-Galerkin法（Aが対称正定値行列の場合）
・ x(n) の誤差を e(n) = x(n) – x* とすると，Ae(n) = Ax(n) – b = r(n)
・ よって，
e(n)TAe(n) = r(n)TA–1r(n) = bTn(A)A–1n(A)b
= bTA–1n(A)An(A)A–1b
= e(0)Tn(A)An(A)e(0) = ∥n(A)e(0)∥A
・ すなわち，Ritz-Galerkin 法は次の最小化問題を解いている
minn ∈Pn ∥n(A)e(0)∥A
（29/61）


# Page. 31

![Page Image](https://bcdn.docswell.com/page/G752YYVL74.jpg)

代表的なクリロフ部分空間法
（30/61）


# Page. 32

![Page Image](https://bcdn.docswell.com/page/9J24GG337R.jpg)

代表的なクリロフ部分空間法
• 残差最小化に基づく解法
・ GMRES法（Generalized Minimum Residual，一般化最小残差法）
・ MINRES法（Minimum Residual，最小残差法）
• Ritz-Galerkin 法に基づく解法
・ CG法（Conjugate Gradient，共役勾配法）
・ FOM法（Full Orthogonalizaion Method，完全直交化法）
• Petrov-Galerkin 法に基づく解法
・ Bi-CG法（Bi-Conjugate Gradient，双共役勾配法）
・ QMR法（Quasi Minimal Residual，準最小残差法）
・ CGS法（Conjugate Gradient Squared，自乗共役勾配法）
・ Bi-CGSTAB法（Stabilized Bi-CG，安定化双共役勾配法）
（31/61）


# Page. 33

![Page Image](https://bcdn.docswell.com/page/DEYMVV18EM.jpg)

GMRES法
• 原理
・ Kn(A; b) の中で残差 r(n) = Ax(n) – b を最小にする x(n) を 選ぶ
・ x(n)∈Kn(A; b) より x(n) = Qny(n) （y(n)∈Rn）と書くと，
～
・ ただし，第2の等号では Arnoldi 分解の式 AQn = Qn+1Hn を用いた
・ 第3の等号では，q1 = b /∥b∥2 であることを用いた
・ 第4の等号では，Qn+1 が列直交行列であることから，Qn+1 を除いても
2ノルムが変わらないことを用いた
・ 最後の式の∥・∥2 の中は n+1 次元ベクトル，y(n) は n 次元ベクトルであ
るから，残差の最小化は y(n) に関する最小2乗問題となる
（32/61）


# Page. 34

![Page Image](https://bcdn.docswell.com/page/VJNWMM29J8.jpg)

GMRES法（続き）
• 最小2乗問題の解法
～
～
・ 最小2乗問題 min∥Hny(n) –∥b∥2 e1∥2 を解くには，係数行列 Hn を Hn=
～
Un Rn （Un∈R(n+1)×n, Rn∈Rn×n）とQR分解し，連立1次方程式 Rny(n) =
UnT∥b∥2 e1 を解けばよい
～
・ Hn は (n+1)×n ヘッセンベルグ行列なので，そのQR分解はギブンス
回転を用いて O(n2) の計算量で行える
～
～
・ さらに，Hn–1 のQR分解が与えられたとき，Hn のQR分解は O(n) で計
算できる
# 追加された第 n 列のみについて，ギブンス回転を作用させればよい
（33/61）


# Page. 35

![Page Image](https://bcdn.docswell.com/page/YE9XNN1373.jpg)

GMRES法（続き）
• アルゴリズム
Arnoldi 過程による Kn(A; b) の
正規直交基底の生成
～
ギブンス回転による Hn のQR分解
～
（Hn のQR分解の結果を用いる）
UnT∥b∥2 e1 の計算
後退代入により
Rny(n) = UnT∥b∥2 e1 を解く
x(n) = Qny(n) の計算
（34/61）


# Page. 36

![Page Image](https://bcdn.docswell.com/page/GE82KK4LED.jpg)

GMRES法（続き）
• GMRES法の特徴
・ 計算の途中で破綻が起こることはない
・ 残差が n とともに単調に減少
# 非対称行列向けのクリロフ部分空間法としては，最も頑健な方法
・ 所要メモリ量と1反復あたりの演算量が n に比例して増加
# qn を q1, q1, …, qn–1 に対して直交化する必要があるため
# 所要メモリ・演算量を削減する方法としてリスタートがある
• MINRES法
・ 対称行列（必ずしも正定値でない）のための残差最小化法
～
・ Hn が3重対角行列になることを利用し，所要メモリ量と1反復あたりの
演算量を n に依存しない定数にできる
・ 丸め誤差に弱く，x(n) 中に ((A))2 のオーダーで丸め誤差が混入*
* H. A. van der Vorst： “Iterative Krylov Methods for Large Linear Systems”, Cambridge Univ. Press, 2003.
（35/61）


# Page. 37

![Page Image](https://bcdn.docswell.com/page/LEL2ZZ6Q7R.jpg)

CG法
• 原理
・ A は対称正定値行列とする
・ 残差 r(n) = Ax(n)–b ∈ Kn+1(A; b)が Kn(A; b) に直交するよう x(n) を選ぶ
# Ritz-Galerkin 法
• Lanczos過程との関係
・ Lanczos過程では，基底ベクトル qn+1∈ Kn+1(A; b) を，Kn(A; b) に直交
するよう選ぶ
・ したがって，Lanczos過程を実行し，qn+1 を適当にスケーリングすれば，
CG法の残差 r(n) が得られるはずである
・ 残差 r(n) の計算法が決まれば，解 x(n) の計算法も決まる（p. 14）
（36/61）


# Page. 38

![Page Image](https://bcdn.docswell.com/page/4JM8VVGKJW.jpg)

CG法（続き）
• Lanczos過程における基底ベクトルのスケーリング
・ qn+1 は3項漸化式を満たすから，それをスケーリングした r(n) も3項漸
化式を満たすはず
・ r(n) が r(n–1)，r(n–2) と直交するという条件より，
・ さらに，両辺の各項を {b, Ab, …, Anb} で展開したときの b の係数を
考えると，左辺は 0，右辺の r(n–1)，r(n–2)，r(n) はそれぞれ –1（p. 14）で
あるから，
・ これらにより r(n) の漸化式が定まる
（37/61）


# Page. 39

![Page Image](https://bcdn.docswell.com/page/PJR522Z6E9.jpg)

CG法（続き）
• CG法のアルゴリズム（3項漸化式）
– r(n) の漸化式に r(n) = Ax(n) – b を代入して x(n) の漸化式を導出
（38/61）


# Page. 40

![Page Image](https://bcdn.docswell.com/page/PEXKDD9DJX.jpg)

CG法（続き）
• アルゴリズム（連立2項漸化式）
ステップ幅の決定
探索方向に向かって xn を動かす
残差の計算
次の探索方向の決定
• 最小化問題の解法としての解釈
・ 関数  (x) = (1/2) xTAx – xTb を最小化
・ pn–1 が第 n ステップでの探索方向
・ pn–1 方向に進んで関数値が極小になるように，ステップ幅をn を決定
（39/61）


# Page. 41

![Page Image](https://bcdn.docswell.com/page/3EK5XXPDJD.jpg)

CG法（続き）
• CG法の性質
・ Kn(A; b) = span{b, Ab, A2b, …, An–1b} = span{x1, x2,…, xn}
= span{r0, r1,…, rn–1} = span{p0, p1,…, pn–1}
・ 残差の直交性： rnTrj = 0 (j &lt; n)
・ 探索方向の A-共役性： pnTApj = 0 (j &lt; n)
# 最小化問題の解法と見たとき，各方向 pn への探索の繰り返しにより，大
域的な最小値が求まることを保証
• CG法の特徴
・ 計算の途中で破綻が起こることはない
～
・ 誤差の A ノルム∥en∥A =∥xn – x*∥A が n とともに単調に減少
・ 所要メモリ量と1反復あたりの演算量は n に依存しない
# 3項間漸化式により計算が行えるため
（40/61）


# Page. 42

![Page Image](https://bcdn.docswell.com/page/L73K55XP75.jpg)

CG法（続き）
• CG法の収束性
・ 誤差の A ノルム∥en∥A = の減少について，次の評価式が成り立つ
ただし， は A の条件数
• 証明の概略
・ 残差多項式で表現した∥en∥A の最小性より，
～
ただし，Pn は定数項が1の n 次多項式，(A) は A の固有値の集合
・ ここで，特に  () をチェビシェフ多項式に取って最右辺を評価すると
上の定理が得られる
（41/61）


# Page. 43

![Page Image](https://bcdn.docswell.com/page/87D3DD23JG.jpg)

CG法（続き）
• CG法の収束性（続き）
・ A の相異なる固有値の個数が m* の場合，CG法は高々m* 回で収束
する
# 前ページの  () として，これら m* 個の固有値で 0 となる m* 次多項式
が取れるから，n = m* のとき，最右辺は 0 となる
・ 固有値が少数のクラスターに固まっている場合も，収束は速い
# これらのクラスター上で小さな値を取る低次の多項式が作れるから
• 前処理
・ 固有値を1個あるいは少数個のクラスターに集中させることができれ
ば，CG法の収束は速くなる
・ この目的のため，Ax = b を別の（対称正定値な係数行列を持つ）連
立1次方程式に変形して解く前処理（後述）が行われる
（42/61）


# Page. 44

![Page Image](https://bcdn.docswell.com/page/VJP4XX1P78.jpg)

Bi-CG法
• 原理
・ 方程式
に対して形式的にCG法を適用する
• アルゴリズム
（43/61）


# Page. 45

![Page Image](https://bcdn.docswell.com/page/2EVXLLYVJQ.jpg)

Bi-CG法（続き）
• Breakdown
・ rn と sn とが直交すると，アルゴリズム中で分母が0になり，計算を続
行できなくなる。これは，双直交基底を生成できないことに相当する。
これを第1種の breakdown と呼ぶ
・ 一方，Tn のLU分解において，ピボットが0になることに相当する
breakdown もある。これを第2種の breakdown と呼ぶ
・ これらは，いずれもアルゴリズム上の工夫によって回避可能
• Bi-CG法の特徴
・ 残差は必ずしも n とともに単調に減少しない
・ 所要メモリ量と1反復あたりの演算量は n に依存しない
# GMRES法と比べた場合の大きな長所
# 3項間漸化式により計算が行えるため
・ A と AT の両方による乗算が必要
# 反復当たりの行列ベクトル積の回数はGMRES法の2倍
（44/61）


# Page. 46

![Page Image](https://bcdn.docswell.com/page/57GVXX917L.jpg)

Bi-CG法から派生した解法
• CGS法
・ Bi-CG法を変形し，rn = (Rn(A))2b, sn = b* となるように工夫した手法
# snTrn は不変なので，様々な係数は Bi-CG 法と同じに計算できる
・ 残差多項式が2乗されるため，収束が速くなると予想される
# ただし，収束の不規則性も2乗で効くため，Bi-CG法より不安定
・ A による乗算が２回となり，AT による乗算は不要
• Bi-CGSTAB法
・ Sn(z) = (1 – 0z)(1 – 1z)・・・(1 – n–1z) という形の n 次多項式 Sn(z)
（安定化多項式）を用いて，rn = Sn(A)R(A)b, sn = b* となるように工夫し
た手法
# この場合も，様々な係数は Bi-CG 法と同じに計算できる
・ 各ステップで局所的に残差を最小化するように n を決めることで，
CGS法よりも収束が安定化される
（45/61）


# Page. 47

![Page Image](https://bcdn.docswell.com/page/4EQ688WNEP.jpg)

Bi-CG法から派生した解法
• GPBi-CG法
・ Bi-CGSTAB法における安定化多項式 Sn(z) として，３項漸化式によ
り生成される多項式を用いた方法
・ 一般に，Bi-CGSTAB 法よりもさらに優れた収束性を持つ
• QMR法
・ Bi-CG法において，残差を擬似的に最小化するようにした方式
・ 一般に，残差の収束は Bi-CG 法よりもスムーズ（振動が少ない）
# 収束は，Bi-CG法に比べて必ずしも速くない
# しかし，残差の振動が少ないほうが丸め誤差の点で有利
（46/61）


# Page. 48

![Page Image](https://bcdn.docswell.com/page/KJ4422K3J1.jpg)

前処理
（47/61）


# Page. 49

![Page Image](https://bcdn.docswell.com/page/LE14MMXZEG.jpg)

前処理の概念
• 前処理
・ 行列 K が，行列 A を何らかの意味で近似しているとする
・ このとき，方程式 K–1Ax = K–1b を考えると，その係数行列 K–1A は A
と比べて単位行列に近くなっていると考えられる。これにより，クリロ
フ部分空間法の収束はより速くなると期待できる
・ このような方程式の変形を前処理（左前処理）と呼ぶ
・ 各反復では，Ap の代わりに K–1Ap を計算する
• 右前処理，左前処理，
・ 右前処理： AK–1 y = b， x = K–1 y
・ 両側前処理： K1–1AK2–1 y = K1–1b， x = K2–1 y
• 前処理行列 K に求められる性質
・ K–1A がなるべく単位行列に近いこと
・ K の生成が容易であること
・ K–1x の計算が容易であること
（48/61）


# Page. 50

![Page Image](https://bcdn.docswell.com/page/GEWX33N672.jpg)

疎行列に対するLU分解とフィルイン
• フィルインのレベル
– 元々の非ゼロ要素 aij に対して，そのレベルを l(aij) = 0 と定義する
– フィルインにより aij≠0 となったとき，aij のレベルを次式で定義する
l(aij) = max(l(aik), l(akj)) + 1
レベル 0
レベル 1
レベル 2
レベル 3
レベル 4
レベル 5
（49/61）


# Page. 51

![Page Image](https://bcdn.docswell.com/page/47Z6VV3R73.jpg)

様々な前処理
• 不完全LU分解
・ 行列 A のLU分解の際に，フィルインを無視して近似的に A≒LU と分
解し，K = LU とおく。これを ILU(0) 前処理と呼ぶ
# １段階のみのフィルインを許す ILU(1) 前処理，２段階のフィルイン（フィ
ルインによって生じたフィルイン）までを許すILU(1) 前処理も定義される
・ 対称正定値行列の場合は，不完全コレスキー分解型前処理（IC(0)
前処理）が使われる
・ 幅広いクラスの行列に対して効果がある
・ 一般に並列化は困難
# K–1x の計算で，前進消去 L–1x ，後退代入 U–1y が必要となるため
# 並列化向きの不完全LU分解についても，多くの研究がある
・ フィルインの絶対値の大きさによって無視するかどうかを決める
ILUT前処理，U–1 の列ノルムの大きさによって無視するかどうかを決
める ILUC 前処理もあり，悪条件の問題で有効
（50/61）


# Page. 52

![Page Image](https://bcdn.docswell.com/page/YJ6255X17V.jpg)

様々な前処理（続き）
• 近似逆行列前処理（SPAI; SParse Approximate Inverse）
・ 行列 A そのものでなく，A–1 を直接近似する行列 M を求める
・ 具体的には，M の非零要素位置を適当に決めた上で，最小化問題
minM∥I – AM∥F を最小２乗法により解く
・ 不完全LU分解に比べ，並列化が容易
・ ILU 分解と同程度の効果を出すには，一般により多くの非零が必要
• 対角スケーリング
・ 対角行列 D1, D2により A → D1–1AD2–1 と変換し，対角要素を１とする
・ 簡単で並列化も容易だが，問題によっては有効
• 行と列の並べ替え
・ 適当な並べ替えを行った後に不完全LU分解を行うと，収束性が改善
することがある
（51/61）


# Page. 53

![Page Image](https://bcdn.docswell.com/page/GJ52YYXLJ4.jpg)

波束シミュレーション
への応用
（52/61）


# Page. 54

![Page Image](https://bcdn.docswell.com/page/9E24GGX3JR.jpg)

時間依存シュレーディンガー方程式とその離散化
• 波束シミュレーションの基礎方程式
– 時間依存シュレディンガー方程式
– H： エルミート演算子（ハミルトニアン）
– 以下では簡単のため，H が , t に依存しない場合を考える
• 空間離散化
– 差分法
以下ではこちらを扱う
– Ritz-Galerkin法（有限要素法 / 基底関数による展開）
（B ： 重なり行列）
（53/61）


# Page. 55

![Page Image](https://bcdn.docswell.com/page/D7YMVV887M.jpg)

クランク・ニコルソン法
• 導出
– (1,1)次パデ近似： exp(iHt) ≒ (I – (1/2)iHt) –1(I + (1/2)iHt )
– 標準的な導出法： 時刻 t + (1/2)t での中心差分
• 特徴
– 時間発展のために連立1次方程式を解くことが必要
– 時間に関して2次精度
– 波動関数のノルムを厳密に保存
• 時間発展演算子の固有値の絶対値は
シュレディンガー方程式の時間発展の計算に最適
（54/61）


# Page. 56

![Page Image](https://bcdn.docswell.com/page/VENWMMQ9E8.jpg)

クランク・ニコルソン法に現れる連立1次方程式
• 方程式
• 係数行列の性質
– 大規模・疎行列・非エルミート
– i 倍すると A + I （A： 正定値エルミート行列， ： 複素数）の形
• 数値解法
– 直接法： ガウスの消去法（帯行列向け解法，スパースソルバなど）
– 反復法： クリロフ部分空間法（GMRES法，BiCG法など）
係数行列の構造を利用した効率的なアルゴリズムの設計を考える
（55/61）


# Page. 57

![Page Image](https://bcdn.docswell.com/page/Y79XNND3J3.jpg)

A + I に対するGMRES法
• A のエルミート性の利用
– GMRES法は，クリロフ部分空間 Kn(A; b) = span{b, Ab, …, An–1b}
（n = 1, 2, …）の中で残差を最小にする解を探索
– クリロフ部分空間のシフト不変性 Kn(A+I; b) = Kn(A; b) より，エルミ
ート行列 A に対するクリロフ部分空間が使える
– ランチョス法と同様，直近の3本の基底ベクトルのみを用いて
Kn(A+I; b) の正規直交基底を生成可能
• 得られる解法の特徴
– 3項間漸化式に基づく解法 → 演算量・メモリ量の負担軽減
– 残差最小性を持つ → スムーズな収束性
• 問題点
– 前処理を行うと，シフト不変性が崩れてしまう
– 前処理なしで，どの程度の収束性が得られるかが鍵
（56/61）


# Page. 58

![Page Image](https://bcdn.docswell.com/page/G782KK3L7D.jpg)

数値実験結果
• GMRES法は COCG法（Bi-CG法系統の解法）に比べてスムーズな収束
を実現．計算量は同程度．
[1] H. Seito, T. Hoshi and Y. Yamamoto: On Using the Shifted Minimal Residual Method for Quantummechanical Wave Packet Simulation, JSIAM Letters, Vol. 11, pp. 13-16 (2019).
[2] S. Hidaka, S. Kudo, T. Hoshi and Y. Yamamoto, Performance of the shifted minimal residual
method for multiply shifted linear systems with real symmetric or complex Hermitian coefficient
matrices, Computer Physics Communications, to appear.
（57/61）


# Page. 59

![Page Image](https://bcdn.docswell.com/page/L7L2ZZVQJR.jpg)

終わりに
（58/61）


# Page. 60

![Page Image](https://bcdn.docswell.com/page/4EM8VVNKEW.jpg)

終わりに
• 本発表では，まずクリロフ部分空間とその基底の生成法を説
明し，方程式 Ax = b をクリロフ部分空間に射影する３種類
の方法を述べた
• これに基づき，GMRES法，CG法，Bi-CG法の３つのアルゴ
リズムと，その変種を説明した
• また，収束を加速させるための前処理法について，代表的な
手法を説明した
• 実際の問題にクリロフ部分空間法を適用するには，係数行
列の分類を検討した後，Lis や PETScなどのパッケージを
使って様々な解法・前処理を試してみるのがよい
（59/61）


# Page. 61

![Page Image](https://bcdn.docswell.com/page/PER522D679.jpg)

参考文献
（60/61）


# Page. 62

![Page Image](https://bcdn.docswell.com/page/P7XKDD1DEX.jpg)

参考文献
•
L. N. Trefethen and D. Bau III: “Numerical Linear Algebra”, SIAM,
Philadelphia, 1997.
•
H. A. van der Vorst: “Iterative Krylov Methods for Large Linear Systems”,
Cambridge University Press, Cambridge, 2003.
•
Y. Saad: “Iterative Methods for Sparse Linear Systems”, 2nd Ed., SIAM,
Philadelphia, 2003.
•
R. Barrett et al.: “Templates for the Solution of Linear Systems: Building
Blocks for Iterative Methods”, SIAM, Philadelphia, 1994.
•
藤野清次, 張紹良: “反復法の数理”, 朝倉書店, 1996.
•
杉原正顯, 室田一雄: “線形計算の数理”, 岩波書店, 2009.
•
櫻井鉄也, 松尾宇泰, 片桐孝洋： “数値線形代数の数理とHPC”, 共立出版, 2018.
•
張紹良編著： “計算科学のための基本数理アルゴリズム”, 共立出版, 2019.
（61/61）


