数理解析研究所講究録 1441 巻 2005 年 89-10289特異な行列のための改良 ORTHOMIN(m) 法
張 紹良 (Shao-Liang ZHANG) **阿部 邦美 (Kuniyoshi ABE) **岐阜聖徳学園大学経済情報学部 (Gifu Shotoku university)**1

東京大学大学院工学系研究科 (University of Tokyo)

はじめに$n\cross n$の特異な行列

$A$

を係数,

$n$

次ベクトル

$b$

を右辺項とする線形方程式(1.1)$Ax=b$

を Krylov 部分空間法によって解くことを考える. このような特異な線形方程式を解く場合,解が存在しない可能性があるので, 式 (1.1) の代わりに最小二乗問題(1.2)$\min||b-Ax||_{2}$$X\in \mathbb{R}^{n}$は常に一意に決を解くことにする. 式 (1.2) に対する解を とするとき, 残差が成り立つ. この r- を最小残差と呼ぶことにする.
まり,法) [6], Biこのような問題を解く場合, Conjugate Gradient method (共役勾配法,
Conjugate Gradient method (双共役勾配法, Bi-CG 法) [4] などを適用すると, 残差は一$\overline{7}\cdot=b-A\overline{x}$$\overline{x}$$A^{\mathrm{T}}\overline{r}=0$$\mathrm{C}\mathrm{G}$

旦減少する傾向を示すが, 最小残差に達することなぐ発散してしまうことが知られて る法などのように残差が振動する解法ではな[8]. したがって, 最小二乗解を求める場合,$\mathrm{t}\mathrm{a}$$\mathrm{C}\mathrm{G}$く,

残差ノルムの単調減少性が数学的に保証されている解法を適用することが望まし帆 そのような性質をもつ解法の一つとして, Generalized Conjugate Residual method (一般化共
役残差法, GCR 法) [3] がある. GCR 法は反復回数の増加にともな , 演算量, 記憶容量が増加するため, $m$ 回反復した後に得られた近似解をあらためて初期値として, 再度 反復を行なうリスタート版 (restarted version), または反復過程で求められた最新の $m$ 個のベクトルなどの情報を記憶しておく切断版 (truncated version) が使用される. リスタート版
は GCR(m) 法, 切断版は ORTHOMIN(m) 法 [7] として知られて る. これらの解法の残差(直は長い漸化式 (直前から $m$ 戦前の項まで) を用いて更新されるのに対し, 短 ‘漸化式法) [3] もよmethod (共役残差法,
新される Conjugate
前の 2 項) によっ
法
く知られている解法である. これらの解法のうちで切断版である ORTHOMIN(m) 法や$\mathrm{A}\mathrm{a}$$\mathrm{t}\backslash$$\mathrm{b}$$\mathrm{C}\mathrm{R}$$-c\Leftrightarrow$$\mathrm{R}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{d}\iota\iota \mathrm{a}\mathrm{l}$$\mathrm{C}\mathrm{R}$を特異な線形方程式に適用した場合, 解法の収束性は数学的な理論と異なることが報告され

ている.

すなわち,

残差は最小残差ノルムに収束することが数学的に保証されて ‘る $[5, 8]$にもかかわらず, 実際の有限桁演算においては,$\mathrm{t}$残差が最小残差ノルムに収束した後に, 丸

80
め誤差の影響によって最小残差ノルム以下になるという現象があらわれる. そのため, 近年
法と数学的に同値で理論通りに収束するアルゴリズムが提案されている [1]. しかでは,し, ORTHOMIN(m) 法については, 理論通りに収束するアルゴリズムが開発されていな .
$\mathrm{C}\mathrm{R}$$\mathrm{A}\backslash$

そこで, 本論文では, 従来の ORTHOMIN(m) 法を改良することによって新たなアルゴリ

ズムを提案する. その解法は, ORTHOMIN(m) 法と数学的に同値であるが, 新たな漸化式の係数, および補助ベクトルを定義するため, ORTHOMIN(m) 法とは異ったアルゴリズムとなる. また, 提案するアルゴリズムと従来の ORTHOMIN(m) 法の一回反復当たりの計算量は同じである. さらに, 数値実験において, 従来の ORTHOMIN(m)’ 法の残差が理論と矛盾する収束振舞いであるのに対し, 提案するアルゴリズムの残差は理論通りに収束することを示す.

まず本論文 2 節では, ORTHOMIN(m) 法のアルゴリズム, および数学的な収束性について述べる. 次に, 3 節では, 新たな漸化式の係数, および補助ベクトルを定義することによって改良版アルゴリズムを導く. さらに, これらのアルゴリズムの反復一回当たりの計算量を比較する. 4 節では, 理論的な収束条件を満たす特異な行列と収束条件を満たさない特異な行列を取り上げ, 残差の収束特性を比較する. そして, われわれが提案したアルゴリズムが従来の ORTHOMIN(m) 法より有効であることを示す.
本論文では以下の記号を用いる.
任意の正方行列を $X$ とするとき,

...

2

.—

Orthomin(m) 法
本節では, ORTHOMIN(m) 法のアルゴリズム, および特異な線形方程式に対するOR-

THOMIN (m) 法の収束定理を述べる.

まず, ORTHOMIN(m) 法のアルゴリズムを記述する [7].

を用意し, $r_{0}=b$ -Ax。を計算する. $\beta_{-1,j}=0$
アルゴリズム 1 (Orthomin(m)) :初期値に対して残差が収束するまで以下を繰り返す.
とする. $k=0,1,$$x_{0}$$\ldots$$q_{k}=Ar_{k}+ \sum_{j=k-m}^{k-1}\beta_{k-1,j}q_{j}$pk=rk+j=\Sigma k - lm\beta k-l,jpj フ$\alpha_{k}=\frac{(r_{k},q_{k})}{(q_{k},q_{k})}$

,

(2.1)
(2.2)

,

xk+l=xk+\mbox{\boldmath $\alpha$}kpk ラ

(2.3)

,

(2.4)

$r_{k+1}=r_{k}-\alpha_{k}q_{k}$91$(k-m+1\leq j\leq k)$$\beta_{k,j}=-\frac{(Ar_{k+1},q_{j})}{(q_{j},q_{j})}$.

次に, 特異な線形方程式に対する ORTHOMJN(m) 法の収束定理について述べる. まず, 係の列べクトルを部分空間の行列
数行列 $A$ の階数を $r(n>r)$ とする. ここで,を並べたものとする. このとき, 次のような収束定
range (A) の正規直交基底 ,$Q_{1}$$n\mathrm{x}r$91$q_{2},$

理が明らかにされている $[5, 8]$.

定理 1 (S.-L. Zhang, et al.)$\ldots,$$q_{r}$

CI-CII b よ同値である.

(CI) 任意の右辺項 , および初期近似解
かつ収束する.
$b$

(CII)
注意 1

$M(Q_{1}^{\mathrm{T}}AQ_{1})$null(A)かつ rankAが定値, かつ null(A)定理 2(K. Hayami)に対して, ORTHOMIN(m) 法は破綻せず,$=\mathrm{n}\mathrm{u}11(A^{\mathrm{T}})$である.

が定値」であることと 「 $Q(A)$ が半定値,」 であることは同値である.

$=\mathrm{n}\mathrm{u}\mathrm{U}(A^{\mathrm{T}})$$=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}M(A)$$x_{0}$のとき, 「$b\in \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}(A)$$M(Q_{1}^{\mathrm{T}}AQ_{1})$のとき,CI-C2 は同値である.

(C1) 任意の初期近似解 x。に対して, ORTHOMIN(m) 法は破綻せず, かつ残差が 0 に収束する.
が定値である.
(C2)$M(Q_{1}^{\mathrm{T}}AQ_{1})$定理 1, 2 は, ORTHOMIN(m) 法を特異な線形方程式に適用した場合に残差が最小残差ノルムに収束することを示している. しかし, 実際の有限桁演算においては, 丸め誤差の影響
で最小残差ノルム以下になる現象が現われる. すなわち, 定理 1, 2 で示されている数学的な結果と矛盾することがある. そこで, 次節では, 最小残差ノルムに収束するようなアルゴ
リズムを提案する.

3

Orthomin (m) 法の改良
本節では, アルゴリズム 1 における漸化式, およびパラメータを異なった補助ベクトル,

パラメータを用いて書き変えることによって, 新たなアルゴリズムを導く.
ORTHOMIN(m) 法におけるパラメータ$\alpha_{k},$$\alpha_{k}=$$\beta_{k-1,j}=$$\beta_{k-1,j}$を$\zeta_{k},$\mbox{\boldmath $\zeta$}kプ$\frac{\zeta_{j}}{\zeta_{k}}\eta_{k,j+1}$と定義する. また, 新たな補助ベクトル$y_{k}$$:=\zeta_{k-1}q_{k-1}$$z_{k}$$:=\zeta_{k-1}p_{k-1}$$\eta_{k,i+1}$

によって

82
を導入する. すると, 残差$r_{k}$

を更新するための漸化式 (2.1), (2.4) は次のように書き変えられる.

$y_{k+1}$$r_{k+1}$

さらに, 近似解$oe_{k}$$= \zeta_{k}Ar_{k}+\sum_{j=k-m+1}^{k}\eta_{k,\tilde{g}}y_{j}$

(3.1)

,

(3.2)

$=r_{k}-y_{k+1}$ .

を更新するための漸化式 (2.2) $-(2.3)$ は$z_{k+1}$$= \zeta_{k}r_{k}+\sum_{j=k-rn+1}^{k}\eta_{k,j^{Z}j}$$x_{k+1}$$=x_{k}+p_{k+1}$

(3.3)

,

(3.4)

と書き変えられる.
が満たす性質, およびパラメータ , \eta k,\iota - の計算方法を導く. 係と同様, 残差ノルムを最小化するよは ORTHOMIN(m) 法のパラメータ

次に, ベクトル$r_{k},$$\zeta_{k}$$y_{k}$

数
うに決められる.
$\zeta_{k},$$\eta_{k,j}$$\alpha_{k},$$\beta_{k-1,j}$$\min||r_{k+1}||_{2}=\min_{\zeta_{k},\eta k,j}||r_{k}-\zeta_{k}Ar_{k}-\sum_{j=k-m+1}^{k}\eta_{k,j}y_{j}||_{2}$

したがって, ベクトル$r_{k},$$y_{k}$(3.5).

は, 次のような性質を満たす.

$(r_{k+1}, Ar_{k})$$=0$ ,$(r_{k+1}, y_{j})$

$=0$

(3.6)
$(j=k-m+1, \ldots, k)$

.

(3.7)

$y_{j})=(y_{k+1}, y_{j})(j=k-m+1, \ldots, k)$ を満たすことがわ
また, 式 (3.2), (3.7) より
かる. それゆえ, 式 (3.1), (3.6)-(3.7) を利用すれば, 次のような性質を導くことができる.
$\langle$$r_{k},$$(r_{k},y_{j})=(y_{k+1}, y_{j})=0$(3.8)$(j=k-m+1, \ldots, k)$ .

さらに, 内積 $(y_{k}, y_{k})$ は, 式 (3.1), (3.2) を用いて$(y_{k},y_{k})$$(r_{k-1}-r_{k}, \zeta_{k-1}Ar_{k-1}+\sum_{j=k-m}^{k-1}\eta_{k-1,j}y_{j})$$=$,と変形でき, 性質 (3.6)-(3.8) を利用して$(y_{k},y_{k})$と整理できる. 以上, ベクトル$y_{k}$(3.9)

は性質 (3.6)-(3.9) を満たす.

の具体的な計算方法について述べる. 式 (3.5) より, パラメータは残差ノルムを最小化するように決められるので, ORTHOMIN(m) 法のパラメータ

パラメータ$\eta_{k,j}$$r_{k},$$=\zeta_{k-1}(r_{k-1}, Ar_{k-1})$$\zeta_{k\prime}\eta k,j$$\zeta_{k}$$\alpha_{k}$

,
,

88
$\beta_{k-1,j}$と同様,$\frac{\partial}{\partial\zeta_{k}}||r_{k}-\zeta_{k}Ar_{k}-\sum_{j=k-m+1}^{k}\eta_{k,j}y_{j}||_{2}^{2}=0$

,

$=0$

$\frac{\partial}{\partial\eta_{k,j}}||r_{k}-\zeta_{k}Ar_{k}-\sum_{j=k-m+1}^{k}\eta_{k,j}y_{j}||_{2}^{2}$を計算することによって求められる. これらを陽に記述するとき, 性質 (3.6)-(3.9) を用いて整理すれば, パラメータ$\zeta_{k}$$\eta_{k,j}$$\zeta_{k},$$\eta_{k,j}$は$=$$\frac{(Ar_{k},Ar_{k})}{(Ar_{k},r_{k})-\sum_{j=k-m+1}^{k}\frac{1}{\nu_{\dot{\mathit{3}}}}(Ar_{k},y_{j})(Ar_{k},y_{j})}$$=$$- \frac{\zeta_{k}}{\nu_{j}}(y_{j}, Ar_{k})$と書き表すことができる. ただし,,$(j=k-m+1, \ldots, k)$$\nu_{k}.--\zeta_{k-1}(Ar_{k-1}, r_{k-1})$(3.10)(3.11).

われわれの提案するアルゴリズムは次

以上, 式 (3.1)-(3.4), (3.10)-(3.11) をまとめると,のように表される.
を準備し, $r0=b-Ax_{0}$ を計算する. $y_{0}=-r_{0\prime}$
アルゴリズム 2(改良版) : 初期値に対して残差が収束するまで以下を繰り返す.
$l/0=(y_{0}, y_{0})$ とする.
$k=0,1,$$x_{0}$$\ldots$$\zeta_{k}=$$\eta_{k,j}$$=$$\frac{(Ar_{k_{7}}Ar_{k})}{(Ar_{k},r_{k})-\sum_{j=k-m+1}^{k}\frac{1}{\mathcal{U}j}(Ar_{k},y_{j})(Ar_{k},y_{j})}$$- \frac{\zeta_{k}}{\nu_{j}}(y_{j}, Ar_{k})$(for $k=0$ ,$\nu_{k+1}$$(k-m+1\leq j\leq k)$ ,$\zeta_{k}=\frac{(Ar_{k},r_{k})}{(Ar_{k},Ar_{k})}$$=\zeta_{k}(Ar_{k}, r_{k})$,$\eta_{k}=0$),$z_{k+1}= \zeta_{k}r_{k}+\sum_{j=k-m+1}^{k}\eta_{k,j^{Z}j}$,$x_{k+1}=x_{k}+z_{k+1}$ ,$y_{k+1}= \zeta_{k}Ar_{k}+\sum_{j=k-m+1}^{k}\eta_{k,j}y_{j}$$r_{k+1}=r_{k}-y_{k+1}$,.

このアルゴリズムを ORTHOMIN(m) 法の$\mathrm{A}\mathrm{Z}$版と名付け, AZ-ORTHOMIN(m) と略す.

AZ-ORTHOMIN(m) 法の導出過程から明らかなように, AZ-ORTHOMIN(m) 法と従来の ORTHOMIN(m) 法は数学的に同値である. したがって, 特異な線形方程式に対する AZ-ORTHOMIN(m)法の収束定理として定理」, 2 が成り立つ.

84
次に, 従来の ORTHOMIN(m) 法と AZ-ORTHOMIN(m) 法に関する反復一回当たりの行列とベクトルとの積, 内積 ベクトル和やスカラー倍の演算量を Table 1 に示す. 従来の ORTHOMIN m) 法と AZ-ORTHOMIN (m) 法の演算量は同じであることがわかる.
$\langle$

Table 1. Computational costs per iteration

add or scaling: adding a vector to another vector or scalar multiplication of a vector.

4

数値実験

本節では, 特異な非対称行列の問題を取り上げ, 提案したアルゴリズムが ORTHOMIN(m) 法)
PC-AT 互換機 (Penlum Xeon
よりも有効であることを示す. すべての数値実験$3.06\mathrm{G}\mathrm{H}\mathrm{z}$$l\lambda$

において富士通 Fortran コンパイラの倍精度演算によって実行された. また, ORTHOMIN(m)とした. さらに, ORTHOMIN(m) 法
法と AZ-ORTHOMIN(m) 法の初期ベクトルは,と AZ-ORTHOMIN(m) 法の切断係数 $m$ を 50 とし, 実用的な値より大きく設定した. その理由は, 提案したアルゴリズムの丸め誤差に対する影響を調べるためである.
$x_{\text{。}}=0$

4.1

モデル問題

特異な非対称行列を係数にもつモデル問題として, 正方領域 $\Omega=(0,1)\mathrm{x}(0,1)$ において,
次の移流拡散方程式の離散近似解を求める [2].
$\Delta u+d\frac{\partial u}{\partial x}=f(x, y)$$0<x,$ $y<1$ .

,

(4.1)

ただし, 境界条件は 41.1, 4.12 小簾で述べられる 2 種類を与え, 定理 1 の収束条件を満たす問題, および定理 1, 2 の収束条件を満たさない問題を扱う.

4.1.1

null(A)$=\mathrm{n}\mathrm{u}11(A^{\mathrm{T}})$の場合

式 (4.1) に周期境界条件$u(x, 0)=u(x, 1)$ ,$u(0,y)=u(1, y)$

を $x,$ 方向ともに $m$ 等分した差分格子を考え, 式を課す [2]. この境界値問題に対して,(4.1) を 5 点中心差分で離散近似する. こうして得られた係数行列 $A$ はの行列,$\Omega$$y$$m\mathrm{x}m$$T_{m}$$I_{m}$

85を用いて, 以下のような $n\mathrm{x}n(n=m^{2})$ の非対称行列で表せる.

$..I_{m}...\cdot$$T_{m}I_{m}I_{m}\ovalbox{\tt\small REJECT}$$A:= \frac{1}{h^{2}}\ovalbox{\tt\small REJECT} I_{m}T_{m}I_{m}$

(4.2)

.

$..I_{m}^{\cdot}..\cdot$

ただし,$\alpha_{-}\alpha_{+}-4\ovalbox{\tt\small REJECT}$$T_{m}:=\ovalbox{\tt\small REJECT}-\alpha_{-}4\alpha_{+}\alpha_{+}..\cdot.\cdot$.

$I_{m}l\mathrm{h}m\mathrm{x}m\text{の^{}\backslash \backslash }R\Phi’l^{\sim}3^{\cdot}P^{\mathrm{I}\mathrm{J}\text{を表わす}}\cdot \text{ま}.arrow,$$m=100\text{と}\ovalbox{\tt\small REJECT} \text{定し},$$|-‘-,$$\text{て_{}\backslash }\backslash ,$,$\cdot.\alpha_{-}...\cdot$$h:— \frac{1}{m’}\alpha_{\pm}:=1\pm\underline{d}hT\text{とする}$.

$\text{数}l_{\mathrm{L}}^{\mathrm{g}}\text{実}\Re^{\text{て_{}\backslash }\backslash }\iota\mathrm{h}$$d=0.5,1.5 \text{として}2\ovalbox{\tt\small REJECT}\ovalbox{\tt\small REJECT} \text{の異}rx\text{る}\Re_{\backslash }\text{数}\mathrm{q}_{\overline{\mathrm{J}}}fi\mathrm{I}\rfloor l^{\mathrm{r}}-backslash \mathrm{f}^{\text{し}}-\tau^{\simeq}\frac{-}{\beta}+\text{算を_{}l\overline{\mathrm{J}}\mathrm{J}\vee\supset}^{\prime\gamma_{-}}\cdot$$e:=(1, \ldots, 1, \ldots, 1)^{\mathrm{T}}\text{と定}\Leftrightarrow \text{する}$$A^{\mathrm{T}}e=0\text{と}r_{\zeta}\mathrm{r}\text{るのて_{}\backslash }\backslash ,$$\lambda\backslash \mathrm{f}\mathrm{f}\mathrm{f}^{\mathrm{A}};\iota_{\mathrm{O}}\beta\mathrm{B}^{\mathrm{I}^{\backslash }}\not\simeq \mathrm{j}\mathrm{E}\text{定}\{_{\llcorner}^{\mathrm{g}\text{て_{}\iota}\backslash },$.

$arrow \text{の}\text{とき},$$\text{式}(4.2)\text{の}\mathrm{q}_{\overline{\mathrm{T}}}p_{\mathrm{I}}\rfloor Af.\mathrm{x}_{\iota}\mathfrak{x}\text{し}arrow\zeta,$$\mathrm{n}\mathrm{u}\mathbb{I}(A)=\mathrm{n}\mathrm{u}11(A^{\mathrm{T}})=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\{e\}\mathrm{B}^{1^{\tau}}ffi\text{り}\Phi\cdot\supset$

rankA

.

$\text{さら}\}^{arrow}.,$$=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}M(A)=n-1\text{て^{}\tau^{\backslash }}\text{ある}arrow’ \text{と}\mathrm{B}^{1^{\backslash }}\Phi^{\frac{-}{\overline{w}}\ovalbox{\tt\small REJECT}_{\theta\backslash }-\subset^{\backslash \text{き}\mathrm{g})}}\backslash$

.

.

,

$Ae=$

$\tau’\tilde{\tau}F^{1\mathrm{J}}(4.2)\text{の}$.

$\text{し}^{\wedge}.\mathrm{B}^{1^{\backslash }}\cdot\supset \text{て}$,$\mathrm{f}’\overline{\mathrm{T}}F|\mathrm{J}(4.2)\mathrm{F}\mathrm{h}\text{定理}1\text{の}(\mathrm{C}\mathrm{I}\mathrm{I})\text{を}\backslash \text{満}\vee-\text{す}$4L2$\mathrm{n}\mathrm{u}11(A)\neq \mathrm{n}\mathrm{u}11(A^{\mathrm{T}})$の場合

式 (4.1) に Nuemann 境界条件$\frac{\partial u}{\partial n}|_{\partial\Omega}=0$

を $x,$ 方向ともに $m$ 等分した差分格子を考え, 式を課す [2]. この境界値問題に対して,の行列 &,(4.1) を 5 点中心差分で離散近似する. こうして得られた係数行列 はを用いて, 以下のような $n\mathrm{x}n(n=m^{2})$ の非対称行列で表せる.
$\Omega$$y$$A$$\cdot I_{m}I_{m}^{\cdot}$$2^{\cdot}I_{m}\mathit{3}_{m}^{\cdot}$$\mathit{3}_{m}I_{m}$

ただし,

.

$S_{m}:=\ovalbox{\tt\small REJECT}\alpha_{-}-4.-.42$.

$\alpha_{+}.\alpha_{-}.$$\ovalbox{\tt\small REJECT}$

(4.3)

.

.

.

$A:= \frac{1}{h^{2}}\ovalbox{\tt\small REJECT}^{\mathit{3}_{m}}I_{m}2.I_{m}S_{m}.$$m\mathrm{x}m$$\cdot-\cdot 42^{\cdot}\alpha_{+}-4\ovalbox{\tt\small REJECT}$.

88
4.11 小節と同様,$I_{m}$

は $m\mathrm{x}m$ の単位行列,$\alpha\pm=1\pm-dhT$$h= \frac{1}{m},$とする. 数値実験では,$d=0.5$ として計算を行なった.

$m=1\mathrm{O}\mathrm{O},$

ここで,$D_{m}:=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(1,2/\alpha_{-}, 2\alpha_{+}/\alpha_{-}^{2}, \ldots, 2\alpha_{+}^{m-3}/\alpha_{-}^{m-2}, \alpha_{+}^{m-2}/\alpha_{-}^{m-2})\in \mathbb{R}^{m\mathrm{X}m},$ $D:=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}$と定義する. このとき, 行列 (4.3) に対して, $Ae=0$ ,$\{e\},$
および $A^{\mathrm{T}}De=0$ となるので, null(A){De} が成り立つ. す$(D_{m}, 2D_{m}, \ldots, 2D_{m}, D_{m})\in \mathbb{R}^{n\cross n}$$\mathrm{n}\iota\iota \mathrm{U}(A^{\mathrm{T}})=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}$$=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}$なわち, 行列 (4.3) は, 定理 1 の (CII) を満たさない.

4.13

計算条件

行列 (4.2), (4.3) を係数にもつ線形方程式では, 右辺項 b はと分解して考えることができる. すると, 右辺項 $b=b_{R}+b_{N}$ において,$b=b_{R}+b_{N}(b_{R}\in \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}(A),$$\mathrm{n}\mathrm{u}11(A^{\mathrm{T}}))$$b_{N}\in$$b_{N}\neq 0$のとき, 反復解法を適用した場合,$||r_{k}||_{2}^{2}=||b-Ax_{k}||_{2}^{2}=||b_{N}||_{2}^{2}+||b_{R}-Ax_{k}\ovalbox{\tt\small REJECT}$に安定収束することが望ましい. こ
は最小残差ベクトルのような右辺ベクトル をもつ問題を数値実験の対象とする.
一様乱数によって金を生成し, $b_{R}=A\hat{x}$ によって右辺項をつくり, 人工的に最小残差ベク
トルを加える. すなわち, 行列 (4.2) に対しては式 (4.4), 行列 (4.3) に対しては式 (4.5)となるので, 残差ベクトル$b_{N}$$r_{k}$$b$となる.
$Ax=b_{R}+ \delta\frac{e}{||e||_{2}}$A2,$=b_{R}+ \delta\frac{De}{||De||_{2}}$

(4.4)

.

(4.5)

として数値実験を行なう. このとき, 収束停止条件は残差ノルムに達した時である. しかし, 実際の応用分野で現れる問題では, 最小残差が最小残差が明らかでないことが多い. このような場合, 通常, 停止条件を設けず数千回の反復を継続する. そこで, われわれは反復を 3000 回まで継続して行なうことにする.
ただし,$\delta=10^{-6}$$||r_{k}||_{2}$$10^{-6}$

4.2

計算結果

線形方程式 (4.4), (4.5) に従来の ORTHOMIN(m) 法, AZ-ORTHOMIN(m) 法を適用する.
式 (4.4) の場合, 定理 1 の収束条件を満たすので, 残差が最小残差ノルムに収束することが理論的に保証される. 一方, 式 (4.5) の場合, 定理 1, 2 の収束条件を満たさないので,理論的に残差が最小残差ノルムに収束するとは言えない. 本小節では, これらの問題に対する ORTHOMIN(m) 法,-ORTHOMIN(m) 法の残差ノルムの収束振舞いを考察する.
$d=0.5,1.5$ とした行列 (4.2) を係数にもつ線形方程式 (4.4) に ORTHOMIN(m) 法, AZORTHOMIN(m) 法を適用したときの収束特性をそれぞれ Figs. 1, 2 に示す. グラフの横軸
は反復回数縦軸は残差ノルム $(\log_{10}(||r_{k}||_{2}))$ を表す. ただし, ベクトル はアルゴリズ”
は, それムにおいて漸増的に計算された残差である. グラフ中の” ”,
ぞれ AZ-ORTHOMIN (m) 法, 従来の ORTHOMIN (m) 法を意味する.
$\mathrm{A}\mathrm{Z}$$r_{k}$$\mathrm{A}\mathrm{Z}$$7’ \mathrm{O}\mathrm{R}\mathrm{T}\mathrm{H}\mathrm{O}\mathrm{M}\mathrm{I}\mathrm{N}(\mathrm{m})$87Fig. 1. The residual 2-norm history of AZ-ORTHOMIN(50) and ORTHOMIN(50) for thecoefficient matrix (4.2) with $d=0.5$ .

FIigg. 2. The residual 2-norm history of AZ-ORTHOMIN(50) and ORTHOMIN(50) for thecoefficient matrix (4.2) with $d=1.5$ .

88に到達したときの反復回数, および 10
り小さくなり始めたときの反復回数を Table 2 に示す.
また, 従来の ORTHOMIN(m) 法の残差が 10

Table 2. The iteration counts at which$||rk||_{2}$becomes 10

よor less in Figs. 1-2.

さらに, $d=0.5$ とした行列 (4.3) を係数にもつ線形方程式 (4.5) に ORTHOMIN(m) 法,

AZ-ORTHOMIN(m) 法を適用した場合の収束特性を Fig, 3 に示す.

Fig. 3. The residual 2-norm history of AZ-ORTHOMIN(50) and ORTHOMIN(50) for thecoefficient matrix (4.3) with $d=0.5$ .
察].

の場合の考察を行なう. 従来の ORTHOMIN(m) 法の残差は最小
残差ノルム 10 以下となり, 理論と反する結果となった. 一方, われわれが提案したアル
ゴリズムの残差は, 理論的な最小残差ノルムに収束し, その値を維持することができた. したがって, 従来の ORTHOMIN(m) 法より AZ-ORTHOMIN(m) 法が有効であると言える.
次に,の場合の考察を行なう. 理論的には, 残差は最小残差ノルムに収束することが保証されていない. しかし, 従来の ORTHOMIN(m) 法, AZ-ORTHOMlN(m)法の残差はともに最小残差ノルムに収束し, その値を維持することができた. また, ORTHOMIN(m) 法と AZ-ORTHOMIN(m) 法の収束振舞いに差が生じなかった.
まず, nu 坦 A)$=\mathrm{n}\mathrm{u}\mathrm{U}(A^{\mathrm{T}})$$\mathrm{n}\iota 111(A)\neq \mathrm{n}\mathrm{u}11(A^{\mathrm{T}})$88
次に, 従来の ORTHOMIN(m) 法と AZ-ORTHOMIN(m) 法から得られた近似解の精度を調べるために, 真の残差ノルム $(||b-Ax_{k}||_{2})$ を比較する. ORTHOMIN(m) 法, および AZORTHOMIN(m) 法の真の残差ノルムは, 理論的には 10 に収束し, その値を維持しなければならない.

$d=0.5,1.5$ とした行列 (4.2) を係数にもつ線形方程式 (4.4) に ORTHOMIN(m) 法, AZORTHOMIN(m) 法を適用した場合の収束特性をそれぞれ Figs. 4, 5 に示す. さらに, $d=0.5$とした行列 (4.3) に対する収束特性を Fig. 6 に示す. グラフの横軸は反復回数, 縦軸は真の残差ノルム $(\log_{10}(||b-Ax_{k}||_{2}))$ を表す.

Fig. 4. The explicitly computed residual -norm history of AZ-ORTHOMIN(50) and
ORTHOMJN(50) for the coefficient matrix (4.2) with $d=0.5$ .
$2$

[考察]
まず, null(A)ルムはおよそ 10なかった. また,$=\mathrm{n}\mathrm{u}11(A^{\mathrm{T}})$または 10の$\#\mathrm{B}\mathrm{A}-$の考察を行なう. 従来の ORTHOMIN(m) 法の真の残差ノまで増加し,

理論的な最小残差ノルムを維持することができと真の残差ノルム $(||b-Ax_{k}||_{2})$漸化的に求められた残差ノルムとの収束振舞いに大きな差が生じた. 一方, AZ-ORTHOMIN(m) 法の真の残差ノルムは理論的な最小残差ノルムに収束し, その値をおおよそ維持することができた. したがって, AZORTHOMIN(m) 法から得られた近似解は従来の ORTHOMIN(m) 法より精度が良いと言える.
の場合の考察を行なう. 従来の ORTHOMIN(m) 法, AZ-ORTHOMIN(m)
次に,法の真の残差は, 共に理論的な最小残差ノルムに収束し, その値を維持することができた.
$(||r_{k}||2)$$\mathrm{n}\mathrm{u}11(A)\neq \mathrm{n}\mathrm{u}11(A^{\mathrm{T}})$

すなわち,両解法から適切な近似解が得られていると考えられる.

100$.\mathrm{A}\mathrm{Z}$$1\circ\prime \mathrm{t}\mathrm{h}\mathrm{o}\mathrm{m}i\mathrm{n}(n\urcorner)$.

’$\circ$$\mathrm{x}$

$0$

-2

.

..

$.j.\cdot\cdot\cdot \mathrm{x}-\cdot--\cdot--\kappa---\cdots- \mathrm{x}\cdot-\cdot-\cdot--\kappa---$

$i|$

.

$\acute{.}$$.\cdot.\sim’\cdot-- \mathrm{x}----\cdot\cross\cdot$

.

.

-

$\mathrm{r}$

-8

. x.
$\cdot$$\cdot$I

0500$\mathrm{t}000$

1500

2000

2500

3000

-ORTHOMIN(50)and
Fig. 5. The explicitly computed residual 2-norm history of
ORTHOMIN(5O) for the coefficient matrix (4.2) with $d=1.5$ .
$\mathrm{A}\mathrm{Z}$.

$[mathring]_{\mathrm{x}}$$,\mathrm{O}\hslash \mathrm{h}\mathrm{o}\mathfrak{m}\mathrm{i}\mathrm{n}.\langle\iota \mathfrak{n}\}\mathrm{A}\mathrm{Z}^{\cdot}$

$0$

-2

-4
$\{$

$-\cdot-\cdot-\cdot$

-8

0

500

1000

1500

2000

2500

3000

Fig. 6. The explicitly computed residual -norm history of AZ-ORTHOMIN(50) and
ORTHOMIN(50) for the coefficient matrix (4.3) with $d=0.5$ .
$2$

101
5

まとめ
われわれは 従来の ORTHOMIN(m) 法のアルゴリズムを変形することによって AZ-ORTHOMIN(m)法を導いた. その解法は従来の ORTHOMIN(m) 法と数学的に同値ではあるが, 異なる漸化式の係数や補助ベクトルを使用するため, アルゴリズムは異なる. その結果, 特異な線形方
程式に対して, 従来の ORTHOMIN(m) 法よりも次のような点で利点をもつことがわかった.
(1) ORTHOMIN(m) 法と AZ-ORTHOMIN(m) 法に関する一回反復当たりの計算量は, 同じである.

(2) m 徂$(A)=\mathrm{n}$面$(A^{\mathrm{T}})$の場合: 従来の ORTHOMIN(m) 法を使用した場合, 丸め誤差の影響

で最小残差ノルム以下になるという理論と反する現象が起こる. しかしながら, われわれが提案した AZ-ORTHOMIN(m) 法の残差は, 最小残差ノルムを維持することができた. さらに, AZ-ORTHOMlN(m) 法から得られた近似解は従来の ORTHOMIN(m) 法
より精度が良かった. したがって, AZ-ORTHOMIN(m) 法は従来の ORTHOMIN(m) 法

より有効であると言える.
(3)

いの場合 :ORTHOMIN(m) 法, AZ-ORTHOMIN(m) 法は, 収束に差を生じることなく, 共に最小残差ノルムに収束した. また, 両解法から適切な近似解が得られていた. したがって, 理論的に収束することが保証されて t‘な V‘問題に
従来の ORTHOMIN(m) 法, AZ-ORTHOMIN(m) 法を適用した場合であっても, それら$\ovalbox{\tt\small REJECT}^{-}\ovalbox{\tt\small REJECT}$$\mathrm{n}\mathrm{u}\mathrm{U}(A)\neq \mathrm{n}\mathrm{u}\mathrm{U}(A^{\mathrm{T}})$の解法の残差は収束する可能性があると言える.
謝辞

本研究の一部は科学技術振興機構 CREST の助成を受けて行なったものである.

参考文献
型の三項漸化式に基づく非対称行列のた[1] 阿部邦美, 張紹良, 三井斌友, MRT 法:
めの反復解法, 日本応用数理学会論文誌, 7(1997), 37-50.
$\mathrm{T}\mathrm{R}$$\mathrm{C}\mathrm{G}$[2] BROWN, P. N. and WALKER, H. F., GMRES on (nearly)
Matrix Anal. and APPI., 18 (1997), 37-51.

$\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{g}\iota\iota \mathrm{l}\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{y}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{s}_{?}8\mathrm{I}\mathrm{A}\mathrm{M}$J.

H. C. and SCHULTZ, M. H., Variational iterative methods[3] EISENSTAT, S. C.,
for nonsymmetric systems of linear equations, SIAM J. Numer. Anal., 20(1983), 345$\mathrm{E}\mathrm{L}\mathrm{M}\mathrm{A}\mathrm{N}_{7}$357.
[4] FLETCHER, R., Conjugate gradient methods for indefinite systems, in Numerical
G., Lecture Notes in Mathematics, 506 (1976),
Analysis Dundee 1975, ed. by$\mathrm{W}\mathrm{a}\mathrm{t}\mathrm{s}\mathrm{o}\mathrm{n}_{7}$Springer-Verlag, 73-89.
[5] 速水謙, 特異な系に対する共役残差法の収束性について, 日本応用数理学会論文誌 13(2003), 1-33.

102[6] HESTENES, M. R. and STIEFEL, E., Methods of Conjugate Gradients for Solving. Nat. $Bur$. Standards, 49 (1952), 409-435. Springer-Verlag,
Linear Systems, J.
${\rm Res}$

128-139.
[7] VJNSOM, P. K. W., Orthomin, An iterative method for solving sparse sets of simultaneous lnear equations, in Proc. Fourth Symposium on Reservoir Siwvulation, Societyof Petroleum Engineers of AIME, (1976), 149-159.
[8] ZHANG, S.-L., OYANAGI, Y. and SUGIHARA, M., Necessary and sufficient conditions
for the convergence of Orthomin(A) on singular and inconsistent linear systems, Numer.
Math., 87 $(2000),$ $391- 405$ .