連立一次方程式 \(Ax = b\) を解くとき、あらかじめ係数行列 \(A\) を下三角行列 \(L\) (lower) と上三角行列 \(U\) (Upper) に分解しておくと、簡単に解くことができます。\(Ax = b\) は \(LUx = b\) と表すことができるので、\(Ux = y\) とおくと \(Ly = b\) を解くことで \(y\) を求めることができます。次に、求めた \(y\) を使って \(Ux = y\) を解けば、\(Ax = b\) の解 \(x\) を求めることができます。簡単な例を示しましょう。
連立方程式
x + y + z = 10
2x + 4y + 6z = 38
2x + 4z = 14
[ 1. 1. 1.; [ 1. 0. 0.; [1. 1. 1.;
2. 4. 6.; = 2. 1. 0.; * 0. 2. 4.;
2. 0. 4.] 2. -1. 1. ] 0. 0. 6. ]
(1) LU 分解
[ 1. 0. 0. 10.; y1 = 10
2. 1. 0. 38.; y2 = 38 - 2 * 10 = 18
2. -1. 1. 14. ] y3 = 14 - 2 * 10 + 18 = 12
(2) 前進代入
[ 1. 1. 1. 10.; x1 = 10 - x2 - x3 = 3
0. 2. 4. 18.; x2 = (18 - 4 * x3) / 2 = 5
0. 0. 6. 12. ] x3 = 12 / 6 = 2
(3) 後退代入
LU 分解はガウスの消去法と同じアルゴリズムで求めることができます。U は前進消去で得た行列で、L の要素は前進消去するとき係数を 0 にするための倍率になります。実際に L * U を計算すると元の行列に戻ります。Ly = b を解く場合、L は下三角行列なので前進代入で y を求めることができます。次に Ux = y を解く場合、U は上三角行列なので後退代入で x を求めることができます。
ガウスの消去法の計算量は \(O(N^3)\) ですが、前進代入と後進代入の計算量は \(O(N^2)\) です。同じ係数行列を何度も使う場合、あらかじめ係数行列を LU 分解しておけば、ガウスの消去法よりも効率的に連立方程式を解くことができるわけです。
Julia には LU 分解を行う関数 lu() が用意されているので、実用的にはそれを使えばいいのですが、Julia とアルゴリズムの学習ということで、あえてプログラムを作ってみましょう。
実際には \(A\) を \(L\) と \(U\) の二つに分けるのではなく、\(U\) の 0 の部分に \(L\) の要素をセットした行列 \(A'\) を作ることにします。プログラムは次のようになります。
リスト : LU 分解
function mylu(xs::Matrix{Float64})
n = size(xs, 1)
zs = copy(xs)
for i = 1 : n
for j = i + 1 : n
temp = zs[j, i] / zs[i, i]
for k = i + 1 : n
zs[j, k] -= temp * zs[i, k]
end
zs[j, i] = temp
end
end
zs
end
function mylu_solver(xs::Matrix{Float64}, ys::Vector{Float64})
n = size(xs, 1)
zs = copy(ys)
# 前進代入
for i = 2 : n, j = 1 : i - 1
zs[i] -= xs[i, j] * zs[j]
end
# 後退代入
for i = n : -1 : 1
for j = i + 1 : n
zs[i] -= xs[i, j] * zs[j]
end
zs[i] /= xs[i, i]
end
zs
end
関数 lu() は引数の係数行列 xs を LU 分解します。アルゴリズムはガウスの消去法と同じですが、係数を 0 にするための倍率 temp を zs[j, i] にセットするところが異なります。
関数 lu_solver() は LU 分解した行列 xs と右辺の定数を格納したベクトル ys を受け取り、その解を求めます。前進代入は後退代入とは逆に前から順番に値を求めていきます。行列 L の対角成分 L[i, i] は 1 なので、割り算する必要はありません。後退代入はガウスの消去法のプログラムと同じです。
それでは実際に試してみましょう。以下に簡単なテストと実行結果を示します。
リスト : 簡単なテスト
a1 = [1. 1.; 2. 4.]
a2 = [1. 1. 1.; 2. 4. 6.; 2. 0. 4.]
a3 = [1. 1. 1. 1.; -1. 1. -1. 1.; 8. 4. 2. 1.; -8. 4. -2. 1.]
a4 = [1. -1. 1. -1. 1.; 12. -6. 2. 0. 0.; 1. 1. 1. 1. 1.; 12. 6. 2. 0. 0.; 4. 3. 2. 1. 0.]
b1 = [100., 272.]
b2 = [10., 38., 14.]
b3 = [-5., -7., -31., -35.]
b4 = [1., 0., 8., 0., 1.]
for (a, b) = zip([a1, a2, a3, a4], [b1, b2, b3, b4])
xs = mylu(a)
println(xs)
println(mylu_solver(xs, b))
println("--------")
end
[1.0 1.0; 2.0 2.0] [64.0, 36.0] -------- [1.0 1.0 1.0; 2.0 2.0 4.0; 2.0 -1.0 6.0] [3.0, 5.0, 2.0] -------- [ 1.0 1.0 1.0 1.0; -1.0 2.0 0.0 2.0; 8.0 -2.0 -6.0 -3.0; -8.0 6.0 -1.0 -6.0] [0.0, -9.0, 1.0, 3.0] -------- [ 1.0 -1.0 1.0 -1.0 1.0; 12.0 6.0 -10.0 12.0 -12.0; 1.0 0.333333 3.33333 -2.0 4.0; 12.0 3.0 6.0 -12.0 -3.55271e-15; 4.0 1.16667 2.9 0.266667 -1.6] [0.3125, 0.0, -1.875, 3.5, 6.0625] --------
正常に動作していますね。
LU 分解にピボット選択を適用することも簡単です。次のリストを見てください。
リスト : LU 分解 (ピボット選択バージョン)
# ピボット選択
function select_pivot(xs::Matrix{Float64}, idx::Vector{Int}, i::Int)
_, k = findmax(abs.(xs[i:end, i]))
k += i - 1
if k != i
temp = xs[i, :]
xs[i, :] = xs[k, :]
xs[k, :] = temp
idx[i], idx[k] = idx[k], idx[i]
end
end
function mylu_pv(xs::Matrix{Float64})
n = size(xs, 1)
zs = copy(xs)
idx = collect(1 : n)
for i = 1 : n
select_pivot(zs, idx, i)
for j = i + 1 : n
temp = zs[j, i] / zs[i, i]
for k = i + 1 : n
zs[j, k] -= temp * zs[i, k]
end
zs[j, i] = temp
end
end
zs, idx
end
function mylu_solver_pv(xs::Matrix{Float64}, idx::Vector{Int}, ys::Vector{Float64})
mylu_solver(xs, ys[idx])
end
select_pivot() の引数 idx に交換した行の情報をセットします。係数行列の大きさを n とすると、idx の初期値は 0 から n - 1 までの整数を格納したベクトルになります。i 行と k 行を交換したときは idx[i] と idx[k] の値も交換します。関数 mylu_pv() はピボット選択しながら LU 分解を行い、L と U を格納した行列 zs とベクトル idx を返します。
関数 mylu_solver_pv() は、lu_pv() が返した行列 xs とインデックスリスト idx を受け取り、引数 ys のベクトルの要素を idx に合わせて並べ替えてから mylu_solver() を呼び出すだけです。
それでは実際に試してみましょう。以下に簡単なテストと実行結果を示します。
リスト : 簡単なテスト
a1 = [1. 1.; 2. 4.]
a2 = [1. 1. 1.; 2. 4. 6.; 2. 0. 4.]
a3 = [1. 1. 1. 1.; -1. 1. -1. 1.; 8. 4. 2. 1.; -8. 4. -2. 1.]
a4 = [1. -1. 1. -1. 1.; 12. -6. 2. 0. 0.; 1. 1. 1. 1. 1.; 12. 6. 2. 0. 0.; 4. 3. 2. 1. 0.]
a5 = [0. 2. 4.; 1. 1. 1.; 4. 2. 6.]
b1 = [100., 272.]
b2 = [10., 38., 14.]
b3 = [-5., -7., -31., -35.]
b4 = [1., 0., 8., 0., 1.]
b5 = [14., 10., 38.]
for (a, b) = zip([a1, a2, a3, a4, a5], [b1, b2, b3, b4, b5])
xs, idx = mylu_pv(a)
println(xs)
println(idx)
println(mylu_solver_pv(xs, idx, b))
println("--------")
end
[2.0 4.0; 0.5 -1.0] [2, 1] [64.0, 36.0] -------- [2.0 4.0 6.0; 1.0 -4.0 -2.0; 0.5 0.25 -1.5] [2, 3, 1] [3.0, 5.0, 2.0] -------- [ 8.0 4.0 2.0 1.0; -1.0 8.0 0.0 2.0; 0.125 0.0625 0.75 0.75; -0.125 0.1875 -1.0 1.5] [3, 4, 1, 2] [0.0, -9.0, 1.0, 3.0] -------- [12.0 -6.0 2.0 0.0 0.0; 1.0 12.0 0.0 0.0 0.0; 0.333333 0.416667 1.33333 1.0 0.0; 0.0833333 -0.0416667 0.625 -1.625 1.0; 0.0833333 0.125 0.625 -0.230769 1.23077] [2, 4, 5, 1, 3] [0.3125, 0.0, -1.875, 3.5, 6.0625] -------- [4.0 2.0 6.0; 0.0 2.0 4.0; 0.25 0.25 -1.5] [3, 1, 2] [5.0, 3.0, 2.0] --------
正常に動作していますね。
ここで、Julia の関数 lu() について簡単に説明しておきましょう。lu(A) は行列 \(A\) を LU 分解します。返り値を F とすると、F.L に下三角行列、F.U に上三角行列、F.p にピボット選択で交換した行の位置を表すベクトルが格納されています。次の例を見てください。
julia> a1 = [1. 1.; 2. 4.]
2×2 Matrix{Float64}:
1.0 1.0
2.0 4.0
julia> lua1 = lu(a1)
LU{Float64,Matrix{Float64}}
L factor:
2×2 Matrix{Float64}:
1.0 0.0
0.5 1.0
U factor:
2×2 Matrix{Float64}:
2.0 4.0
0.0 -1.0
julia> lua1.L
2×2 Matrix{Float64}:
1.0 0.0
0.5 1.0
julia> lua1.U
2×2 Matrix{Float64}:
2.0 4.0
0.0 -1.0
julia> lua1.p
2-element Vector{Int64}:
2
1
julia> (lua1.L * lua1.U)[lua1.p, :]
2×2 Matrix{Float64}:
1.0 1.0
2.0 4.0
行列 a1 を LU 分解して変数 lua1 にセットします。lua1.L と lua1.U を乗算して行を交換すると、元の配列 a1 に戻ります。そして、lua1 に対して演算子 \ や関数 det(), inv() を適用すると、a1 に適用したことと同じ結果を得ることができます。
julia> lua1 \ [100., 272.]
2-element Vector{Float64}:
64.0
36.0
julia> a1 \ [100., 272.]
2-element Vector{Float64}:
64.0
36.0
julia> det(lua1)
2.0
julia> det(a1)
2.0
julia> inv(lua1)
2×2 Matrix{Float64}:
2.0 -0.5
-1.0 0.5
julia> inv(a1)
2×2 Matrix{Float64}:
2.0 -0.5
-1.0 0.5
詳細は Julia のマニュアル Linear Algebra をお読みくださいませ。
Julia には行列 \(A\) の行列式 (determinant) を求める関数 det(A) が用意されていますが、\(A\) を LU 分解するときに計算することもできます。次に示す行列式と三角行列の性質を使います。以下では \(A\) の行列式を \(|A|\) と表しています。
ここでは数学的な証明は行わずに簡単な例を示すだけにとどめます。
(1)
[[a, b], @ [[e, f], = [[ae + bg, af + bh],
[c, d]] [g, h]] [ce + dg, cf + dh]]
右辺の行列式の乗算 (ad - bc) * (eh - fg) = adeh - adfg - bceh + bcfg
左辺の行列式 (ae + bg)(cf + dh) - (af + bh)(ce + dg)
= acef + adeh + bcfg + bdgh - acef - adfg - bceh - bdgh
^^^^ ++++ ^^^^ ++++
= adeh - adfg - bceh + bcfg
(2) [[a, b], [[c, d], [c, d]] [a, b]] ad - bc cb - da = - (ad - bc)
(3) [[a, b], [[a, 0], [0, d]] [c, d]] ad - b0 = ad ad - 0c = ad [[a, b, c], [d, e, f], => aei + bfg + chd - ceg - bdi - ahf [g, h, i]], 上三角行列は d, h, g が 0 なので行列式は aei になる 下三角行列は b, c, f が 0 なので行列式は aei になる
前回作成したプログラムでは行列 \(A\) を LU に分解したとき、\(L\) の対角成分が 1 になるので、\(U\) の対角成分の積が \(A\) の行列式になります。ピボット選択を行う場合、行列式の符号を記憶しておいて、行を交換したときは符号を反転すればいいでしょう。
プログラムは次のようになります。
リスト : 行列式の計算
# ピボット選択
function select_pivot(xs::Matrix{Float64}, idx::Vector{Int}, i::Int)
_, k = findmax(abs.(xs[i:end, i]))
k += i - 1
if k != i
temp = xs[i, :]
xs[i, :] = xs[k, :]
xs[k, :] = temp
idx[i], idx[k] = idx[k], idx[i]
-1
else
1
end
end
# LU 分解
function mylu_pv(xs::Matrix{Float64})
n = size(xs, 1)
zs = copy(xs)
idx = collect(1 : n)
d = 1
for i = 1 : n
d *= select_pivot(zs, idx, i)
for j = i + 1 : n
if zs[i, i] == 0 break end
temp = zs[j, i] / zs[i, i]
for k = i + 1 : n
zs[j, k] -= temp * zs[i, k]
end
zs[j, i] = temp
end
end
d * prod(diag(zs)), zs, idx
end
select_pivot() は行を交換したら -1 を、交換しなければ 1 を返すように変更します。mylu_pv() は行列式、LU 分解した行列、行の交換を記録したベクトルを返します。行列式の符号は変数 d にセットし、select_pivot() の返り値と乗算します。これで行を交換したら符号を反転することができます。対角成分 zs[i, i] が 0 ならば LU 分解できないので、break でループを脱出します。行列式の計算は関数 diag() で zs の対角成分を取り出し、prod() で要素の乗算を求めるだけです。対角成分に 0 の要素があるならば行列式は 0 になります。
それでは実際に試してみましょう。
リスト : 簡単なテスト
a1 = [1. 1.; 2. 4.]
a2 = [1. 1. 1.; 2. 4. 6.; 2. 0. 4.]
a3 = [1. 1. 1. 1.; -1. 1. -1. 1.; 8. 4. 2. 1.; -8. 4. -2. 1.]
a4 = [1. -1. 1. -1. 1.; 12. -6. 2. 0. 0.; 1. 1. 1. 1. 1.; 12. 6. 2. 0. 0.; 4. 3. 2. 1. 0.]
a5 = [0. 2. 4.; 1. 1. 1.; 4. 2. 6.]
a6 = [2. 4. 2. 2.; 4. 10. 3. 3.; 2. 6. 1. 1.; 3. 7. 1. 4.]
for a = [a1, a2, a3, a4, a5, a6]
d, xs, idx = mylu_pv(a)
println(d, " ", det(a))
println(xs)
println(idx)
println("--------")
end
2.0 2.0 [2.0 4.0; 0.5 -1.0] [2, 1] -------- 12.0 12.0 [2.0 4.0 6.0; 1.0 -4.0 -2.0; 0.5 0.25 -1.5] [2, 3, 1] -------- 72.0 72.0 [ 8.0 4.0 2.0 1.0; -1.0 8.0 0.0 2.0; 0.125 0.0625 0.75 0.75; -0.125 0.1875 -1.0 1.5] [3, 4, 1, 2] -------- 384.0000000000001 384.0000000000001 [12.0 -6.0 2.0 0.0 0.0; 1.0 12.0 0.0 0.0 0.0; 0.333333 0.416667 1.33333 1.0 0.0; 0.0833333 -0.0416667 0.625 -1.625 1.0; 0.0833333 0.125 0.625 -0.230769 1.23077] [2, 4, 5, 1, 3] -------- -12.0 -12.0 [4.0 2.0 6.0; 0.0 2.0 4.0; 0.25 0.25 -1.5] [3, 1, 2] -------- 0.0 0.0 [4.0 10.0 3.0 3.0; 0.5 -1.0 0.5 0.5; 0.75 0.5 -1.5 1.5; 0.5 -1.0 -0.0 0.0] [2, 1, 4, 3] --------
最後は行列式が 0 になる場合です。この場合、解を一意的に定めることはできません。解は無数にあるか存在しないかのどちらかになります。
連立一次方程式の解の存在条件は行列の「階数 (rank)」を使って判定することができます。階数の定義はいくつか方法がありますが、ここでは行列の基本変形を使うことにしましょう。以下に示す操作を行の基本変形といいます。
同様に列の基本変形を定義することができます。これらの基本変形は、ガウスの消去法やガウス・ジョルダン法で用いる操作と同じです。一般に、任意の行列は基本変形を何回か適用すると「階段行列」に変形することが可能です。階段行列は左下半分に 0 が階段状に並んだ行列のことです。
階段行列に変形したあと、零ベクトル (要素がすべて 0 のベクトル) を除いた行ベクトルの本数が階数になります。上図でいうと、(1) と (2) の行列には零ベクトルはないので階数は 3 になりますが、(3) の行列は零ベクトルが一つあるので階数は 2 になります。
[[a, b, c, d], [[a, b, c, d], [[a, b, c, d],
[0, e, f, g], [0, e, f, g], [0, e, f, g],
[0, 0, h, i]] [0, 0, 0, h]] [0, 0, 0, 0]]
(1) rank = 3 (2) rank = 3 (3) rank = 2
図 : 階段行列
ここで、上図の行列が拡大係数行列を表しているとしましょう。(1) の場合、係数行列と拡大係数行列の階数はどちらも 3 になり、拡大係数行列の行数と一致します。この場合、連立一次方程式の解は一意的に定まります。(2) の場合、3 行 3 列の係数行列としてみると、その階数は 2 になり、拡大係数行列の階数よりも少なくなります。この場合は解がありません。係数がすべて 0 なのに定数が残っているので、連立方程式として矛盾しているわけです。
(3) の場合、係数行列と拡大係数行列の階数はともに 2 になりますが、拡大係数行列の行数とは一致しません。このとき、零ベクトルに対応する変数に適当な値を代入すると、残りの変数を決定することができます。つまり、連立一次方程式の解は存在するが、一意的には定めることができないわけです。まとめると次のようになります。
簡単な例を示しましょう。
(1) 2a + 4b + 2c + 2d = 8 4a + 10b + 3c + 3d = 17 2a + 6b + c + d = 9 3a + 7b + c + 4d = 11 行の基本変形 [ 4. 10. 3. 3. 17. ; 0. -1. 0.5 0.5 -0.5 ; 0. 1. -0.5 -0.5 0.5 ; 0. -0.5 -1.25 1.75 -1.75 ] [ 4. 10. 3. 3. 17. ; 0. -1. 0.5 0.5 -0.5; 0. 0. 0. 0. 0. ; 0. 0. -1.5 1.5 -1.5 ] [ 4. 10. 3. 3. 17. ; 0. -1. 0.5 0.5 -0.5; 0. 0. -1.5 1.5 -1.5; 0. 0. 0. 0. 0. ] 係数行列と拡大係数行列の階数は 3 d = k とおくと c = (-1.5 - 1.5d) / -1.5 = k + 1 b = (-0.5 - 0.5c - 0.5d) / -1 = k + 1 a = (17 - 10b - 3c - 3d) / 4 = -4k + 1 k は任意の定数なので解は無数にある
(2) 2a + 4b + 2c + 2d = 8 4a + 10b + 3c + 3d = 17 2a + 6b + c + d = 10 3a + 7b + c + 4d = 11 行の基本変形 [ 4. 10. 3. 3. 17. ; 0. -1. 0.5 0.5 -0.5 ; 0. 1. -0.5 -0.5 1.5 ; 0. -0.5 -1.25 1.75 -1.75 ] [ 4. 10. 3. 3. 17. ; 0. -1. 0.5 0.5 -0.5; 0. 0. 0. 0. 1. ; 0. 0. -1.5 1.5 -1.5 ] [ 4. 10. 3. 3. 17. ; 0. -1. 0.5 0.5 -0.5; 0. 0. -1.5 1.5 -1.5; 0. 0. 0. 0. 1. ] 係数行列の階数は 3 で拡大係数行列の階数は 4 なので解はない
(3) 2a + 4b + 2c + 2d = 8 4a + 10b + 3c + 3d = 17 2a + 6b + c + 2d = 9 3a + 7b + c + 4d = 11 行の基本変形 [ 4. 10. 3. 3. 17. ; 0. -1. 0.5 0.5 -0.5 ; 0. 1. -0.5 0.5 0.5 ; 0. -0.5 -1.25 1.75 -1.75 ] [ 4. 10. 3. 3. 17. ; 0. -1. 0.5 0.5 -0.5; 0. 0. 0. 1. 0. ; 0. 0. -1.5 1.5 -1.5 ] [ 4. 10. 3. 3. 17. ; 0. -1. 0.5 0.5 -0.5; 0. 0. -1.5 1.5 -1.5; 0. 0. 0. 1. 0. ] 係数行列と拡大係数行列の階数は 4 d = 0 / 1 = 0 c = (-1.5 - 1.5*0) / -1.5 = 1 b = (-0.5 - 0.5*1 - 0.5*0) / -1 = 1 a = (17 - 10*1 - 3*1 - 3*0) / 4 = 1 解は一意的に定まる
拙作のページで説明した「ガウスの消去法」は、拡大係数行列を階段行列に変形することで解を求めています。このような方法を「直接法 (Direct Method)」といいます。これに対し、適当な初期解から始めて、繰り返し計算することで真の解に収束させていく方法が考えられます。これを「反復法 (Iterative Method)」といいます。
直接法は厳密解を求めることができますが、係数行列が大きくなると解くのに時間がかかるようになります。厳密解を求めるのが難しい場合、反復法を使うと現実的な時間で近似解を求めることが可能です。今回は基本的な反復法である「ヤコビ法」と「ガウス・ザイデル法」について簡単に説明します。
連立方程式 \(Ax = b\) の係数行列 \(A\) を対角行列 \(D\) とそれ以外の要素を持つ行列 \(A'\) に分解します。すると、方程式は次のように変形することができます。
ここで、最後の式を漸化式 \(x_{i+1} = D^{-1}(b - A'x_i)\) と考えて、反復処理により解を求めるのが「ヤコビ法 (Jacobi Method)」です。たとえば、三元連立方程式の場合、漸化式は下図のようになります。
a1 * x + a2 * y + a3 * z = d1
b1 * x + b2 * y + b3 * z = d2
c1 * x + c2 * y + c3 * z = d3
A' = [[ 0, a2, a3], D = [[a1, 0, 0], D-1 = [[1/a1, 0, 0],
[b1, 0, b3], [ 0, b2, 0], [ 0, 1/b2, 0],
{c1, c2, 0]] [ 0, 0, c3]] [ 0, 0, 1/c3]]
xi+1 = (d1 - a2 * yi - a3 * zi) / a1
yi+1 = (d2 - b1 * xi - b3 * zi) / b2
zi+1 = (d3 - c1 * xi - c2 * yi) / c3
ヤコビ法で解が収束した場合、それが元の連立方程式の解となります。なお、行列 A によっては解が収束しない (発散する) 場合があります。収束の十分条件ですが、『係数行列 \(A\) が対角優位な行列である場合に収束する』ことが知られています。行列 \(A\) の対角成分 \(a_{ii}\) の絶対値が、他の行の成分 \(a_{ij}\) の絶対値の合計よりも大きいことを「対角優位」といいます。
上図の三元連立方程式でいえば、|a1| > |b1| + |c1|, |b2| > |a2| + |c2|, |c3| > |a3| + |b3| を満たすとき、ヤコビ法の解は収束することになります。
収束の判定ですが、一般的には \(x_{i+1}\) と \(x_i\) の差分が許容誤差 \(\epsilon\) に収まったときに収束と判定します。具体的な方法を以下に示します。
1 は差分の合計値、2 は差分を解で割った値の合計値、3 は差分を解で割った値の最大値が ε 以下になったとき収束と判定します。今回は 2 の方法を使うことにします。
それではプログラムを作りましょう。Julia を使うとヤコビ法はとても簡単にプログラムすることができます。
リスト : ヤコビ法
function jacobi(a::Matrix{Float64}, b::Vector{Float64}; max_iter=512, eps=1e-6)
d = diag(a)
n = length(b)
x = zeros(n)
for i = 1 : max_iter
y = (b - (a * x - d .* x)) ./ d
if sum(abs.((y - x) ./ y)) < eps
return y, i
end
x = y
println("$i $x") # Debug
end
end
引数 a が係数行列、b が右辺値を格納した配列、キーワード引数 max_iter は繰り返しの最大回数、eps は許容誤差を表します。最初に、a の対角成分を diag() で取り出して変数 d にセットします。変数 x には解を格納する配列をセットします。x は 0 で初期化します。あとは for ループで漸化式 \(y = D^{-1}(b - A'x)\) を繰り返し計算します。実際には、A * x を計算して、そこから d .* x を引き算しています。Julia の場合、演算子 ./ は要素同士の割り算になるので、対角行列 d の逆行列を求める必要はありません。
収束条件のチェックも簡単です。sum(abs.((y - x) ./ y)) で合計値を求め、それが eps よりも小さくなったならば解は収束したので x を返します。for ループが終了した場合、解は収束していないので nothing を返します。
それでは実行してみましょう。次に示す方程式をヤコビ法で解きます。なお、変数 x の途中経過を表示するようにプログラムを修正しています。
9x + y + 2z = 9 x + 9y + z = 18 2x + y + 9z = -5
julia> a, b = jacobi([9. 1. 2.; 1. 9. 1.; 2. 1. 9.], [9., 18., -5.])
1 [1.0, 2.0, -0.555556]
2 [0.901235, 1.95062, -1.0]
3 [1.00549, 2.01097, -0.972565]
4 [0.992684, 1.99634, -1.00244]
5 [1.00095, 2.00108, -0.997968]
6 [0.999428, 1.99967, -1.00033]
7 [1.00011, 2.0001, -0.999836]
8 [0.999952, 1.99997, -1.00004]
9 [1.00001, 2.00001, -0.999986]
10 [0.999996, 2.0, -1.0]
11 [1.0, 2.0, -0.999999]
12 [1.0, 2.0, -1.0]
13 [1.0, 2.0, -1.0]
([1.0, 2.0, -1.0], 14)
julia> a
3-element Vector{Float64}:
0.9999999667086195
1.9999999761533003
-1.0000000318590117
厳密解は [1, -2, -1] ですが、ヤコビ法だと 14 回で収束しています。ガウス・ザイデル法を使うと、これよりも少ない回数で解を求めることができます。
ヤコビ法で漸化式を計算するとき、右辺式は xi, yi zi の値を使いますが、1 行目から順番に計算していくと、2 行目の yi+1 を計算するときには xi+1 の値がすでに求まっています。3 行目の zi+1 を計算するときは、xi+1 と yi+1 の値が求まっています。これらの値を使って漸化式を計算する方法を「ガウス・ザイデル法 (Gauss-Siedel Method)」といいます。
たとえば、三元連立方程式の場合、ガウス・ザイデル法の漸化式は下図のようになります。
a1 * x + a2 * y + a3 * z = d1
b1 * x + b2 * y + b3 * z = d2
c1 * x + c2 * y + c3 * z = d3
A' = [[ 0, a2, a3], D = [[a1, 0, 0], D-1 = [[1/a1, 0, 0],
[b1, 0, b3], [ 0, b2, 0], [ 0, 1/b2, 0],
{c1, c2, 0]] [ 0, 0, c3]] [ 0, 0, 1/c3]]
xi+1 = (d1 - a2 * yi - a3 * zi) / a1
yi+1 = (d2 - b1 * xi+1 - b3 * zi) / b2
zi+1 = (d3 - c1 * xi+1 - c2 * yi+1) / c3
ガウス・ザイデル法の解の収束条件や収束の判定方法はヤコビ法と同じです。
それではプログラムを作りましょう。次のリストを見てください。
リスト : ガウス・ザイデル法
function gauss_seidel(a::Matrix{Float64}, b::Vector{Float64}; max_iter=512, eps=1e-6)
n = length(b)
x = zeros(length(b))
for i = 1 : max_iter
err = 0.0
for j = 1 : n
s = b[j]
for k = 1 : n
if j != k
s -= a[j, k] * x[k]
end
end
s /= a[j, j]
err += abs((x[j] - s) / x[j])
x[j] = s
end
if err < eps
return x, i
end
println("$i $x") # Debug
end
end
ポイントは 2 番目の for ループで配列 x の値を逐次的に更新していくところです。これで漸化式を計算するときに計算済みの値を用いることができます。変数 err には誤差の合計値を格納します。2 番目の for ループが終了したとき、err が eps よりも小さい場合、解は収束したので x を返します。
それでは実際に試してみましょう。
9x + y + 2z = 9 x + 9y + z = 18 2x + y + 9z = -5
julia> a, b = gauss_seidel([9. 1. 2.; 1. 9. 1.; 2. 1. 9.], [9., 18., -5.])
1 [1.0, 1.88889, -0.987654]
2 [1.0096, 1.99756, -1.00186]
3 [1.00068, 2.00013, -1.00017]
4 [1.00002, 2.00002, -1.00001]
5 [1.0, 2.0, -1.0]
([1.0, 2.0, -1.0], 6)
julia> a
3-element Vector{Float64}:
0.9999999187463209
2.0000000119619115
-0.9999999832727282
ヤコビ法だと収束するまで 14 回かかっていたのが、ガウス・ザイデル法を使うと 6 回ですみました。
それでは簡単なテストとして、ガウスの消去法、ヤコビ法、ガウス・ザイデル法の実行時間を比較してみましょう。次のリストを見てください。
リスト : 簡単なテスト
# テストデータの作成
function makedata(n)
a = rand(n, n)
for i = 1 : n
a[i, i] = sum(a[:, i]) * 2
end
x = rand(n)
a, x, a * x
end
function test()
for n = [500, 1000, 1500, 2000]
a, x, b = makedata(n)
println("-----", n, "-----")
for (f, m) = [(gauss, "gauss"), (jacobi, "jacobi"), (gauss_seidel, "gauss_seidel")]
println(m)
@time ((x1, i) = f(a, b))
println(sum(abs.(x1 - x)), " $i")
end
end
end
乱数でテストデータを生成します。これを関数 makedata() で行っています。連立方程式 Ax = b の係数行列 A と解 x を乱数で生成して、右辺値の b を計算します。係数行列は対角優位になるよう修正しています。makedata() は A, x, b を返すので、A と b から x を求め、その誤差と実行時間を表示します。
実行結果は次のようになりました。
julia> test() -----500----- gauss 0.370358 seconds (280.63 k allocations: 23.013 MiB, 81.47% compilation time) 2.772894761027267e-13 0 jacobi 0.222565 seconds (223.23 k allocations: 15.954 MiB, 99.03% compilation time) 1.3331930535996454e-8 34 gauss_seidel 0.066094 seconds (10.85 k allocations: 744.867 KiB, 91.86% compilation time) 9.877097024383469e-10 12 -----1000----- gauss 1.684524 seconds (2.01 k allocations: 15.554 MiB, 0.37% gc time) 8.471478726845838e-13 0 jacobi 0.020860 seconds (245 allocations: 1.861 MiB) 2.756374565789435e-8 34 gauss_seidel 0.051413 seconds (6 allocations: 8.094 KiB) 2.0114432242208935e-9 12 -----1500----- gauss 14.205831 seconds (3.01 k allocations: 34.788 MiB, 0.05% gc time) 1.5606604008254477e-12 0 jacobi 0.052806 seconds (252 allocations: 2.865 MiB) 2.100354414957123e-8 35 gauss_seidel 0.115528 seconds (6 allocations: 12.031 KiB) 3.0632216258441765e-9 12 -----2000----- gauss 41.330912 seconds (4.01 k allocations: 61.651 MiB, 0.24% gc time) 2.3809471375816488e-12 0 jacobi 0.087211 seconds (259 allocations: 3.907 MiB) 1.4529886723505707e-8 36 gauss_seidel 0.370988 seconds (6 allocations: 15.906 KiB) 4.167792542320416e-9 12 実行環境 : Julia ver 1.10.5, Ubuntu 22.04 (WSL2), Intel Core i5-6200U 2.30GHz
変数の個数が増えるにしたがい、ガウスの消去法の実行時間は大幅に増加します。反復法でも遅くはなりますが、ガウスの消去法よりはずっと高速です。ヤコビ法とガウス・ザイデル法を比較した場合、ヤコビ法のほうが速くなる場合もありました。反復回数はガウス・ザイデル法のほうが少ないので、逐次的に計算する処理で時間がかかっていると思われます。Julia の場合、for ループはそれほど遅くはないのですが、行列で計算したほうが速くなるようです。興味のある方はいろいろ試してみてください。
リスト : 簡単なテスト
using LinearAlgebra
# ピボット選択
function select_pivot(xs::Matrix{Float64}, i::Int)
_, k = findmax(abs.(xs[i:end, i]))
k += i - 1
if k != i
temp = xs[i, :]
xs[i, :] = xs[k, :]
xs[k, :] = temp
end
end
# ガウスの消去法
function gauss(xs::Matrix{Float64}, ys::Vector{Float64})
zs = hcat(xs, ys)
n = size(xs, 1)
# 前進消去
for i = 1 : n - 1
# ピボット選択
select_pivot(zs, i)
for j = i + 1 : n
temp = zs[j, i] / zs[i, i]
for k = i : n + 1
zs[j, k] -= temp * zs[i, k]
end
end
end
# 後退代入
for i = n : -1 : 1
for j = i + 1 : n
zs[i, n + 1] -= zs[i, j] * zs[j, n + 1]
end
zs[i, n + 1] /= zs[i, i]
end
zs[:, n + 1], 0
end
# ヤコビ法
function jacobi(a::Matrix{Float64}, b::Vector{Float64}; max_iter=512, eps=1e-6)
d = diag(a)
n = length(b)
x = zeros(n)
for i = 1 : max_iter
y = (b - (a * x - d .* x)) ./ d
if sum(abs.((y - x) ./ y)) < eps
return y, i
end
x = y
end
end
# ガウス・ザイデル法
function gauss_seidel(a::Matrix{Float64}, b::Vector{Float64}; max_iter=512, eps=1e-6)
n = length(b)
x = zeros(length(b))
for i = 1 : max_iter
err = 0.0
for j = 1 : n
s = b[j]
for k = 1 : n
if j != k
s -= a[j, k] * x[k]
end
end
s /= a[j, j]
err += abs((x[j] - s) / x[j])
x[j] = s
end
if err < eps
return x, i
end
end
end
# テストデータの作成
function makedata(n)
a = rand(n, n)
for i = 1 : n
a[i, i] = sum(a[:, i]) * 2
end
x = rand(n)
a, x, a * x
end
function test()
for n = [500, 1000, 1500, 2000]
a, x, b = makedata(n)
println("-----", n, "-----")
for (f, m) = [(gauss, "gauss"), (jacobi, "jacobi"), (gauss_seidel, "gauss_seidel")]
println(m)
@time ((x1, i) = f(a, b))
println(sum(abs.(x1 - x)), " $i")
end
end
end