分解 (1)LU 連立一次方程式の解法

連立一次方程式の 解法 (1)

LU 分解

• 杉原,室田: 線形計算の数理, 岩波書店, 2009

• 斎藤: 数値解析入門, 東京大学出版会, 2012

• 久保田: 工学基礎 数値解析とその応用,数理工学 社, 2010

• 伊理: 線形代数汎論, 浅倉書店, 2009

• 伊理,藤野: 数値計算の常識, 共立出版, 1996

• 渡部: 連立1次方程式の基礎知識, 九州大学大型 計算機センター広報, Vol.28, No.4, pp.291–349, 1995.

• 数学を使って応用上の問題を解くということ は, 方程式を立てて解を求めることに相当. .

• 多くの場合, 少なくとも局所的には, 線形近 似が有効.

• 線形の数学モデルは非線形の数学モデルより 圧倒的に取り扱いやすい.

• 必要に応じて線形近似してから問題を解くと いうことがよく行われる.

• このような場合には, 最終的に線形方程式を 解けばもとの問題の(近似?)解が得られる.

• 今回および次回の講義では, 方程式に微分演 算が含まれない場合について考える.

• 微分が含まれる場合については第11回から 第14回で取り扱う.

• 変数が1個の線形方程式の解法について悩む ことは何もないので・・・

• はじめから多変数の場合を考える.

• 多変数の微分を含まない線形方程式を連立一 次方程式ともいう.

• Aをm行n列の行列, xをn次のベクトル, bをm次のベクトルとする.

• 数学的には, 連立一次方程式を解くとは, Ax =b

を満たすxをすべて求めることを意味する.

• 連立一次方程式の解: 一意解,不定解,不能解

• Scilabで連立一次方程式を解くには演算子\(環

境によっては￥記号)を使う.

• 行列AおよびbがScilabにおいて変数Aおよ びbに格納されているとき, ScilabでAx =b の解を求めるには,

x=A\bのようにする.

• 一意解の場合には, Ax = bの解とScilabの x=A\bは数値計算の誤差を除いて一致.

• 不定解の場合は, Scilabのx=A\bはAx = b を満たす解のひとつを返す.

• 不能解の場合, ScilabはkAx−bkが最小と なるxを近似解として返す.

一意解の例：

0 3

1

x₂ =

5 ^の解は(x₁, x₂) = (2/3,5/3)であり, Scilabで計算すると以下のようにな る.

-->A=[1 2;0 3];b=[4;5];x=A\b x =

0.6666667 1.6666667

不定解の例: 1 0 ¹

x₂ = 2の解は不定で, (x₁, x₂) = (2,0)はひとつの解であるが, Scilabで計算すると以下 のようになる.

-->A=[1 0];b=2;x=A\b x =

2.

0.

1 x =

3 ^{は解を持たないが}, Scilabは最小二乗近 似解を与える.

-->A=[1;1];b=[2;3];x=A\b x =

2.5

• Aをm行n列の行列, xをn次のベクトル, bをm次のベクトルとする.

• 数学的には, 連立一次方程式を解くとは, Ax =b

を満たすxをすべて求めることを意味する.

• 行列Aとベクトルbを結合して作った行列 B= (A|b)を,連立一次方程式Ax=bの拡 大係数行列という.

• rankA<rankBなら不能解.

• rankA= rankB = dimxなら一意解

• rankA= rankB <dimxなら不定解.

• 行列に行基本変形を施すことで,以下のよう な階段行列が得られる.











0 · · · 0 1 ∗ · · ·

0 · · · 0 1 ∗ · · ·

0 · · · 0 1 ∗ · · ·

. . .











• 階段行列の特徴は以下の通り: 1の左はすべ て零, 1の右側は任意, 1を含まない行はすべ て零だけ, 行列全体を見ると1の系列は右斜 め下に進む(真下不可)

• 第1行左端に零があるかないかは行列によっ て変わる.

• 一意解に対応する拡大係数行列を階段行列に変形 すると(今回の講義ではこの場合を取り扱う)











1 ∗ · · · ∗

0 1 ∗ · · · ∗

... . .. ... ∗ ∗

0 · · · 0 1 ∗











• 拡大係数行列を階段行列に変換する手順は行 基本変形.

• 行基本変形は基本行列を拡大係数行列に左か ら掛けることに対応.

• 行基本変形によって,係数行列Aが上三角行 列(後述)に変形されていることになる.

• 上三角行列とはこういう行列:













∗ · · · ∗ 0 . .. ... ... . .. ... ...

0 . . . 0 ∗













ただし∗は任意の数(零でもよい).

• 連立一次方程式の解法は, 大別すると, 直接 法, 反復法の2種類に分類される.

• 直接法は数値計算の誤差がない場合に有限回 の演算で解を与える手法.

• 反復法は繰り返しによって近似解の系列を生 成する方法.

• 上記以外に, 共役勾配法と呼ばれる方法があ る.

• 消去法は直接法の一種. 今回の講義では直接 法を取り扱う.

• 以下では, 行列Aはn行n列の正則行列と する.

• もっとも素朴なガウスの消去法は,行列Aが正方 かつ正則で,行列Aが行の入れ換えなしに階段行 列に変形できる場合に相当.

• この場合に相当する, 行列のLU分解と呼ばれる ものを,これから導出する.

• A⁽¹⁾ =Aとし,以下,繰り返し計算によって, こ れをA⁽²⁾,A⁽³⁾,. . . のように変形してゆく.

• A=A⁽¹⁾の各成分を以下のように書く.

A⁽¹⁾ =













a⁽¹⁾₁₁ a⁽¹⁾₁₂ · · · a⁽¹⁾_1n a⁽¹⁾₂₁ a⁽¹⁾₂₂ · · · a⁽¹⁾_2n

... ...

a⁽¹⁾_n1 a⁽¹⁾_n2 · · · a⁽¹⁾nn













• a⁽¹⁾₁₁ 6= 0と仮定し,

l₁ =















 1

a⁽¹⁾₂₁ a⁽¹⁾₁₁

...

a⁽¹⁾n1

a⁽¹⁾₁₁

















, u₁ =

a⁽¹⁾₁₁ a⁽¹⁾₁₂ · · · a⁽¹⁾_1n とす

ると・・・

• A⁽¹⁾ =l₁u^T₁ +A⁽²⁾,

A⁽²⁾ =













0 0 · · · 0

0 a⁽²⁾₂₂ · · · a⁽²⁾_2n

... ...

0 a⁽²⁾_n2 · · · a⁽²⁾nn













という形にな

る(各成分の式は略).

• a⁽²⁾₂₂ 6= 0と仮定し,以下のようにおくと:

l₂=



















 0 1

a⁽¹⁾₃₂ a⁽²⁾₂₂

...

a⁽¹⁾n2

a⁽²⁾₂₂





















,u₂=

0 a⁽²⁾₂₂ a⁽²⁾₂₃ · · · a⁽²⁾_2n

• A⁽²⁾ =l₂u^T₂ +A⁽³⁾,

A⁽³⁾ =













0 0 · · · 0 0 0 · · · 0 0 0 a⁽³⁾₂₂ · · ·

... ...













という形になる (各成分の式は略).

• 同様にして一回計算するごとに行列の左およ び上側の零列と零行が1ずつ増えるから・・・

• A⁽ⁿ⁾ = l_nu^T_n +A⁽ⁿ⁺¹⁾ とすると, A⁽ⁿ⁺¹⁾ = 0(零行列)である.

• 以上をまとめるとA=l₁u^T₁ +· · ·+l_nu^T_nと なるが・・・

• L=

l₁ · · · l_n

,U =





 u^T₁

... u^T_n





とおくと・・・

• A = LU となる. これを行列AのLU分 解という. Lが下三角行列(後述), U が上三 角行列になっていることに注意.

• 下三角行列とはこういう行列:













∗ 0 · · · 0 ... . .. ... ...

... . .. 0

∗ · · · ∗













ただし∗は任意の数(零でもよい).

• LU分解を用いて連立一次方程式Ax=LU x= bを解くには, Ly =b, U x =yという2個 の連立一次方程式を順に解けばよい.

• A = LU というLU分解が得られていると き, さらに行列U を対角行列Dと対角要素 が1の上三角行列U^′の積であらわし, A = LDU^′と書き直すことがある. これをLDU 分解という.

L₁D₁U₁ = L₂D₂U₂であったと仮定すると, L⁻¹₂ L₁D₁ = D₂U₂U⁻¹₁ であるが, L⁻¹₂ L₁が対角要素が1の下三角行列, U₂U⁻¹₁ が対角要素が1の上三角行列であることに注意すると, まずD₁=D₂が得られ,続いてL⁻¹₂ L₁D₁=D₂U₂U⁻¹₁ の左 辺が下三角行列,右辺が上三角行列であることから,L⁻¹₂ L₁= I_n,U₂U⁻¹₁ =I_nとなり,よってL₂ =L₁,U₂=U₁となる からである.

• Aが対称行列でLDU分解できるとき, A = LDU とすると, A = A^T とLDU分解の一 意性から, L=U^T,U =L^T が導かれる. し たがって, A = LDL^T と書ける. これを対 称行列AのLDL^T 分解と呼ぶ.

• AがLDU分解できる正定対称行列である場 合には,Dの各要素は正だから,Dの対角要 素の正の平方根を対角要素とする行列を G とすると, D = GG^T であり, したがって A = LGG^TL^T と書ける. C = LGとお くと, A = CC^T である. この表現を, Aの Cholesky分解と呼ぶ.

• Gaussの消去法をLU分解から導く.

• ベクトルl_kの第k成分が1であったことを思い 出し, 第k+ 1成分以降をまとめたn−k−1次 のベクトルを¯l_kと書くことにする.

• L₁ =

1 0

−¯l₁ In−1

とおく(In−1はn−1次の 単位行列, 0はその部分が零であることを示す)

• A⁽¹⁾ =l₁u^T₁+A⁽²⁾の両辺にL_{1を左から乗じ},L₁ の構造を利用して整理すると(詳細は略),L₁A⁽¹⁾ = e₁u^T₁ +A⁽²⁾となる(ただしe₁は第1番目の単位 ベクトル).

• 次に,A⁽²⁾=l₂u^T₂ +A⁽²⁾を考える.

• L₂ =





1 0 0

0 1 0

0 −¯l₂ I_{n 2}



_とおく.

• この行列の構造に注意して計算すると, L₂A⁽²⁾ = e₂u^T₂ +A⁽³⁾となる.

• L₂e₁ =e₁ であることに注意すると, L₂L₁A=e₁u^T₁ +e₂u^T₂ +A⁽³⁾となる.

• 以下同様にして,

L_n· · ·L₁A=e₁u^T₁ +· · ·e_nu^T_n =U となる(A⁽ⁿ⁺¹⁾= に注意).

• 以上によって得られた式のL_n· · ·L₁は基本 行列の積であり,右辺は階段行列になってい る. 階段行列を求める手順がGaussの消去法 だったから, LU分解からGaussの消去法が 導かれたことになる.

• 上述のように階段行列を求める手順を前進消 去という.

• 行交換が必要ない場合には, Gaussの消去法 も, LU分解も, l₁, . . . ,l_n, u₁, . . . ,u_nを求め ることに相当するので, Gaussの消去法とLU 分解は本質的に同じ.

• Ax=bは, ¯b=L_n· · ·L₁bとおくと,

L_n· · ·L₁Ax = U x = ¯b と変形されるから, U x= ¯bという連立一次方程式を解くことに より, 解xが求められる.

• 具体的には, ¯bの各成分を¯b1, . . . ,¯bnとし, U の第(i, j)成分をuij とし, U が正則な上三 角行列であったことに注意すると,まずxn =

¯bn/unnが得られ,次にun−1,n−1xn−1+un−1,nxn=

¯bn−1にこれを代入してxn−1が得られ,という ふうに,逐次的に解xの全成分が得られる. こ の操作を後退代入という.

• 行列Aが正則であっても, a11 = 0であると いうことはあり得る.

• a11 6= 0であっても,絶対値が零に近い場合に は,a11を使ってLU分解あるいはGaussの消 去法によって連立一次方程式の解を求めると, 数値計算の誤差が大きくなる可能性がある.

• 以下では, 仮に, 零あるいは零に近い要素を

「条件が悪い」と呼び,そうでない要素を「条 件が良い」と呼ぶ.

• 計算不能あるいは数値的な条件悪化を防ぐに は, 行を入れ換えて, 条件がよいak1,1が行列 の一番上に来るようにすればよい.

• これは, 見方を変えると,a11のかわりにak1,1

に着目して, LU分解あるいはGaussの消去 法を遂行していることになる.

• 第kステップについても同様に,数値的な条 件が良いajk,k に着目してLU分解あるいは

Gaussの消去法を遂行してゆく.

• LU分解あるいはGaussの消去法の各ステッ プで着目している条件が良い要素の番号(jk, k) のことをピボットあるいは枢軸と呼ぶ.

• 対応するa^(k)_jk,kのことをピボット要素あるい は枢軸要素と呼ぶ.

• 数値的な条件が良いように適切にピボットを 選ぶ操作のことを, ピボット選択あるいは枢 軸選択と呼ぶ.

• ピボット選択の選択法のひとつは,{a^(k)_j,k :j = k, k+ 1, . . . , n}の中で絶対値が最大の要素の 添字を選ぶ方法(行交換によるピボット操作).

• 列交換によるピボット操作と呼ばれる第k列 以降の列について同様の操作をする手法や, 完全ピボット操作と呼ばれる第(p, q)要素(た だしp, q ≥ k)すべての中から絶対値最大の ものを選ぶ手法もある.

• 列交換によるピボット操作や完全ピボット操 作では,変数の順番もこれに対応して入れ換 える必要がある.

• 数値計算の誤差はベクトルbの成分の大きさ にも依存するので, ピボット操作だけで数値 計算の誤差の低減が保証されるとは限らない.

• 方程式Ax =bを解くためにピボット選択付

きのGaussの消去法を使うことは,P をそれ

に対応する置換行列としたとき, P AをLU 分解することに相当する.

• ScilabにおいてLU分解を求める関数はluで ある.

• [L,U,P]=lu(A)とすることで,P A=LUと なる行列L, U, P を求めることができる.

• 実行例は次ページの通り.

0.1428571 1. 0.

0.5714286 0.5 1.

-->U U =

7. 8. 9.

0. 0.8571429 1.7142857

0. 0. 1.110D-16

-->P P =

0. 0. 1.

1. 0. 0.

0. 1. 0.

7. 8. 9.

1. 2. 3.

4. 5. 6.

-->L*U ans =

7. 8. 9.

1. 2. 3.

4. 5. 6.

• Gaussの消去法の終了後に,さらに列基本変 形を施して,行列U を単位行列に変換する方 法もある. これをGauss-Jordan法という.

• 計算量という観点から言うとGauss-Jordan 法にはあまりメリットはないが,並列計算機 には適しているという指摘もある.

• 要素の大部分が零の行列を疎行列という.

• 応用であらわれる大規模行列は多くの場合疎 行列である.

• 疎行列で零要素をメモリに格納することは無 駄であるので,必要な要素だけをメモリに記 憶する方法が工夫されている.

• 疎行列に対する演算は,行列の疎性を破壊し ないことが望ましい.

• Scilabで疎行列を扱うための組み込み関数は

sparse.