M.Hiroi's Home Page

Julia Language Programming

お気楽 Julia プログラミング超入門

[ Home | Light | Julia ]

線形代数編

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

julia> using LinearAlgebra

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

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

julia> a + b
3-element Array{Float64,1}:
 0.9456292131041033
 1.6182823861431488
 1.1222049978783177

julia> a - b
3-element Array{Float64,1}:
 -0.7585074013902424
 -0.2661712969240555
  0.7286559878143941

julia> a - a
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

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

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

julia> dot(a, b)
0.8988190903480853

julia> norm(a)
1.14988101105807

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

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

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

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

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

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

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

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

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

julia> a / 10
2×3 Array{Float64,2}:
 0.0186229  0.00531527  0.0526371
 0.0462431  0.00348494  0.00496523
julia> c = rand(3, 2)
3×2 Array{Float64,2}:
 0.0876367  0.424824
 0.865688   0.593409
 0.91316    0.827241

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

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

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

julia> d * e == d
true

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

julia> e * d == d
true

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

julia> e * e == e
true

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

julia> f ^ 2
2×2 Array{Int64,2}:
 2  1
 1  1

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

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

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

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

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

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

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

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

julia> a1 \ [4, 6]
2-element Array{Float64,1}:
  5.0
 -2.0
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 Array{Int64,2}:
 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)

●特殊な行列とその性質

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

julia> a'
2×2 LinearAlgebra.Adjoint{Int64,Array{Int64,2}}:
 1  3
 2  4

julia> transpose(a)
2×2 Transpose{Int64,Array{Int64,2}}:
 1  3
 2  4

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

julia> b'
3×3 LinearAlgebra.Adjoint{Int64,Array{Int64,2}}:
 0  1  2
 3  4  5
 6  7  8

julia> transpose(b)
3×3 Transpose{Int64,Array{Int64,2}}:
 0  1  2
 3  4  5
 6  7  8

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

julia> c'
2×3 LinearAlgebra.Adjoint{Int64,Array{Int64,2}}:
 0  1  2
 3  4  5

julia> transpose(c)
2×3 Transpose{Int64,Array{Int64,2}}:
 0  1  2
 3  4  5
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 Array{Int64,2}:
 5  6
 7  8

julia> (a * d)'
2×2 Adjoint{Int64,Array{Int64,2}}:
 19  43
 22  50

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

julia> inv(a')
2×2 Array{Float64,2}:
 -2.0   1.5
  1.0  -0.5

julia> inv(a)'
2×2 Adjoint{Float64,Array{Float64,2}}:
 -2.0   1.5
  1.0  -0.5
julia> a = [1 3; 3 2]
2×2 Array{Int64,2}:
 1  3
 3  2

julia> a' == a
true

julia> b = [1 4 5; 4 2 6; 5 6 3]
3×3 Array{Int64,2}:
 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 Array{Int64,2}:
 1  0  0
 0  2  0
 0  0  3

julia> c' == c
true
julia> a = [1 2 3; 4 5 6; 7 8 10]
3×3 Array{Int64,2}:
 1  2   3
 4  5   6
 7  8  10

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

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

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

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

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

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

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

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

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

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

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

julia> diagm(-1 => [1, 2, 3])
4×4 Array{Int64,2}:
 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 Array{Int64,2}:
 1  0  0  0
 1  1  0  0
 0  2  1  0
 0  0  3  1

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

julia> det(a)
6

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

julia> prod(diag(a))
6

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

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

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

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

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

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

julia> eigvals(a)  # 行列 a の固有値を求める
3-element Array{Int64,1}:
 1
 2
 3
julia> a = [1 2 3; 4 5 6; 7 8 10]
3×3 Array{Int64,2}:
 1  2   3
 4  5   6
 7  8  10

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

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

julia> la'
3×3 Adjoint{Int64,LowerTriangular{Int64,Array{Int64,2}}}:
 1  4   7
 0  5   8
 0  0  10

julia> ua'
3×3 Adjoint{Int64,UpperTriangular{Int64,Array{Int64,2}}}:
 1  0   0
 2  5   0
 3  6  10

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

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

julia> det(la)
50

julia> det(ua)
50

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

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

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

julia> eigvals(ua)
3-element Array{Int64,1}:
  1
  5
 10
回転行列 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 Array{Int64,2}:
 0  1
 1  0

julia> a'
2×2 Adjoint{Int64,Array{Int64,2}}:
 0  1
 1  0

julia> det(a)
-1.0

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

julia> b = [0 1 0; 1 0 0; 0 0 1]
3×3 Array{Int64,2}:
 0  1  0
 1  0  0
 0  0  1

julia> b'
3×3 Adjoint{Int64,Array{Int64,2}}:
 0  1  0
 1  0  0
 0  0  1

julia> det(b)
-1.0

julia> inv(b)
3×3 Array{Float64,2}:
 -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 Array{Int64,2}:
 1  2  3
 4  5  6
 7  8  9

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

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

●連立一次方程式の解法

今回は連立一次方程式を解く基本的なアルゴリズムを取り上げます。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] の絶対値がなるべく大きくなるように行を交換します。これを「部分ピボット選択」といいます。なお、参考文献 1 によると、『行だけでなく列も交換する完全ピボット選択という方法もあるが、通常は部分ピボット選択で十分である』 とのことです。

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

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

# ピボット選択
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-1b で解を求めることができます。実をいうと、逆行列はガウス・ジョルダン法を使って簡単に求めることができるのです。次の図を見てください。

逆行列の定義 : 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. 奥村晴彦,『C言語による最新アルゴリズム事典』, 技術評論社, 1991
  2. ガウスの消去法 - Wikipedia

Copyright (C) 2018 Makoto Hiroi
All rights reserved.

[ Home | Light | Julia ]