統計数理（2005）
第 53 巻 第 2 号 331–348
c 2005 統計数理研究所


特集「計算推論―モデリング・数理・アルゴリズム―」［原著論文］GMRES 法による最小二乗問題の解法
速水 謙

1,2

・伊藤 徳史3（受付 2005 年 3 月 24 日；改訂 2005 年 9 月 20 日）

要旨

大規模疎な m × n 行列 A を係数行列とする最小二乗問題に対する主流の反復法は，正規方程式に対して（前処理付き）共役勾配法を適用する CGLS 法である．
本論文ではまず，代案として元の最小二乗問題を，n × m 行列 B を用いて，正方行列 AB ま
たは BA を係数行列にもつ等価な最小二乗問題に変形し，非対称正方行列を係数行列とする連
立一次方程式用のロバストなクリロフ部分空間反復法である一般化最小残差法（GMRES）を適用する手法を提案する．
次に，m ≥ n（優決定），m < n（劣決定），およびランク落ちも含めた一般の場合に対して，
これらの手法が任意の右辺項 b に対して破綻することなく最小二乗解を与えるための行列 B に関する十分条件を導く．そして，B の例として不完全 QR 分解の一つである IMGS（l）法を提案する．
最後に，フルランクな優決定および劣決定問題に対する数値実験により，条件の悪い問題では（対角スケーリングに相当する）IMGS（0）法を用いた提案手法の方が従来の前処理付き CGLS
法より速く最小二乗解を与えることを示す．
キーワード： 最小二乗問題，クリロフ部分空間反復法，CGLS 法，GMRES 法，特異な系，前処理．
1.

はじめに最小二乗問題(1.1)

min b − Ax2

x∈Rnを考える．ただし，A は m × n の実行列で，m ≥ n と m < n の両方の場合を考えるものとす
る．また，rankA ≤ min(m, n) の等号が成り立たない，いわゆるランク落ちの場合も許すものとする．
式（1.1）は正規方程式(1.2)

AT Ax = AT bと等価である．
1
2
3

国立情報学研究所 情報学基礎研究系：〒101–8430 東京都千代田区一ツ橋 2–1–2; hayami@nii.ac.jp
総合研究大学院大学 複合科学研究科情報学専攻

株式会社 ビジネスデザイン研究所：〒100–0006 東京都千代田区有楽町 1–7–1 有楽町電気ビル南館 9F;
ito@business-design.co.jp

332

統計数理 第 53 巻

第 2 号 2005

1.1 直接法
上記の問題（で m ≥ n でランク落ちのない場合）に対する従来法としては，まず直接法がある．
代表的な方法としては，まず A を，m × n の直交行列 Q (QTQ = In ）と，n × n の上三角行列 Rの積に A = QR と分解する．具体的には，Householder 法，修正 Gram-Schmidt 法，Givens 法などを用いる．すると，（1.2）式は RT Rx = RT QT b と等価である．また，rankA = n なら R はT
正則だから，Rx = Q b から後退代入により，式（1.1）の最小二乗解 x が求まる．
ただし，A が大規模で疎な場合は，メモリーと計算時間を節約する工夫が必要となる（Björck,1996）．
1.2 正規方程式を用いた反復法
大規模疎な問題に対しては，メモリーと計算時間の節約のため反復法が有効となる．反復法
による最小二乗問題の従来の解法に関しては例えば Björck（1996）を参照されたい．なお，本研究に関連する先駆的な研究として Tanabe（1974, 1975）がある．
正規方程式（1.2）の係数行列 AT A は対称な正方行列であり，rankA = n なら定値なので，正規方程式に共役勾配法を適用する下記の Conjugate Gradient Least Squares（CGLS）法が従来主流の方法である．
CGLS 法Choose x0 ．
r̃ 0 = AT (b − Ax0 )p0 = r̃ 0
for i = 0, 1, 2, . . . until r̃ i 2 < (r̃i , r̃ i )αi =(pi , AT Api )xi+1 = xi + αi pir̃ i+1 = r̃ i − αi AT Api(r̃ i+1 , r̃ i+1 )βi =(r̃ i , r̃ i )pi+1 = r̃ i+1 + βi pi
endfor

ところで，一般に長方形行列 A の一般逆行列を A† ，rankA = r ，最大，最小特異値を各々σ1 , σr とすると，A の条件数は
σ1
κ(A) := A2 A† 2 =σrで与えられる（例えば Björck, 1996 参照）．よって，
κ(AT A) = σ 2
1

σr= κ(A)2となる．
CGLS 法の収束の速さは（A がランク落ちしている場合も含めて）A の条件数 κ(A) に依存することが知られている（Björck, 1996 参照）．従って，悪条件問題に対しては CGLS 法そのものの収束は遅いので，適切な前処理が必要となる．簡便な前処理法としては，AT A の対角項を用いたスケーリングがある．その他に，例えば不完全コレスキー分解（Meijerink and van der Vorst,
1977），不完全 QR 分解（Jennings and Ajiz, 1984; Saad, 1988），不完全 Givens 直交化（Bai etal., 2001），ロバストな不完全分解（Benzi and Tuma, 2003）などの手法が提案されている．

GMRES 法による最小二乗問題の解法

333

例えば不完全 QR 分解の例として Jennings and Ajiz（1984）による下記の不完全修正 GramSchmidt（IMGS）法がある．
(1)

A = [a1 , . . . , an ], ai = ai (i = 1, . . . , n), τ を閾値パラメタとして，
for i = 1, 2, . . . , n
(i)
a
(i)
rii = ai 2 ,qi = irii
for j = i + 1, . . . , n(i)rij = q T
i aj
if |rij | < τ ai 2 then rij = 0
(i+1)
(i)
= aj − rij q iaj
endfor
endfor

また，Saad（1988）は下記のような不完全修正 Gram-Schmidt（IMGS）法を提案している．
pQ , pR をパラメタとして，
for i = 1, 2, . . . , n
(i)
a
(i)
rii = ai 2 ,qi = irii
Determine the pQ largest elements of q i and assign 0 to the other elements.
for j = i + 1, . . . , n(i)rij = q T
i aj
Determine the pR largest rij ’s for i + 1 ≤ j ≤ n and assign 0 to the others.
(i+1)
(i)
aj
= aj − rij q i
endfor
endfor

このようにして得られた上三角行列 R = (rij ) を用いて正規方程式（1.2）をÃx̃ = b̃,

(1.3)

ただし，
Ã = R−T AT AR−1 ,x̃ = Rx,b̃ = R−T AT b,と前処理し，式（1.3）に CGLS 法を適用する．
1.3 正規方程式に基づかない反復法
前処理をしたとしても，式（1.3）の条件数は AR−1 の条件数の二乗なので悪条件問題に対しては収束がまだ十分速くない可能性がある．これに対して，張（1989），Zhang and Oyanagi（1990，
1991）は正規方程式に基づかないで，A を直接扱う反復法として，n × m の写像行列 B を用いて行列 AB を係数行列とする系に対して，非対称正方行列を係数とする連立一次方程式のクリ
ロフ部分空間反復解法の一つである Orthomin（k）法を適用する，CR-LS（k）法を提案している．
本論文では，Orthomin（k）法よりも破綻しにくい GMRES 法に同様の考えを適用することから始める（伊藤・速水, 2002, 2003; 伊藤, 2003）．次に，行列 BA を係数行列とする系に対して
GMRES 法を適用する方法を提案する（Ito and Hayami, 2004; 伊藤・速水, 2004; 速水・伊藤,
2004, 2005）．そして，m ≥ n の場合だけでなく，m < n の場合や，A がフルランクでない場合も含めて，両手法が任意の b に対して破綻することなく最小二乗解を与えるための行列 B に関する十分条件を導く．最後に，フルランクな優決定および劣決定問題に対する数値実験によ334

統計数理 第 53 巻

第 2 号 2005

り両手法の有効性を示す．
2.

GMRES 法による最小二乗問題の解法

Saad and Schultz（1986）の GMRES 法（Generalized Minimal Residual method）は，対称とは限らない正則な行列を係数行列とする連立一次方程式のロバストなクリロフ部分空間法として知られる．そのままだと，反復数の増加に伴いメモリーと計算量が膨大になるので，下記のよ
うに k 反復毎にそのときの近似解を初期解としてアルゴリズムを再開する（リスタートする）
GMRES（k）法を使うことが多い．
（k = ∞ としたのが GMRES 法に相当する．
）
GMRES
（k）
法
Choose x0 ．
∗ r 0 = b − Ax0
r0v1 =r 0 2for i = 1, 2, . . . , khj,i = (Av i , v j )

(j = 1, 2, . . . , i)iv̂ i+1 = Avi −hj,i v jj=1

hi+1,i = v̂ i+1 2
v i+1 =v̂ i+1
hi+1,iFind y i ∈ Ri which minimizes r i 2 =  r 0 2 ei − H̄i y 2 .
if r i 2 <  thenxi = x0 + [v 1 , . . . , v i ]y i
stop
endif
endforx0 = xk
Go to ∗.

ただし，上記で H̄i = (hpq ) ∈ R(i+1)×i , ei = (1, 0, . . . , 0)T ∈ Ri+1 とする．
上記で hi+1,i = 0 となるとき，アルゴリズムが破綻するという．丸め誤差がなければ，A ∈RN×N が正則な場合は GMRES 法は破綻することなく N 反復以内に連立一次方程式 Ax = b の真の解を与える（Saad and Schultz, 1986）．A が特異な場合の GMRES 法の振る舞いについては，例えば Brown and Walker（1997）や 速水（2003）で論じられている．
さて，GMRES 法を直接式（1.1）の最小二乗問題に適用しようとしても，A は m × n 行列で，初期近似解 x0 に対する残差ベクトル r 0 = b − Ax0 は m 次元ベクトルなので，r 0 に A をかけてクリロフ部分空間を構成することはできない．そこで，n × m 行列 B を使ってこの問題を解決するには二つの方法が考えられる．
2.1 方法 1
まず最初の方法は張（1989），Zhang and Oyanagi（1990, 1991）のように A の右から B をかけて m × m 行列 AB を構成し，m 次元空間の中のクリロフ部分空間  r0 , ABr0 , . . . , (AB)i−1 r 0 を用いるものである．

GMRES 法による最小二乗問題の解法

335

ここで，行列 M に対して R(M ) を M の像空間，N (M ) を M の核空間，部分空間 V に対して V ⊥ を V の直交補空間として以下の準備をする．
定理 2.1. 任意の b ∈ Rm に対して minn b − Ax2 = minm b − ABz2 が成り立つための必
x∈R

z∈R

要十分条件は R(A) = R(AB) である．
証明．
（⇐=）: 自明．
（=⇒）:

R(A) = R(AB) =⇒ R(A) ⊃ R(AB) =⇒ ∃ b̃ ∈ R(A)\R(AB)⇐⇒ ∃ x̃ ∈ Rn ;b̃ = Ax̃， しかしb̃ = ABz∀z ∈ Rm⇐⇒ 0 = minn b̃ − Ax2 < minm b̃ − ABz2 .
x∈R

補題 2.2.

z∈R

R(AAT ) = R(A).

証明．




 

R(AAT ) = Ax | x ∈ R(AT ) = Ax | x ∈ N (A)⊥

= Ax | x ∈ N (A)⊥ ∪ N (A) = R(A).

（例えば，ラオ, 1977, p. 26 の別証も参照．
）
補題 2.3.

R(AT ) = R(B) =⇒ R(A) = R(AB).

証明．補題 2.2 より明らか．
例えば，rankA = rankB = n ならば R(AT ) = R(B) = Rn が成り立つ．
そこで R(AT ) = R(B) を仮定すると，補題 2.3 より R(A) = R(AB) が成り立ち，任意のx0 ∈ Rn に対して Ax0 = ABz 0 となるような z 0 ∈ Rm が存在し，r 0 = b − Ax0 = r 0 − ABz 0 となる．
そこで GMRES（k）法を使って最小二乗問題
min b − ABz2 = minn b − Ax2

z∈Rm

x∈R

を解くことを考える．その際，初期近似解 z = z 0 を ABz 0 = Ax0 となるようにとると，下記のアルゴリズムを得る．
GMRES
（k）
-LS 法 1
Choose x0 (Ax0 = ABz 0 ).
∗ r 0 = b − Ax0 ( = b − ABz 0 )r0v1 =r 0 2

336

統計数理 第 53 巻

for i = 1, 2, . . . , khj,i = (ABvi , v j )v̂ i+1 = ABvi −i


第 2 号 2005

(j = 1, 2, . . . , i)hj,i v jj=1

hi+1,i = v̂ i+1 2
v i+1 =v̂ i+1
hi+1,iFind y i ∈ Ri which minimizes r i 2 =  r 0 2 ei − H̄i y 2 .
xi = x0 + B[v 1 , . . . , v i ]y i( ⇐⇒ z i = z 0 + [v 1 , . . . , v i ]y i )r i = b − Axi
if AT r i 2 <  stop
endfor( ⇐⇒ z 0 = z k )x0 = xk
Go to ∗.

ところで，上記のアルゴリズムで GMRES（k）法の代わりに GMRES 法を用いた方法（k = ∞の場合）を GMRES-LS 法 1 と呼ぶことにする．この GMRES-LS 法 1 が破綻することなく式（1.1）の最小二乗解を与えるための必要十分条件を考える．
まず，下記の定理（Brown and Walker, 1997; 速水, 2003）に注意する．
定理 2.4. Ã ∈ Rm×m に対して下記が成り立つ．
GMRES 法が任意の b ∈ Rm および任意の初期近似解 z 0 ∈ Rm に対して破綻することなくminm b − Ãz2 の最小二乗解を与えるための必要十分条件は R(Ã) = R(ÃT ) である．

z∈R

次に下記が成り立つ．
「R(AB) = R(B T AT ) ⇐⇒ R(A) = R(B T ) 」が成り立つ．
定理 2.5. R(AT ) = R(B) ならば，証明．R(AT ) = R(B) ならば，補題 2.2 よりR(AB) = R(AAT ) = R(A),

R(B T AT ) = R(B T B) = R(B T ).

従って，定理 2.4 で Ã = AB とおくと，定理 2.5 より，下記を得る．
定理 2.6. R(AT ) = R(B) ならば下記が成り立つ．
GMRES-LS 法 1 が任意の b ∈ Rm および任意の x0 ∈ Rn に対して破綻することなく minn b−
Ax2 の最小二乗解を与えるための必要十分条件は R(A) = R(B T ) である．

x∈R

また，系として下記を得る．
系 2.7. R(AT ) = R(B) かつ R(A) = R(B T ) ならば GMRES-LS 法 1 は任意の b ∈ Rm および任意の x0 ∈ Rn に対して破綻することなく minn b − Ax2 の最小二乗解を与える．
x∈R

GMRES 法による最小二乗問題の解法

337

ここで，R(A) = R(B T ) は，B がある正則な行列 C を用いて
B = CAT(2.1)と表されることと等価である．
ところで，Calvetti et al.（2000），Reichel and Ye（2005）は GMRES 法を用いて m ≥ n の場
合の最小二乗問題を解く方法として，行列 A の右側に 0 列ベクトルを m − n 本付け足して特異行列 Ã = [A, 0] を構成し，GMRES 法を
minm b − Ãz2z∈R= minn b − Ax2

x∈Rに適用させることを提案している．
これは，上記の GMRES-LS 法 1 において B = [ In , 0 ], Ã = AB とおくことと等価である．ただし，In は n × n の単位行列である．
従って，rankA = n ならば R(AT ) = R(B) = Rn は成立するが，R(A) = R(B T ) は成立するとは限らないので，定理 2.6 より彼らの方法は最小二乗解に到達する前に破綻する可能性がある．
2.2 方法 2
もう一つの方法として，同じ n × m の写像行列 B を用いて下記のように初期残差ベクトルr 0 ∈ Rm を r̃ 0 = Br 0 ∈ Rn に写像しておいて，n 次元空間内のクリロフ部分空間 r̃ 0 , BAr̃0 , . . . ,(BA)i−1 r̃ 0  を構成し GMRES（k）法を適用する次の方法も考えられる．
GMRES
（k）
-LS 法 2Choose x0 .
∗ r̃ 0 = B(b − Ax0 )v 1 = r̃ 0 /r̃ 0 2
for i = 1, 2, . . . , khj,i = (BAvi , v j )v̂ i+1 = BAvi −i(j = 1, 2, . . . , i)hj,i v jj=1

hi+1,i = v̂ i+1 2
v i+1 = v̂ i+1 /hi+1,iFind y i ∈ Ri which minimizes r̃ i 2 =  r̃ 0 2 ei − H̄i y 2 .
xi = x0 + [v 1 , . . . , v i ]y ir i = b − Axi
if AT r i 2 <  stop
endforx0 = xk
Go to ∗.

この手法は
(2.2)

BAx = Bbに GMRES（k）法を適用するのと等価である．
ここでまず下記が成り立つ．
（この定理および証明は査読者によって大幅に改善された結果で
ある．
）

338

統計数理 第 53 巻

第 2 号 2005

定理 2.8. 任意の b ∈ Rm に対して方程式 BAx = Bb が解をもち，その解が minn b − Ax2
x∈Rの解となるための必要十分条件は R(A) = R(B T ) である．

証明．（=⇒）まず，任意の b ∈ Rm に対して方程式 BAx = Bb が解をもつための必要十分条
件はR(BA) = R(B).

(2.3)

つぎに，方程式 BAx = Bb が解をもつとして，その解が minn b − Ax2 の解となる，すなx∈R

わち，任意の b ∈ Rm に対して，B(b − Ax) = 0 ならば，AT (b − Ax) = 0 となるための必要十分条件はN (B) ⊆ N (AT ).

これを像空間で書けば，(2.4)

R(A) ⊆ R(B T ).

ここで
rankBA ≤ min{rankA, rankB},

すなわちdim R(BA) ≤ min{dim R(B), dim R(A)}に注意すれば，（2.3）より，dim R(B) ≤ dim R(A) を得る．これは，dim R(B T ) ≤ dim R(A) と同値であり，（2.4）と合わせて，R(A) = R(B T ) を得る．
（⇐=）
R(A) = R(B T ) =⇒ R(BA) = R(BB T ) = R(B).
また，
R(A) = R(B T ) =⇒ R(A) ⊆ R(B T ).

次に，定理 2.5 と双対的に下記が成立する．
定理 2.9. R(A) = R(B T ) ならば，
「R(BA) = R(AT B T ) ⇐⇒ R(AT ) = R(B)」が成り立つ．
証明．R(A) = R(B T ) ならば，補題 2.2 より
R(BA) = R(BB T ) = R(B),

R(AT B T ) = R(AT A) = R(AT ).

さらに，定理 2.4 で Ã = BA とおき，定理 2.9 を用いると下記を得る．
定理 2.10. R(A) = R(B T ) ならば下記が成り立つ．
GMRES-LS 法 2 が任意の b ∈ Rm および任意の x0 ∈ Rn に対して破綻することなく minn b−
Ax2 の最小二乗解を与えるための必要十分条件は R(AT ) = R(B) である．

x∈R

GMRES 法による最小二乗問題の解法

339

また，下記の系が成り立つ．
系 2.11. R(A) = R(B T ) かつ R(AT ) = R(B) ならば GMRES-LS 法 2 は任意の b ∈ Rm および任意の x0 ∈ Rn に対して破綻することなく minn b − Ax2 の最小二乗解を与える．
x∈R

ところで，Reichel and Ye（2005）は GMRES 法を用いて m ≤ n の場合の最小二乗問題を解く方法として，行列 A の下に 0 行ベクトルを n − m 本付け足して特異行列



Ã =

を構成し，さらに


b̃ =

A
0

b
0

として，GMRES 法を
min b̃ − Ãx2

x∈Rnに適用させることを提案している．
これは，上記の GMRES-LS 法 2 において，Im を m × m の単位行列として，



B=

Im
0

とおいて，Ã = BA,b̃ = Bb とおくことと等価である．
従って，rankA = m ならば R(A) = R(B T ) = Rm は成立するが，R(AT ) = R(B) は成立するとは限らないので，定理 2.10 より彼らの方法は最小二乗解に到達する前に破綻する可能性がある．
2.3 B の充たすべき条件
以上をまとめると，ランク落ちの場合も含めた rankA ≤ min(m, n) の一般の場合には下記が
成り立つ．
条件：
R(A) = R(B T ),(2.5)

R(AT ) = R(B)が成り立てば，系 2.7 および系 2.11 より，GMRES-LS 法 1 および 2 は任意の b ∈ Rm および任意の x0 ∈ Rn に対して minn b − Ax2 の最小二乗解を与える．ここで，条件（2.5）は B = ATx∈Rとおけば成り立つ．
次に，rankA = min(m, n) のフルランクの場合を考える．まず，m ≥ n = rankA （優決定問題）の場合を考える．任意の正則な行列 C ∈ Rn×n を用いて B = CAT とおくことを考える．
B, AT ∈ Rn×m に対して以下が成り立つ．
C ∈ Rn×n : 正則

B = CAT ,
B T = AC T ,C T : 正則

R(A) = R(B T ).

340

統計数理 第 53 巻

第 2 号 2005

従って，n = rankAT = rankA = rankB T = rankB より R(AT ) = R(B) = Rn も成り立つ．
ここで，AB ∈ Rm×m , BA ∈ Rn×n , m ≥ n だから，GMRES-LS 法 2（BAx = Bb）を用いる
方が GMRES-LS 法 1 を用いるより反復当たりの計算量が少なくてすむ．rankA = n ならば
AT A，そして diag(AT A) は正則なので，C の簡単な例として Zhang and Oyanagi（1991）のよ
うに C := {diag(AT A)}−1 とおくことが考えられる．
次に，rankA = m ≤ n（劣決定）の場合は，任意の正則な行列 C ∈ Rm×m を用いて B = AT C
とおくことを考える．B, AT ∈ Rn×m に対して以下が成り立つ．
C ∈ Rm×m : 正則

B = AT C,

T

R(A ) = R(B).

従って，m = rankA = rankA = rankB = rankB T より R(A) = R(B T ) = Rm も成り立つ．
ここで，AB ∈ Rm×m , BA ∈ Rn×n , m ≤ n だから，GMRES-LS 法 1（ minm b − ABz2 ）をTz∈R用いる方が GMRES-LS 法 2 を用いるより反復当たりの計算量が少なくてすむ．rankA = m な
らば AAT ，そして diag(AAT ) は正則なので，C の簡単な例として C := {diag(AAT )}−1 が考えられる．
3.

B の選び方R(A) = R(B T ), R(AT ) = R(B) を充たす以外に，収束を速めるには，B は AB ≈ Im またはBA ≈ In を充たすことが望ましい．
簡単な候補としては，先にも述べたように，m ≥ n = rankA の場合は C := {diag(AT A)}−1 とおき，rankA = m ≤ n の場合は C := {diag(AAT)}−1 とおくことがまず考えられる．
3.1 不完全 QR 分解
ここでは，m ≥ n = rankA（優決定）の場合に，もう少し複雑な前処理として，A = QR + E で与えられる不完全 QR 分解の適用を検討する．ただし，Q ∈ Rm×n は直交行列の近似，R ∈ Rn×n
は上三角行列，E は誤差行列とする．
従来はこの R を式（1.3）の前処理に用いて CG 法を適用していたが，ここでは B = R−1 QT とおいて，GMRES（k）-LS 法 1，2 を適用することを考える．
例えば GMRES（k）-LS 法 2 の場合，R−1 QT Ax = R−1 QT b に GMRES（k）法を適用する．ここ−1 T
で R Q A は CGLS 法に対する式（1.3）の Ã = R−T AT AR−1 よりも条件がよいと期待される．
3.2 IMGS（l）法
不完全 QR 分解の一つとして，修正 Gram-Schmidt 法（例えば Björck, 1996）の近似として行列 Q の現在の列ベクトルを過去の l 本の列ベクトルとのみ直交化させる，下記の IMGS
（l）法を
考える．
(1)

ai = ai(i = 1, . . . , n)for i = 1, 2, . . . , n

(i)

ai
rii
for j = i + 1, . . . , min(i + l, n)(i)rij = q Ti aj
(i)rii = ai 2 ,

(i+1)

aj
endfor
endforqi =(i)= aj − rij q i

GMRES 法による最小二乗問題の解法

341

上記の各ステップにおいて，各 q i は a1 , . . . , an の線形結合で，q 1 , . . . , q n  = a1 , . . . , an  が充たされるので，C̃ を正則な行列として，Q = AC̃ が成り立つ．従って，QT = C̃ T AT が成り立つ．
また，rii = 0 (i = 1, . . . , n) より，R = (rij ) も正則である．
従って，C を正則な行列として，B = R−1 QT = R−1 C̃ T AT = CAT が成り立つ．よって，IMGS（l）法は条件 R(A) = R(B T ) および R(AT ) = R(B) を充たす．
ここで，下記が成り立つ．
補題 3.1. IMGS（0）法は B = {diag(AT A)}−1 AT とおくことに相当する．
4.

数値実験

最後にフルランクな問題に対する数値実験結果を示す．
4.1 優決定問題
まず，優決定な最小二乗問題
min b − Ax2 ,x∈Rn

A ∈ Rm×n(m ≥ n)において，係数行列 A は MATLAB 6 のコマンド sprandn で生成した．その条件数や密度（非零要素の割合）は指定し，非零要素の値は正規分布に従う乱数で発生し，非零要素のパタンも乱数を用いて定めた．具体的には，小さめの行列として，m = 1, 000, n = 320，非零要素の密度が 4.9% の行列 RANDS2（条件数：2 × 102 ）および RANDS6（条件数：1 × 106 ）と，大きめの行列として，m = 10, 000, n = 1, 000，非零要素の密度が 1.5%で，条件数が表 1 のような行列
RANDL1–RANDL7 を用いた．
b についても同様に，各々m = 1, 000 および m = 10, 000 に対して，その各成分を正規分布に従う乱数で発生させ，b ∈/ R(A) となるようにした．
手法としては，下記の（1）–（5）を比較した．
（1）CGLS 法：正規方程式 AT Ax = AT b を CG 法で解く．
（2）正規方程式に Jennings and Ajiz（1984）による IMGS 法の前処理を施してから CGLS 法を適用する．
（3）正規方程式に提案した IMGS（l）法の前処理を施してから CGLS 法を適用する．
（4）IMGS（l）法により B = R−1 QT を求め，GMRES（k）-LS 法 1 を適用する．
（5）IMGS（l）法により B = R−1 QT を求め，GMRES（k）-LS 法 2 を適用する．
表 1.

テスト行列の条件数（優決定問題）．

342

統計数理 第 53 巻

第 2 号 2005

ただし，初期近似解は x0 = 0 とし，収束の判定には，r = b − Ax を残差として，AT r2 /AT b2< 10−6 を用いた．
各手法のプログラムは MATLAB 6 で書き，計算機は NEC の PC（AMD Athlon XP 2000+，
1.66GHz，736MB RAM）を用いた．
（2）の Jennings らの IMGS 法における閾値としては τ = 1（正規方程式の行対角項スケーリン
グに相当）が一番（計算時間の上で）速かったが，下記の実験では τ = 0.1 とおいた．Saad（1988）の IMGS 法も試みたが，Jennings らの方法よりさらに遅かった．
また，提案した IMGS（l）法において l の値の最適値を調べるために，RANDS2，RANDS6 に対して，（5）の IMGS（l）法を用いた GMRES-LS 法 2（リスタートなし）
で l の値を変化させたところ，
各々表 2，表 3 のように l = 0（正規方程式の行対角項スケーリングに相当）が一番速かったので下記の実験では l = 0 とおいた．これは，補題 3.1 にあるように，B = CAT , C = {diag(AT A)}−1
とおくことに相当する．
さらに，表 1 の行列のうち，RANDL2，RANDL3 に対して（IMGS（0）法を用いた）GMRES（k）-LS 法 2 において，リスタートの周期 k を変化させたところ，表 4 のようになった．
表 4 において，k はリスタート周期，iter，time は各々収束に要する反復数と計算時間（秒）を示す．また，∗ は各問題で一番速かったものを示す．ここでは，反復数，計算時間とも，リス
タートを行わない方が少なかったので，以下の実験でもリスタートなしの GMRES 法を用いた．
まず，図 1，2 に各々上記の RANDS2，RANDS6 に対する各手法の反復数 vs. AT r2 /AT b2のグラフを示す．悪条件問題 RANDS6 では（5）の GMRES 法 2 は（2）の前処理 CGLS 法や（4）
表 2.

IMGS（l）法を用いた GMRES-LS 法 2 での l の影響（RANDS2）．収束判定：
AT r2 /AT b2 < 10−6 ．

表 3.

IMGS（l）法を用いた GMRES-LS 法 2 での l の影響（RANDS6）．収束判定：
AT r2 /AT b2 < 10−6 ．

GMRES 法による最小二乗問題の解法

表 4.

GMRES（k ）-LS 法 2 でのリスタート周期 k の影響（優決定問題）．上段：リスタート
周期，中段：反復数，下段: 計算時間（秒），収束判定：AT r2 /AT b2 < 10−6 ．

図 1.

反復数 vs. AT r2 /AT b2（RANDS2）．

図 2.

反復数 vs. AT r2 /AT b2（RANDS6）．

343

344

統計数理 第 53 巻
表 5.

第 2 号 2005

各手法の比較（優決定問題）．上段：反復数，下段：計算時間（秒），収束判定：
AT r2 /AT b2 < 10−6 ．

の GMRES 法 1 よりもよい収束性を示しているのに注意されたい．
また，表 5 に表 1 の各問題に対する各手法の比較結果を示す．各ますの上段は収束に要した反復数を，下段は計算時間（秒）を示す．ただし，† はその反復数では収束しなかったことを示す．
表 5 から下記のことが言える．
条件数が 101 ∼ 105 の問題では（3）の IMGS
（0）前処理付きの CGLS 法が一番速かった．
（2）の
Jennings らの IMGS 前処理付きの CGLS 法は反復数は少ないが，計算時間がかかった．これは，前処理や前進後退代入の部分の影響と思われる．
条件数が 106 , 107 の条件の悪い問題では，（1）の CGLS 法，（2），（3）の前処理付き CGLS 法，（4）の IMGS（0）前処理付き GMRES-LS 法 1 は反復数，計算時間ともかかりすぎたのに対し，（5）の IMGS（0）前処理付きの GMRES-LS 法 2 は少ない反復数と計算時間で収束した．
（5）が（3）のIMGS（0）
前処理付きの CGLS 法より速く収束した理由としては，GMRES 法は Gram-Schmidtの直交化を陽に行うのに対し，CG 法は三項漸化式により直交化を行うため，特に悪条件問題
において GMRES 法は CG 法より丸め誤差の影響を受けにくいことが考えられる．
また，（5）の GMRES-LS 法 2 は（4）の GMRES-LS 法 1 より反復当たりの計算時間も少なくてすんでいる．これは，方法 2 は n × n 行列 BA に基づくクリロフ部分空間を用いているのに対し，方法 1 は m × m 行列 AB に基づくクリロフ部分空間を用いており，今の場合 n < m の優決定の最小二乗問題を扱っているためである．
4.2 劣決定問題
次に，劣決定な最小二乗問題
min b − Ax2 ,x∈Rn

A ∈ Rm×n(m < n)に対する数値実験結果を示す．係数行列 A は表 1 の各行列 RANDLn を転置したもので，以下では RANDLnT とよぶことにする．従って，m = 1, 000, n = 10, 000 で，非零要素の密度は

GMRES 法による最小二乗問題の解法

345

1.5%である．また，RANDLnT の条件数は対応する元の行列 RANDLn の条件数と等しい．bは x∗ = (1, . . . , 1)T に対して b = Ax∗ とおいた．
劣決定問題では rankA = m であれば，Ax = b となる x が（無数に）存在するので収束はr = b − Ax として r2 /b2 < 10−6 により判定した．初期近似解は x0 = 0 とした．
この問題を下記の（1）–（6）の手法で解いた．

（1）CGLS 法．
（2）C = {diag(AAT )}−1 により対称に前処理した CGLS 法．
（IMGS（0）法により対角項スケー
リングを施した CGLS 法．
）（3）GMRES-LS 法 1 で B = AT とおいたもの．
（4）GMRES-LS 法 1 で B = AT C, C = {diag(AAT )}−1 とおいたもの．
（5）GMRES-LS 法 2 で B = AT とおいたもの．
（6）GMRES-LS 法 2 で B = AT C, C = {diag(AAT )}−1 とおいたもの．
プログラム言語と計算機は優決定問題のものと同じである．
表 6 に各問題に対する各手法の比較結果を示す．各ますの上段は反復数を，下段は計算時間（秒）を示す．ただし，† はその反復数では収束しなかったことを示す．
表 6 より，反復数では（3）か（4）の GMRES-LS 法 1 が一番少ない反復数で収束している．計
算時間では，RANDL1T から RANDL6T の問題では（1）の CGLS 法か（2）の IMGS（0）-CGLS 法が一番速く収束しているが，条件の最も悪い RANDL7T では（3）の GMRES-LS 法 1 が一番速く収束している．また，前処理をしない（1）の CGLS 法と（3）の GMRES-LS 法 1 を比べると，条件数が 105 以上の問題では後者の方が速く収束している．最後に，GMRES-LS 法 1 と 2 を比べると，収束までの反復数は前者の方がやや少ない上に，劣決定問題では前者の方が扱う行列が小さく，反復当たりの計算量が少ないため，速く収束している．
表 6.

各手法の比較
（劣決定問題）．上段：反復数，下段：計算時間（秒），収束判定：r2 /b2 <10−6 ．

346

5.

統計数理 第 53 巻

第 2 号 2005

まとめ

GMRES 法を最小二乗問題に適用する方法として，もとの m × n の係数行列 A に対して，n × mの行列 B を用いて， minm b − ABz2 に GMRES 法を適用する方法 1 と， minn Bb − BAx2
z∈R

x∈R

に GMRES 法を適用する方法 2 を提案した．
また，それらの方法が任意の b に対して破綻することなく最小二乗解を与えるための行列 Bに関する十分条件を導いた．さらに，B の一例として不完全 QR 分解である IMGS（l）法を提案
した．
フルランクな優決定の最小二乗問題に対する数値実験によると，条件の悪い問題では IMGS（0）を用いた GMRES 法 2（正規方程式に行対角項スケーリングを施してから GMRES 法を適用することに相当する）は前処理付きの CGLS 法などの従来法よりも速く収束した．
また，フルランクな劣決定問題に対する数値実験でも，非常に条件の悪い問題では前処理なしの GMRES 法 1 が前処理付きの CGLS 法などの従来法よりも速く収束した．
謝

辞本研究に関して有益な示唆をいただいた，張 紹良先生，杉原厚吉先生，Yimin Wei 先生，
細田陽介先生，そして，多くの有益な指摘をいただいた査読者の方々に感謝いたします．また，本研究に関連する分野も含めた，幅広い分野において先駆的な研究をされてこられ，常日頃より有益な示唆をいただいている田邉國士先生に敬意と感謝を表します．本研究は文部科学省科学研究費補助金の助成を受けています．

参 考 文 献Bai, Z.-Z., Duﬀ, I. S. and Wathen, A. J.（2001）. A class of incomplete orthogonal factorization
methods. I: Methods and theories, BIT, 41（1）, 53–70.
Benzi, M. and Tuma, M.（2003）. A robust preconditioner with low memory requirements for large
sparse least squares problems, SIAM Journal on Scientific Computing, 25（2）, 499–512.
Björck, A.（1996）. Numerical Methods for Least Squares Problems, SIAM, Philadelphia.
Brown, P. N. and Walker, H. F.（1997）. GMRES on（nearly）singular systems, SIAM Journal on
Matrix Analysis and Applications, 18（1）, 37–51.
Calvetti, D., Lewis, B. and Reichel, L.（2000）. GMRES-type methods for inconsistent systems, Linear
Algebra and Its Applications, 316, 157–169.
謙（2003）. GMRES and GCR（k ）on singular systems, 第 32 回数値解析シンポジウム講演予稿集，2003 年 5 月，箱根，71–74.
速水 謙，伊藤徳史（2004）. GMRES 法を最小二乗問題に適用するための前処理行列が充たすべき性
質について，日本応用数理学会 2004 年度年会講演予稿集，212–213.
速水 謙，伊藤徳史（2005）. 前処理付き GMRES 法による最小二乗問題の解法，京都大学数理解析研
究所講究録 1441, 21 世紀における数値解析の新展開，2005 年 7 月，114–128.
伊藤徳史（2003）. GMRES 法の線形最小二乗問題への適用，修士論文，東京大学大学院情報理工学系
研究科数理情報学専攻．
伊藤徳史，速水 謙（2002）. GMRES（k ）法の線型最小二乗問題への適用，日本応用数理学会 2002 年度
年会予稿集，http://pjsiam.jstage.jst.go.jp/ja
伊藤徳史，速水 謙（2003）. 不完全 QR-GMRES（k ）法による線形最小二乗問題の解法，情報処理学会
第 65 回全国大会講演論文集，（1-117）–（1-118）．
Ito, T. and Hayami, K.（2004）. Preconditioned GMRES methods for least squares problems, NII

速水

GMRES 法による最小二乗問題の解法

347

Tech. Report, NII-2004-006E, 1–29, National Institute of Informatics.
伊藤 徳史，速水 謙（2004）. 前処理 GMRES 法による最小二乗問題の解法，2004 年度日本応用数理学会年会予稿集，210–211.
Jennings, A. and Ajiz, M. A.（1984）. Incomplete methods for solving AT Ax = b, SIAM Journal on
Scientific and Statistical Computing, 5, 978–987.
Meijerink, J. A. and van der Vorst, H.（1977）. An iterative method for linear systems of which thecoeﬃcient matrix is a symmetric M-matrix, Mathematics of Computation, 31（137）, 148–162.
ラオ，C. R.（1977）. 『統計的推測とその応用』，東京図書，東京．
Reichel, L. and Ye, Q.（2005）. Breakdown-free GMRES for singular systems, SIAM Journal on
Matrix Analysis and Applications, 26（4）, 1001–1021.
Saad, Y.（1988）. Preconditioning techniques for nonsymmetric and indeﬁnite linear systems, Journalof Computational and Applied Mathematics, 24, 89–105.
Saad, Y. and Schultz, M. H.（1986）. GMRES: A generalized minimal residual algorithm for solving（3）,
nonsymmetric linear systems, SIAM Journal on Scientific and Statistical Computing, 7
856–869.
Tanabe, K.（1974）. Characterization of linear stationary iterative processes for solving a singularsystem of linear equations, Numerische Mathematik, 22, 349–359.
Tanabe, K.（1975）. Neumann-type expansion of reﬂexive generalized inverses of a matrix and its
hyperpower iterative method, Linear Algebra and Its Applications, 10, 163–175.
張 紹良（1989）. 共役残差法の一般化，博士論文，筑波大学大学院博士課程工学研究科電子情報学専攻．
Zhang, S.-L. and Oyanagi, Y.（1990）. A necessary and suﬃcient convergence condition of orthomin（k ）methods for least squares problem with weight, Annals of the Institute of Statistical Mathematics, 42, 805–811.
Zhang, S.-L. and Oyanagi, Y.（1991）. Orthomin（k ）method for linear least squares problem, Journalof Information Processing, 14, 121–125.

348

Proceedings of the Institute of Statistical Mathematics Vol. 53, No. 2, 331–348(2005)

Solution of Least Squares Problems Using GMRES Methods
Ken Hayami1 and Tokushi Ito2
1 National Institute of Informatics

2 Business Design Laboratory Co., Ltd.

The most commonly used iterative method for solving large sparse least squares problems minn b − Ax2 , A ∈ Rm×n is the CGLS method, which applies the (preconditioned)x∈R

conjugate gradient method to the normal equation AT Ax = AT b.
This paper considers alternative methods by using an n × m matrix B and applying
the Generalized Minimal Residual (GMRES) method, which is a robust Krylov subspace
iterative method for solving systems of linear equations with nonsymmetric coeﬃcientmatrix, to minm b − ABz2 or minn Bb − BAx2 .
z∈R

x∈R

Next, we give a suﬃceint condition concerning B for the proposed methods to givea least squares solution without breakdown for arbitrary b, for over-determined, underdetermined and possibly rank-deﬁcient problems. Then, as an example for B, we propose
the IMGS(l) method, which is an incomplete QR decomposition.
Finally, we show by numercial experiments on full-rank over-determined and underdetermined problems that, for ill-conditioned problems, the proposed method using the
IMGS(0) method, which is equivalent to diagonal scaling, gives a least squares solution
faster than previous preconditioned CGLS methods.

Key words: Least squares problems, iterative method, Krylov subspace method, GMRES method,
CGLS method, preconditioning.