Julia Language Programming

線形代数編

●行列とベクトルの基礎知識

• 行列 (matrix) は要素 (element) を長方形または正方形に並べた「表 (table)」みたいなもの
• 要素は数、式、記号など
• 要素のことを成分ともいう
• Julia では 2 次元配列 Array{T, 2} が行列を表す (別名 Matrix{T})
• 行列の横方向の並びを「行 (row)」といい、縦方向の並びを「列 (column)」という
• すべての成分が 0 の行列を「零行列」という
• 行の数と列の数が等しい行列 (n 行 n 列の行列) を「正方行列」という
• 正方行列 \(A\) の左上から右下の対角成分が 1 で、それ以外の成分が 0 の行列を「単位行列」という
• 単位行列は \(E\) や \(I\) と表記されることが多い

  [ 1 2 3;     [ 1 2;     [ 1 2 3:     [ 1 0 0;
    4 5 6 ]      3 4;       4 5 6;       0 1 0;
                 5 6 ]      7 8 9 ]      0 0 1 ]
  
  2 行 3 列   3 行 2 列   3 行 3 列    単位行列
  

• 行列 \(A\) と \(B\) において、行の数 m と列の数 n が等しいとき、\(A\) と \(B\) を同じ型であるという
• 同じ型の行列 \(A\) と \(B\) において、対応する成分がすべて等しいとき、\(A\) と \(B\) は等しい (\(A = B\))
• 1 行だけの行列を「行ベクトル」、1 列だけの行列を「列ベクトル」という

• 大きさ (長さ) と方向を持つ量のことを「ベクトル (vector)」という
• 大きさだけで方向を持たない量のことを「スカラー (scalar)」という
• Julia では 1 次元配列 Array{T, 1} がベクトルを表す (別名 Vector{T})
• ベクトルだけの演算であれば、行ベクトルと列ベクトルの区別は不要
• 行列とベクトルを演算する場合、行ベクトルと列ベクトルを考慮する必要がある
• Julia の場合、1 次元配列は列ベクトルとして扱われるようだ
• 大きさが 0 のベクトルを「ゼロベクトル (零ベクトル)」という
• 大きさ 1 のベクトルを「単位ベクトル」という

• ベクトルの計算
  • 要素同士の演算は演算子の前にドット ( . ) を付けることで行う
  • 和と差は演算子 +, - でもできる
  • ベクトルのスカラー倍は演算子 *, / を使う
  • 以下の関数を使用するときは using LinearAlgebra が必要
  • ベクトル v, w の内積は関数 dot(v, w) を使う
    • 演算子 * は行列の乗算で使用する (行ベクトル * 列ベクトル は演算できる)
  • ベクトル v の大きさ (ノルム) \(|v|\) は関数 norm(v) で求める

julia> using LinearAlgebra

julia> a = rand(3)
3-element Vector{Float64}:
 0.09356090585693044
 0.6760555446095466
 0.9254304928463559

julia> b = rand(3)
3-element Vector{Float64}:
 0.8520683072471729
 0.9422268415336021
 0.1967745050319618

julia> a + b
3-element Vector{Float64}:
 0.9456292131041033
 1.6182823861431488
 1.1222049978783177

julia> a - b
3-element Vector{Float64}:
 -0.7585074013902424
 -0.2661712969240555
  0.7286559878143941

julia> a - a
3-element Vector{Float64}:
 0.0
 0.0
 0.0

julia> a .* b
3-element Vector{Float64}:
 0.07972028267802682
 0.6369976804987324
 0.18210112717132615

julia> a ./ b
3-element Vector{Float64}:
 0.10980446645082148
 0.7175082632003716
 4.702999978051219

julia> dot(a, b)
0.8988190903480853

julia> norm(a)
1.14988101105807

julia> ea = a / norm(a)
3-element Vector{Float64}:
 0.08136572824247251
 0.5879352194775963
 0.8048054398209563

julia> norm(ea)
1.0

• 行列の計算 (1)
  • 要素同士の演算は演算子の前にドット ( . ) を付けることで行う
  • 和と差は演算子 +, - でもできる
  • 行列のスカラー倍は演算子 *, / を使う

julia> a = rand(2, 3)
2×3 Matrix{Float64}:
 0.186229  0.0531527  0.526371
 0.462431  0.0348494  0.0496523

julia> b = rand(2, 3)
2×3 Matrix{Float64}:
 0.169693  0.124246  0.742511
 0.797795  0.963666  0.685439

julia> a + b
2×3 Matrix{Float64}:
 0.355922  0.177399  1.26888
 1.26023   0.998515  0.735092

julia> a - b
2×3 Matrix{Float64}:
  0.0165354  -0.0710936  -0.21614
 -0.335364   -0.928816   -0.635787

julia> a - a
2×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> a .* b
2×3 Matrix{Float64}:
 0.0316017  0.00660402  0.390836
 0.368925   0.0335831   0.0340337

julia> a ./ b
2×3 Matrix{Float64}:
 1.09744   0.427801   0.708906
 0.579636  0.0361633  0.0724387

julia> a * 10
2×3 Matrix{Float64}:
 1.86229  0.531527  5.26371
 4.62431  0.348494  0.496523

julia> a / 10
2×3 Matrix{Float64}:
 0.0186229  0.00531527  0.0526371
 0.0462431  0.00348494  0.00496523

• 行列の計算 (2)
  • 行列の積は *, 累乗は ^, 逆行列 \(A^{-1}\) は関数 inv(A) を使う
    • 累乗 ^ と逆行列 inv() は正方行列であること
  • 正方行列 \(A\) に逆行列 \(A^{-1}\) が存在するとき、\(A\) は「正則行列 (regular matrix)」であるという
  • \(A\) が正則行列であれば、連立方程式 \(Ax = b\) の解 \(x\) は一意に定まる
  • Julia の場合、\(x\) を inv(A) * b で求めるよりも、演算子 \ を使ったほうがよいようだ
  • x = A \ b (A は正方行列、x と b はベクトル)

julia> c = rand(3, 2)
3×2 Matrix{Float64}:
 0.0876367  0.424824
 0.865688   0.593409
 0.91316    0.827241

julia> d = a * c
2×2 Matrix{Float64}:
 0.542995  0.546091
 0.116035  0.258206

julia> e = [1. 0.; 0. 1.]
2×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0

julia> d * e
2×2 Matrix{Float64}:
 0.542995  0.546091
 0.116035  0.258206

julia> d * e == d
true

julia> e * d
2×2 Matrix{Float64}:
 0.542995  0.546091
 0.116035  0.258206

julia> e * d == d
true

julia> e * e
2×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0

julia> e * e == e
true

julia> f = [1 1; 1 0]
2×2 Matrix{Int64}:
 1  1
 1  0

julia> f ^ 2
2×2 Matrix{Int64}:
 2  1
 1  1

julia> f ^ 3
2×2 Matrix{Int64}:
 3  2
 2  1

julia> f ^ 4
2×2 Matrix{Int64}:
 5  3
 3  2

julia> f ^ 10
2×2 Matrix{Int64}:
 89  55
 55  34

julia> f ^ 40
2×2 Matrix{Int64}:
 165580141  102334155
 102334155   63245986

julia> a1 = [2 3; 4 7]
2×2 Matrix{Int64}:
 2  3
 4  7

julia> ra1 = inv(a1)
2×2 Matrix{Float64}:
  3.5  -1.5
 -2.0   1.0

julia> a1 * ra1
2×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0

julia> ra1 * [4, 6]
2-element Vector{Float64}:
  5.0
 -2.0

julia> a1 \ [4, 6]
2-element Vector{Float64}:
  5.0
 -2.0

• 行列の計算 (3)
  • 以下の関数を使用するときは using LinearAlgebra が必要
  • det(A), 行列 A の行列式 (determinant)
  • 行列 A が正則行列であれば、det(A) は 0 にならない
  • tr(A), 行列 A のトレース (trace, 跡)
    • トレースは行列の対角線上にある要素の和のこと (\(a_{11} + a_{22} + \cdots + a_{nn}\))
  • rank(A), 行列 A のランク (rank)
  • 行列 A が n 次の正則行列であれば、rank(A) の値は n になる

julia> det(a1)
2.0

julia> tr(a1)
9

julia> rank(a1)
2

julia> b1 = [1 2 3; 4 5 6; 7 8 9]
3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9

julia> det(b1)
0.0

julia> tr(b1)
15

julia> rank(b1)
2

julia> inv(b1)
ERROR: SingularException(3)

●特殊な行列とその性質

• 転置行列 (transposed matrix)
  • 行列 A の行と列を入れ替えたもので、\(A^{\mathrm{T}}\) と表記する
  • julia の場合、A' または関数 transpose(A) で転置行列を求めることができる

julia> a = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> a'
2×2 adjoint(::Matrix{Int64}) with eltype Int64:
 1  3
 2  4

julia> transpose(a)
2×2 transpose(::Matrix{Int64}) with eltype Int64:
 1  3
 2  4

julia> b = reshape(collect(0:8), 3, 3)
3×3 Matrix{Int64}:
 0  3  6
 1  4  7
 2  5  8

julia> b'
3×3 adjoint(::Matrix{Int64}) with eltype Int64:
 0  1  2
 3  4  5
 6  7  8

julia> transpose(b)
3×3 transpose(::Matrix{Int64}) with eltype Int64:
 0  1  2
 3  4  5
 6  7  8

julia> c = reshape(collect(0:5), 3, 2)
3×2 Matrix{Int64}:
 0  3
 1  4
 2  5

julia> c'
2×3 adjoint(::Matrix{Int64}) with eltype Int64:
 0  1  2
 3  4  5

julia> transpose(c)
2×3 transpose(::Matrix{Int64}) with eltype Int64:
 0  1  2
 3  4  5

• 転置行列の基本的な性質 (A, B は正方行列)
  1. \(tr(A^{\mathrm{T}}) = tr(A)\)
  2. \(|A^{\mathrm{T}}| = |A|\), (\(|A|\) は \(A\) の行列式)
  3. \((AB)^{\mathrm{T}} = B^{\mathrm{T}}A^{\mathrm{T}}\)
  4. \((A^{\mathrm{T}})^{-1} = (A^{-1})^{\mathrm{T}}\)

julia> tr(a)
5

julia> tr(a')
5

julia> det(a)
-2.0

julia> det(a')
-2.0

julia> d = [5 6; 7 8]
2×2 Matrix{Int64}:
 5  6
 7  8

julia> (a * d)'
2×2 adjoint(::Matrix{Int64}) with eltype Int64:
 19  43
 22  50

julia> d' * a'
2×2 Matrix{Int64}:
 19  43
 22  50

julia> inv(a')
2×2 adjoint(::Matrix{Float64}) with eltype Float64:
 -2.0   1.5
  1.0  -0.5

julia> inv(a)'
2×2 adjoint(::Matrix{Float64}) with eltype Float64:
 -2.0   1.5
  1.0  -0.5

• 対称行列と対角行列
  • 正方行列 \(A\) とその転置行列 \(A^{\mathrm{T}}\) が等しい行列のことを「対称行列 (symmetric matrix)」という
  • その中で対角成分 \(a_{ii}\) 以外の要素が 0 の行列のことを「対角行列 (diagonal matrix)」という
  • 単位行列は対角行列の一種

julia> a = [1 3; 3 2]
2×2 Matrix{Int64}:
 1  3
 3  2

julia> a' == a
true

julia> b = [1 4 5; 4 2 6; 5 6 3]
3×3 Matrix{Int64}:
 1  4  5
 4  2  6
 5  6  3

julia> b' == b
true

julia> c = [1 0 0; 0 2 0; 0 0 3]
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

julia> c' == c
true

• 操作関数
  • Symmetric(A, uplo=:U), 行列 \(A\) から対称行列 (view) を生成する
    • :U は upper, 対角線の上半分と対称になる
    • :L は lower, 対角線の下半分と対称になる
    • 対角成分以外の要素は書き換え不可
    • データ型は対称行列を表すものになる
  • Diagonal(A or v), 対角行列の生成
    • 引数に行列 \(A\) を渡すと \(A\) の対角成分と同じ対角行列を生成する
    • 引数にベクトル v を渡すと v の要素が対角成分になる
    • データ型は対角行列を表すものになる
  • diag(A, k = 0), 行列 \(A\) の対角線上の要素を取得する
    • k = 0 は主対角線を取り出す
    • K > 0 は上側の、k < 0 は下側の副対角線を取り出す
  • diagm(k => v), ベクトル v から対角行列を生成する
    • k は diag() と同じ意味
    • k => v は複数指定できる
    • データ型は普通の配列

julia> a = [1 2 3; 4 5 6; 7 8 10]
3×3 Matrix{Int64}:
 1  2   3
 4  5   6
 7  8  10

julia> Symmetric(a)
3×3 Symmetric{Int64,Matrix{Int64}}:
 1  2   3
 2  5   6
 3  6  10

julia> Symmetric(a, :L)
3×3 Symmetric{Int64,Matrix{Int64}}:
 1  4   7
 4  5   8
 7  8  10

julia> Diagonal(a)
3×3 Diagonal{Int64, Vector{Int64}}:
 1  ⋅   ⋅
 ⋅  5   ⋅
 ⋅  ⋅  10

julia> Diagonal([1,2,3,4])
4×4 Diagonal{Int64, Vector{Int64}}:
 1  ⋅  ⋅  ⋅
 ⋅  2  ⋅  ⋅
 ⋅  ⋅  3  ⋅
 ⋅  ⋅  ⋅  4

julia> diag(a)
3-element Vector{Int64}:
  1
  5
 10

julia> diag(a, 1)
2-element Vector{Int64}:
 2
 6

julia> diag(a, 2)
1-element Vector{Int64}:
 3

julia> diag(a, -1)
2-element Vector{Int64}:
 4
 8

julia> diag(a, -2)
1-element Vector{Int64}:
 7

julia> diagm(0 => [1, 2, 3])
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

julia> diagm(1 => [1, 2, 3])
4×4 Matrix{Int64}:
 0  1  0  0
 0  0  2  0
 0  0  0  3
 0  0  0  0

julia> diagm(-1 => [1, 2, 3])
4×4 Matrix{Int64}:
 0  0  0  0
 1  0  0  0
 0  2  0  0
 0  0  3  0

julia> diagm(0 => [1, 1, 1, 1], -1 => [1, 2, 3])
4×4 Matrix{Int64}:
 1  0  0  0
 1  1  0  0
 0  2  1  0
 0  0  3  1

• 対角行列の基本的な性質 (\(A\) は対角行列)
  1. \(|A|\) は対角成分 \(a_{ii}\) を乗算した値になる
  2. \(A^n\) は対角成分 \(a_{ii}\) を n 乗したものになる
  3. \(A^{-1}\) は対角成分 \(a_{ii}\) を逆数にしたものになる
  4. 対角行列の固有値は対角成分 \(a_{ii}\) になる (固有値は後で説明する)

julia> a = Diagonal([1,2,3])
3×3 Diagonal{Int64,Vector{Int64}}:
 1  ⋅  ⋅
 ⋅  2  ⋅
 ⋅  ⋅  3

julia> det(a)
6

julia> diag(a)
3-element Vector{Int64}:
 1
 2
 3

julia> prod(diag(a))
6

julia> a * a
3×3 Diagonal{Int64,Vector{Int64}}:
 1  ⋅  ⋅
 ⋅  4  ⋅
 ⋅  ⋅  9

julia> a * a * a
3×3 Diagonal{Int64,Vector{Int64}}:
 1  ⋅   ⋅
 ⋅  8   ⋅
 ⋅  ⋅  27

julia> a ^ 2
3×3 Diagonal{Int64,Vector{Int64}}:
 1  ⋅  ⋅
 ⋅  4  ⋅
 ⋅  ⋅  9

julia> a ^ 3
3×3 Diagonal{Int64,Vector{Int64}}:
 1  ⋅   ⋅
 ⋅  8   ⋅
 ⋅  ⋅  27

julia> inv(a)
3×3 Diagonal{Float64,Vector{Float64}}:
 1.0   ⋅    ⋅
  ⋅   0.5   ⋅
  ⋅    ⋅   0.333333

julia> Diagonal([1/1, 1/2, 1/3])
3×3 Diagonal{Float64,Vector{Float64}}:
 1.0   ⋅    ⋅
  ⋅   0.5   ⋅
  ⋅    ⋅   0.333333

julia> eigvals(a)  # 行列 a の固有値を求める
3-element Vector{Int64}:
 1
 2
 3

• 三角行列
  • 上三角行列 (Upper triangular matrix) は対角成分よりも下側の成分が 0 である行列のこと
  • 下三角行列 (lower triangular matrix) は対角成分よりも上側の成分が 0 である行列のこと
  • UpperTriangular(A), 行列 \(A\) から上三角行列 (view) を生成する
  • LowerTriangular(A), 行列 \(A\) から下三角行列 (view) を生成する
  • 生成された view は上 (下) 三角行列を表すデータ型になる

• 三角行列の主な性質 (\(A\) は三角行列, \(L\) は下三角行列, \(U\) は上三角行列)
  1. \(L^{\mathrm{T}}\) は \(U\) に, \(U^{\mathrm{T}}\) は \(L\) になる
  2. \(U * U\) は \(U\) に、\(L * L\) は \(L\) になる
  3. \(|A|\) は対角成分 \(a_{ii}\) を乗算した値になる
  4. \(U^{-1}\) は \(U\) に、\(L^{-1}\) は \(L\) になる
  5. \(A\) の固有値は対角成分 \(a_{ii}\) になる (固有値は後で説明する)

julia> a = [1 2 3; 4 5 6; 7 8 10]
3×3 Matrix{Int64}:
 1  2   3
 4  5   6
 7  8  10

julia> la = LowerTriangular(a)
3×3 LowerTriangular{Int64,Matrix{Int64}}:
 1  ⋅   ⋅
 4  5   ⋅
 7  8  10

julia> ua = UpperTriangular(a)
3×3 UpperTriangular{Int64,Matrix{Int64}}:
 1  2   3
 ⋅  5   6
 ⋅  ⋅  10

julia> la'
3×3 UpperTriangular{Int64, Adjoint{Int64, Matrix{Int64}}}:
 1  4   7
 ⋅  5   8
 ⋅  ⋅  10

julia> ua'
3×3 LowerTriangular{Int64, Adjoint{Int64, Matrix{Int64}}}:
 1  ⋅   ⋅
 2  5   ⋅
 3  6  10

julia> la * la
3×3 LowerTriangular{Int64, Matrix{Int64}}:
   1    ⋅    ⋅
  24   25    ⋅
 109  120  100

julia> ua * ua
3×3 UpperTriangular{Int64,Matrix{Int64}}:
 1  12   45
 ⋅  25   90
 ⋅   ⋅  100

julia> det(la)
50

julia> det(ua)
50

julia> inv(la)
3×3 LowerTriangular{Float64,Matrix{Float64}}:
  1.0     ⋅     ⋅
 -0.8    0.2    ⋅
 -0.06  -0.16  0.1

julia> inv(ua)
3×3 UpperTriangular{Float64,Matrix{Float64}}:
 1.0  -0.4  -0.06
  ⋅    0.2  -0.12
  ⋅     ⋅    0.1

julia> eigvals(la)
3-element Vector{Int64}:
  1
  5
 10

julia> eigvals(ua)
3-element Vector{Int64}:
  1
  5
 10

• 直交行列 (orthogonal matrix)
  • \(A^{\mathrm{T}}\) と \(A^{-1}\) が等しい行列のこと
  • 直交行列の簡単な例として、回転行列や置換行列などがある

• 直交行列の基本的な性質
  1. 直交行列の行列式は \(\pm 1\)
  2. 直交行列の逆行列も直交行列
  3. 対称行列は直交行列で対角化できる (あとで説明する)

回転行列 R は直交行列

R = [[cos(x), -sin(x)],
     [sin(x),  cos(x)]]

転置行列
RT = [[cos(x), sin(x)]
      [-sin(x), cos(x)]]

行列式 |R| = cos(x) * cos(x) + sin(x) * sin(x) = 1

逆行列
R-1 = [[cos(x), sin(x)],
       [-sin(x), cos(x)]]          

逆行列の転置 (R に戻るので直交行列)
(R-1)T = [[cos(x), -sin(x)],
          [sin(x), cos(x)]]

julia> a = [0 1; 1 0]
2×2 Matrix{Int64}:
 0  1
 1  0

julia> a'
2×2 adjoint(::Matrix{Int64}) with eltype Int64:
 0  1
 1  0

julia> det(a)
-1.0

julia> inv(a)
2×2 Matrix{Float64}:
 -0.0  1.0
  1.0  0.0

julia> b = [0 1 0; 1 0 0; 0 0 1]
3×3 Matrix{Int64}:
 0  1  0
 1  0  0
 0  0  1

julia> b'
3×3 adjoint(::Matrix{Int64}) with eltype Int64:
 0  1  0
 1  0  0
 0  0  1

julia> det(b)
-1.0

julia> inv(b)
3×3 Matrix{Float64}:
 -0.0  1.0  -0.0
  1.0  0.0  -0.0
  0.0  0.0   1.0

julia> c = [1 2 3; 4 5 6; 7 8 9]
3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9

julia> c * b        # 右から b を掛けると列の交換
3×3 Matrix{Int64}:
 2  1  3
 5  4  6
 8  7  9

julia> b * c        # 左から b を掛けると行の交換
3×3 Matrix{Int64}:
 4  5  6
 1  2  3
 7  8  9

• 各行各列にちょうど一つだけ 1 の要素を持ち、それ以外は全て 0 となるような正方行列を「置換行列」という
• 置換行列を掛け算すると、列または行を入れ替える働きをする

●連立一次方程式の解法

今回は連立一次方程式を解く基本的なアルゴリズムを取り上げます。Julia には連立一次方程式 \(Ax = b\) の \(x\) を求める演算子 \ (x = A \ b) が用意されているので、実用的にはそちらを使えばいいのですが、Julia とアルゴリズムのお勉強ということで、あえてプログラムを作ってみましょう。

●ガウスの消去法

ガウスの消去法 (Gaussian elimination) は人が連立方程式を解くときの方法とほとんど同じです。次の図を見てください。

a1 * x + a2 * y + a3 * z = d1       a1 * x + a2  * y + a3  * z = d1        a1 * x + a2  * y + a3   * z = d1
b1 * x + b2 * y + b3 * z = d2  ==>           b2' * y + b3' * z = d2'  ==>           b2' * y + b3'  * z = d2'
c1 * x + c2 * y + c3 * z = d3                c2' * y + c3' * z = d3'                          c3'' * z = d3''

                                          図 : 前進消去

1 番目の式を b1 / a1 倍して 2 番目の式から引き算すると、x の項を消去することができます。同様の方法で 3 番目の式の x の項を消去することができます。次に、2 番目の式を c2' / b2' 倍して 3 番目の式から引き算すると y の項を消去することができます。この処理を「前進消去」といいます。

前進消去を行うと、最後の式には変数が z しかありません。z の値は d3'' / c3'' に決定することができます。次はひとつ前の式に戻り、 変数 z に d3'' / c3'' を代入すると変数 y の値を求めることができます。これを繰り返して、先頭の式に戻ると変数 x の値を求めることができます。この処理を「後退代入」といいます。

連立一次方程式は行列を使って Ax = b と表すことができます。A を係数行列といい、A と b を連結した行列を拡大係数行列といいます。

[ a1, a2, a3;    [ a1, a2, a3, d1;
  b1, b2, b3;      b1, b2, b3, d2;
  c1, c2, c3 ]     c1, c2, c3, d3 ]

   係数行列         拡大係数行列

ガウスの消去法は、拡大係数行列を使うと簡単にプログラムを作ることができます。Julia でプログラムすると、次のようになります。

リスト : ガウスの消去法

function gauss(xs::Matrix{Float64}, ys::Vector{Float64})
    # 拡大係数行列
    zs = hcat(xs, ys)
    n = size(xs, 1)
    # 前進消去
    for i = 1 : n - 1
        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
        println(zs)    # debug
    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]
end

関数 gauss() の引数 xs が係数行列、ys が右辺の定数を格納したベクトルです。まず最初に、拡大係数行列を生成して変数 zs にセットします。あとはアルゴリズムをそのままプログラムしただけです。Julia の for ループはかなり速いので、配列の indexing や dot() などの機能を使わなくても、このままで十分だと思います。

それでは実際に試してみましょう。次に示す連立方程式を解いてみました。

リスト : 簡単なテスト

#  x +  y = 100
# 2x + 4y = 272
println(gauss([1. 1.; 2. 4.], [100., 272.]))

#  x +  y +  z = 10
# 2x + 4y + 6z = 38
# 2x      + 4z = 14
println(gauss([1. 1. 1.; 2. 4. 6.; 2. 0. 4.], [10., 38., 14.]))

#   a +  b +   c + d = -5
#  -a +  b -   c + d = -7
#  8a + 4b +  2c + d = -31
# -8a + 4b + -2c + d = -35
println(gauss([1. 1. 1. 1.; -1. 1. -1. 1.; 8. 4. 2. 1.; -8. 4. -2. 1.], [-5., -7., -31., -35.]))

#   a -  b +  c - d + e = 1
# 12a - 6b + 2c         = 0 
#   a +  b +  c + d + e = 8
# 12a + 6b + 2c         = 0
#  4a + 3b + 2c + d     = 1
println(gauss([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.], 
[1.,0.,8.,0.,1.]))

実行結果は次のようになります。なお、拡大行列の表示は手作業で整形しています。

[1.0 1.0 100.0;
 0.0 2.0 72.0]
[64.0, 36.0]

[1.0  1.0  1.0  10.0;
 0.0  2.0  4.0  18.0;
 0.0 -2.0  2.0  -6.0]
[1.0  1.0  1.0  10.0;
 0.0  2.0  4.0  18.0;
 0.0  0.0  6.0  12.0]
[3.0, 5.0, 2.0]

[1.0  1.0  1.0  1.0  -5.0;
 0.0  2.0  0.0  2.0 -12.0;
 0.0 -4.0 -6.0 -7.0   9.0;
 0.0 12.0  6.0  9.0 -75.0]
[1.0  1.0  1.0  1.0  -5.0;
 0.0  2.0  0.0  2.0 -12.0;
 0.0  0.0 -6.0 -3.0 -15.0;
 0.0  0.0  6.0 -3.0  -3.0]
[1.0  1.0  1.0  1.0  -5.0;
 0.0  2.0  0.0  2.0 -12.0;
 0.0  0.0 -6.0 -3.0 -15.0;
 0.0  0.0  0.0 -6.0 -18.0]
[0.0, -9.0, 1.0, 3.0]

[1.0 -1.0   1.0      -1.0           1.0           1.0;
 0.0  6.0 -10.0      12.0         -12.0         -12.0;
 0.0  2.0   0.0       2.0           0.0           7.0;
 0.0 18.0 -10.0      12.0         -12.0         -12.0;
 0.0  7.0  -2.0       5.0          -4.0          -3.0]
[1.0 -1.0   1.0      -1.0           1.0           1.0;
 0.0  6.0 -10.0      12.0         -12.0         -12.0;
 0.0  0.0   3.33333  -2.0           4.0          11.0;
 0.0  0.0  20.0     -24.0          24.0          24.0;
 0.0  0.0   9.66667  -9.0          10.0          11.0]
[1.0 -1.0   1.0      -1.0           1.0           1.0;
 0.0  6.0 -10.0      12.0         -12.0         -12.0;
 0.0  0.0   3.33333  -2.0           4.0          11.0;
 0.0  0.0   0.0     -12.0          -3.55271e-15 -42.0;
 0.0  0.0   0.0      -3.2          -1.6         -20.9]
[1.0 -1.0   1.0      -1.0           1.0           1.0;
 0.0  6.0 -10.0      12.0         -12.0         -12.0;
 0.0  0.0   3.33333  -2.0           4.0          11.0;
 0.0  0.0   0.0     -12.0          -3.55271e-15 -42.0;
 0.0  0.0   0.0      -4.44089e-16  -1.6          -9.7]
[0.3125, 0.0, -1.875, 3.5, 6.0625]

●ガウス・ジョルダン法

ガウスの消去法は、係数行列の部分を「上三角行列」に変換することで、連立一次方程式を解きました。ここで、係数行列を単位行列に変換できれば、後退代入することなく解を求めることができます。これを「ガウス・ジョルダン法 (Gauss - Jordan elimination)」といいます。簡単な動作例を下図に示します。

次の連立方程式をガウス・ジョルダン法で解く

  x +  y +  z = 10
 2x + 4y + 6z = 38
 2x      + 4z = 14

拡大係数行列を zs とする

(1) zs[1, 1] を 1 にする
[  1.   1.   1.  10.;    # zs[1, 1:end] /= 1
   2.   4.   6.  38.;
   2.   0.   4.  14.]

(2) zs[2, 1], zs[3, 1] を 0 にする
[  1.   1.   1.  10.;     
   0.   2.   4.  18.;     # zs[2, 1:end] -= zs[1, 1:end] * 2
   0.  -2.   2.  -6.]     # zs[3, 1:end] -= zs[1, 1:end] * 2

(3) zs[2, 2] を 1 にする
[  1.   1.   1.  10.;
   0.   1.   2.   9.;     # zs[2, 2:end] /= 2
   0.  -2.   2.  -6.]

(4) zs[1, 2], zs[3, 2] を 0 にする
[  1.   0.  -1.   1.;     # zs[1, 2:end] -= zs[2, 2:end] * 1
   0.   1.   2.   9.;
   0.   0.   6.  12.]     # zs[3, 2:end] -= zs[2, 2:end] * -2

(5) zs[3, 3] を 1 にする
[  1.   0.  -1.   1.;
   0.   1.   2.   9.;
   0.   0.   1.   2.]     # zs[3, 3:end] /= 6

(6) zs[1, 3], zs[2, 3] を 0 にする
[  1.   0.   0.   3.;     # zs[1, 3:end] -= zs[3, 3:end] * -1
   0.   1.   0.   5.;     # zs[2, 3:end] -= zs[3, 3:end] * 2
   0.   0.   1.   2.]

このように、ガウス・ジョルダン法の基本的な考え方は実にシンプルで、プログラムも簡単になるのですが、速度はガウスの消去法よりも少し遅くなるといわれています。以下にプログラムと実行例を示します。

リスト : ガウス・ジョルダン法

function gauss_jordan(xs::Matrix{Float64}, ys::Vector{Float64})
    # 拡大係数行列
    zs = hcat(xs, ys)
    n = size(xs, 1)
    for i = 1 : n
        # zs[i, i] を 1 にする
        temp = zs[i, i]
        zs[i, i:end] /= temp
        # zs[i] 以外の行の i 番目の係数を 0 にする
        for j = 1 : n
            if i != j
                temp = zs[j, i]
                for k = 1 : n + 1
                    zs[j, k] -= zs[i, k] * temp
                end
            end
        end
        println(zs)   # debug
    end
    zs[:, n + 1]
end

[1.0 1.0 100.0;
 0.0 2.0 72.0]
[1.0 0.0 64.0;
 0.0 1.0 36.0]
[64.0, 36.0]

[1.0  1.0  1.0 10.0;
 0.0  2.0  4.0 18.0;
 0.0 -2.0  2.0 -6.0]
[1.0  0.0 -1.0  1.0;
 0.0  1.0  2.0  9.0;
 0.0  0.0  6.0 12.0]
[1.0  0.0  0.0  3.0;
 0.0  1.0  0.0  5.0;
 0.0  0.0  1.0  2.0]
[3.0, 5.0, 2.0]

[1.0  1.0  1.0  1.0  -5.0;
 0.0  2.0  0.0  2.0 -12.0;
 0.0 -4.0 -6.0 -7.0   9.0;
 0.0 12.0  6.0  9.0 -75.0]
[1.0  0.0  1.0  0.0   1.0;
 0.0  1.0  0.0  1.0  -6.0;
 0.0  0.0 -6.0 -3.0 -15.0;
 0.0  0.0  6.0 -3.0  -3.0]
[1.0  0.0  0.0 -0.5  -1.5;
 0.0  1.0  0.0  1.0  -6.0;
 0.0  0.0  1.0  0.5   2.5;
 0.0  0.0  0.0 -6.0 -18.0]
[1.0  0.0  0.0  0.0   0.0;
 0.0  1.0  0.0  0.0  -9.0;
 0.0  0.0  1.0  0.0   1.0;
 0.0  0.0  0.0  1.0   3.0]
[0.0, -9.0, 1.0, 3.0]

[1.0 -1.0    1.0 -1.0   1.0   1.0;
 0.0  6.0  -10.0 12.0 -12.0 -12.0;
 0.0  2.0    0.0  2.0   0.0   7.0;
 0.0  18.0 -10.0 12.0 -12.0 -12.0;
 0.0  7.0   -2.0  5.0  -4.0  -3.0]
[1.0  0.0 -0.666667  1.0 -1.0 -1.0;
 0.0  1.0 -1.66667   2.0 -2.0 -2.0;
 0.0  0.0  3.33333  -2.0  4.0 11.0;
 0.0  0.0 20.0     -24.0 24.0 24.0;
 0.0  0.0  9.66667  -9.0 10.0 11.0]
[1.0  0.0  0.0   0.6 -0.2   1.2;
 0.0  1.0  0.0   1.0  0.0   3.5;
 0.0  0.0  1.0  -0.6  1.2   3.3;
 0.0  0.0  0.0 -12.0  0.0 -42.0;
 0.0  0.0  0.0  -3.2 -1.6 -20.9]
[1.0  0.0  0.0  0.0  -0.2  -0.9;
 0.0  1.0  0.0  0.0   0.0   0.0;
 0.0  0.0  1.0  0.0   1.2   5.4;
 0.0  0.0  0.0  1.0  -0.0   3.5;
 0.0  0.0  0.0  0.0  -1.6  -9.7]
[1.0  0.0  0.0  0.0  0.0   0.3125;
 0.0  1.0  0.0  0.0  0.0   0.0;
 0.0  0.0  1.0  0.0  0.0  -1.875;
 0.0  0.0  0.0  1.0  0.0   3.5;
 0.0  0.0  0.0  0.0  1.0   6.0625]
[0.3125, 0.0, -1.875, 3.5, 6.0625]

●ピボット選択

ところで、ガウスの消去法とガウス・ジョルダン法は、拡大係数行列の対角成分 zs[i, i] が 0 になると計算できなくなります。たとえば、次の連立方程式は解くことができますが、zs[0, 0] が 0 になっているため、プログラムの実行中に値が無限大 (Inf) になってしまいます。

      2y + 4z = 14   [[0.  2.  4. 14.]
  x +  y +  z = 10    [1.  1.  1. 10.]
 4x + 2y + 6z = 38    [4.  2.  6. 38.]]

    連立方程式        拡大係数行列 zs

このとき、拡大係数行列の行を交換すると連立一次方程式が解ける場合があります。また、zs[i, i] が 0 に等しくなくても 0 に近い値だと解の誤差が大きくなります。そこで、zs[i, i] の絶対値がなるべく大きくなるように行を交換します。これを「部分ピボット選択」といいます。なお、参考文献『Ｃ言語による最新アルゴリズム事典』によると、『行だけでなく列も交換する完全ピボット選択という方法もあるが、通常は部分ピボット選択で十分である』 とのことです。

プログラムは次のようになります。

リスト : 部分ピボット選択

# ピボット選択
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
        println(zs)    # debug
    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]
end

# ガウス・ジョルダン法
function gauss_jordan(xs::Matrix{Float64}, ys::Vector{Float64})
    zs = hcat(xs, ys)
    n = size(xs, 1)
    for i = 1 : n
        # ピボット選択
        select_pivot(zs, i)
        temp = zs[i, i]
        zs[i, i:end] /= temp
        for j = 1 : n
            if i != j
                temp = zs[j, i]
                for k = 1 : n + 1
                    zs[j, k] -= zs[i, k] * temp
                end
            end
        end
        println(zs)   # debug
    end
    zs[:, n + 1]
end

ピボット選択は関数 select_pivot() で行います。abs.(xs[i:, i]) で絶対値を求め、その中から最大値を findmax() で求めます。返り値は最大値と位置を格納したタプルです。実際の位置はそれに i - 1 を加算したものになります。これを変数 k にセットします。k と i が異なる値であれば i 行と k 行と交換します。

あとは select_pivot() を呼び出すだけです。それでは実際に試してみましょう。

リスト : 簡単なテスト

println(gauss([0. 2. 4.; 1. 1. 1.; 4. 2. 6.], [14., 10., 38.]))
println(gauss_jordan([0. 2. 4.; 1. 1. 1.; 4. 2. 6.], [14., 10., 38.]))

[4.0  2.0  6.0 38.0;
 0.0  0.5 -0.5  0.5;
 0.0  2.0  4.0 14.0]
[4.0  2.0  6.0 38.0;
 0.0  2.0  4.0 14.0;
 0.0  0.0 -1.5 -3.0]
[5.0, 3.0, 2.0]

[1.0  0.5  1.5  9.5;
 0.0  0.5 -0.5  0.5;
 0.0  2.0  4.0 14.0]
[1.0  0.0  0.5  6.0;
 0.0  1.0  2.0  7.0;
 0.0  0.0 -1.5 -3.0]
[1.0  0.0  0.0  5.0;
 0.0  1.0  0.0  3.0;
 0.0  0.0  1.0  2.0]
[5.0, 3.0, 2.0]

正しく解を求めることができました。

●逆行列

連立方程式 \(Ax = b\) を解く場合、\(A\) の逆行列 \(A^{-1}\) がわかれば、\(x = A^{-1} b\) で解を求めることができます。実をいうと、逆行列はガウス・ジョルダン法を使って簡単に求めることができるのです。次の図を見てください。

逆行列の定義 : A * A-1 = I

[ a11 a12 a13;       [ b11 b12 b13;       [ 1 0 0;
  a21 a22 a23;   *     b21 b22 b23;    =    0 1 0;
  a31 a32 a33 ]        b31 b32 b33 ]        0 0 1 ]

   係数行列               逆行列         単位行列 I

上記の定義により逆行列 A-1 は次の方程式を解けば求めることができる

A * [b11, b21, b31] = [1, 0, 0]
A * [b12, b22, b32] = [0, 1, 0]
A * [b13, b23, b33] = [0, 0, 1]

これは次の拡大係数行列にガウス・ジョルダン法を適用すること同じ

[ a11 a12 a13  1  0  0;
  a21 a22 a23  0  1  0;
  a31 a32 a33  0  0  1 ]

拡大係数行列の右半分に逆行列が求まる

プログラムも簡単です。次のリストを見てください。

リスト : ガウス・ジョルダン法で逆行列を求める

function matinv(xs::Matrix{Float64})
    n = size(xs, 1)
    # 拡大係数行列
    zs = hcat(xs, diagm(0 => ones(n)))
    for i = 1 : n
        # ピボット選択
        select_pivot(zs, i)
        temp = zs[i, i]
        zs[i, i:end] /= temp
        for j = 1 : n
            if i != j
                temp = zs[j, i]
                for k = 1 : 2n
                    zs[j, k] -= zs[i, k] * temp
                end
            end
        end
    end
    zs[:, n + 1 : end]
end

単位行列は Matrix{Float64}(I, n, n) や diagm(0 => ones(n)) で作成することができます。Diagonal(ones(n)) でも作成できますが、hcat() で xs と連結すると疎行列になるので、今回は diagm() を使いました。あとはガウス・ジョルダン法で解を求め、最後に拡大係数行列の右半分を返します。それでは実行してみましょう。

リスト : 簡単なテスト

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.]
for a = [a1, a2, a3, a4, a5]
    b = matinv(a)
    println(a)
    println(b)
    println(a * b)
    println("----------")
end

[1.0 1.0;
 2.0 4.0]
[ 2.0 -0.5;
 -1.0  0.5]
[1.0 0.0;
 0.0 1.0]
----------
[1.0 1.0 1.0;
 2.0 4.0 6.0;
 2.0 0.0 4.0]
[1.33333  -0.333333  0.166667;
 0.333333  0.166667 -0.333333;
 -0.666667 0.166667  0.166667]
[1.0 2.77556e-17 2.77556e-17;
 0.0 1.0         1.11022e-16;
 0.0 0.0         1.0]
----------
[ 1.0 1.0  1.0 1.0;
 -1.0 1.0 -1.0 1.0;
  8.0 4.0  2.0 1.0;
 -8.0 4.0 -2.0 1.0]
[-0.166667  0.166667  0.0833333 -0.0833333;
 -0.166667 -0.166667  0.166667   0.166667;
  0.666667 -0.666667 -0.0833333  0.0833333;
  0.666667  0.666667 -0.166667  -0.166667]
[ 1.0         0.0          0.0        -2.77556e-17;
  0.0         1.0          0.0         2.77556e-17;
 -1.11022e-16 0.0          1.0        -1.11022e-16;
   0.0       -1.11022e-16 -8.32667e-17 1.0]
----------
[ 1.0 -1.0 1.0 -1.0 1.0;
 12.0 -6.0 2.0  0.0 0.0;
  1.0  1.0 1.0  1.0 1.0;
 12.0  6.0 2.0  0.0 0.0;
  4.0  3.0 2.0  1.0 0.0]
[-0.0625  0.0416667    0.0625  0.0833333 -0.125;
  0.0    -0.0833333    0.0     0.0833333  0.0;
  0.375  -3.46945e-18 -0.375  -0.25       0.75;
 -0.5     0.0833333    0.5    -0.0833333 -5.55112e-17;
  0.1875 -0.0416667    0.8125  0.166667  -0.625]
[ 1.0          2.08167e-17  0.0         2.77556e-17 1.11022e-16;
  1.11022e-16  1.0         -1.11022e-16 1.11022e-16 0.0;
  0.0          6.93889e-18  1.0         5.55112e-17 1.11022e-16;
  1.11022e-16  2.08167e-17 -1.11022e-16 1.0         0.0;
 -5.55112e-17 -1.38778e-17  0.0         5.55112e-17 1.0]
----------
[0.0 2.0 4.0;
 1.0 1.0 1.0;
 4.0 2.0 6.0]
[-0.333333  0.333333  0.166667;
  0.166667  1.33333  -0.333333;
  0.166667 -0.666667  0.166667]
[1.0          0.0  0.0;
 2.77556e-17  1.0  2.77556e-17;
 1.11022e-16  0.0  1.0]

行列と逆行列の積を計算すると、浮動小数点数の計算による誤差がありますが、単位行列になることがわかります。

ところで、実際に連立一次方程式を解くとき、逆行列を使用することはほとんどありません。わざわざ逆行列を求めるよりも、ガウスの消去法を使ったほうがより速く解を求めることができるからです。同じ係数行列の連立方程式を何度も解く場合でも、逆行列を求めるのではなく、次に説明する LU 分解を使うことが多いようです。

●参考文献・URL

1. 奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
2. ガウスの消去法 - Wikipedia