Julia には LinearAlgebra という「標準パッケージ」があります。
標準とはいっても Core や Base パッケージのように組み込みではないのでインストールが必要です。Julia の対話モードから “]” を押して、パッケージモードへ移行し次のコマンドを実行します。
pkg> add LinearAlgebra
線形代数の勉強お勧め動画 (予備校のノリで学ぶ線形代数入門)
- 概観&ベクトル
- 行列 (演算規則)
- 一次変換
- 一次独立と一次従属
- 連立方程式 (掃き出し法)
- 不定と不能
- 階数 (ランク)
- 行列式 (定義と性質)
- 余因子展開
- 逆行列:定義
- 逆行列:掃き出し法
- 固有値・固有ベクトル
- 対角化 (重解がない場合)
- 対角化 (重解がある場合)
ベクトル
ベクトルの初期化
横ベクトル は v1 = [0.5, 1.1, -2.0] 、縦ベクトルは v2 = [3.3; 1.7; 0.2] のようにします。要素の型を指定する場合は次のようにします。
julia> v2 = Int32[1, 2, 4]
3-element Array{Int32,1}:
1
2
4
要素の値が未定義のベクトルは v3 = Vector{Float64}(undef, 3) のようにします。undef の代わりに nothing や missing も使えます。
julia> v3 = Vector{Float64}(undef, 3)
3-element Array{Float64,1}:
1.390662761733164e-309
1.390662959311913e-309
1.39066270598319e-309
ベクトルの絶対値 (ノルム)
ベクトルの絶対値ですが、abs(v) とやるとエラーになります。norm(v) が正しいです。
julia> v2
2-element Vector{Int64}:
1
2
julia> norm(v2)
2.23606797749979
julia> sqrt(5.0) # 1.0^2 + 2.0^2 = 5.0
2.23606797749979
内積と外積
次にベクトルの内積ですが、v1 * v2 とやるとエラーになります。dot(v1, v2) が正しいです。また、外積は cross(v1, v2) で計算できます。
using LinearAlgebra
v1 = [0.5, 1.1, -2.0]
v2 = [3.3; 1.7; 0.2]
# 絶対値
println(norm(v1))
# 内積
println(dot(v1, v2))
# 外積
println(cross(v1, v2))
実行結果
2.3366642891095846
3.12
[3.62, -6.699999999999999, -2.78]
行列
初期化
行列を初期化する場合、行と行の間の区切りに “;” を使います。同じ行の要素の区切りは ” ” を使います。
julia> A = [3 1 3; 0 -1 1; 5 2 -2]
3×3 Array{Int64,2}:
3 1 3
0 -1 1
5 2 -2
行列の型を特に指定しないと64ビット整数または64ビット浮動小数点数になりますが、型を指定したい場合は次のよう行列の前に型名を指定します。
julia> A = Int32[3 1 3; 0 -1 1; 5 2 -2]
3×3 Array{Int32,2}:
3 1 3
0 -1 1
5 2 -2
要素を後から決める場合は、次のように要素が undef の行列を作ります。
julia> A = Matrix{Float64}(undef, 3, 3)
3×3 Array{Float64,2}:
5.0e-323 6.0e-323 6.4e-323
1.0e-323 1.0e-323 1.0e-323
1.39066e-309 1.39066e-309 1.39066e-309
要素に値を代入するには次のように行います。ここで先頭の要素のインデックスは [0, 0]ではなく [1, 1]
になります。
julia> A[1, 1] = 0.0
0.0
julia> A
3×3 Array{Float64,2}:
0.0 6.0e-323 6.4e-323
1.0e-323 1.0e-323 1.0e-323
1.39066e-309 1.39066e-309 1.39066e-309
インデックスの指定ですが、最初の方が行番号、後の方が列番号になります。
julia> A[2, 1] = 5.0
5.0
julia> A
3×3 Array{Float64,2}:
0.0 6.0e-323 6.4e-323
5.0 1.0e-323 1.0e-323
1.39066e-309 1.39066e-309 1.39066e-309
単位行列を作るには次のようにします。
julia> E = Matrix{Int32}(I, 3, 3)
3×3 Array{Int32,2}:
1 0 0
0 1 0
0 0 1
要素がすべて 0 の行列は次のように作ります。
julia> Z = zeros(Int32, 3, 3)
3×3 Array{Int32,2}:
0 0 0
0 0 0
0 0 0
要素がすべて 1 の行列も次のようにして作れます。
julia> O = ones(Int32, 3, 3)
3×3 Array{Int32,2}:
1 1 1
1 1 1
1 1 1
乱数で初期化することもできます。下の例ではシード値等はデフォルト値を使っていますが、もっと細かい指定もできます。
julia> R =rand(3, 3)
3×3 Array{Float64,2}:
0.588734 0.609012 0.680494
0.212418 0.1448 0.87533
0.380341 0.219114 0.570458
行列の演算
行列の演算は足し算、引き算、掛け算がありますが、割り算はありません。ベクトルの掛け算 (内積) では “*” は使えませんでしたが、行列では “*” が使えます。
また、ベクトルも行列の一種なので、行列とベクトルの掛け算も可能です。
足し算、引き算は同じ次元の行列同士でのみ可能です。掛け算は最初の行列の列数と後の行列の行数が同じである必要があります。
次のサンプルは行列の簡単な演算例です。
# 行列の演算
using LinearAlgebra
# 定義
A = Int32[5 9 1; 3 -1 -2; 1 0 8]
println("A = ", A)
B = Int32[1 5 -1; -8 4 4; 2 4 -8]
println("B = ", B)
C = Int32[-5 5 -9; 1 -6 6; -1 1 7]
println("C = ", C)
E = Matrix{Int32}(I, 3, 3)
println("E = ", E)
x = Int32[2; 1; -6]
println("x = ", x)
# 計算
println("A + B = ", A + B)
println("A - B = ", A - B)
println("A * B = ", A * B)
println("A * x = ", A * x)
println("+ 結合法則 ", (A + B) + C == A + (B + C))
println("* 結合法則 ", (A * B) * C == A * (B * C))
println("+ 交換法則 ", A + B == B + A)
println("* 交換法則 ", A * B == B * A) # 等しくない
println("* 分配法則 (スカラー) ", 2.0 * (A + B) == 2.0 * A + 2.0 * B)
println("* 分配法則 (行列) ", C * (A + B) == C * A + C * B)
println("* 単位行列 (交換法則) ", E * A == A * E)
# 正方行列以外の掛け算
D = Int32[2 1 0; 0 1 2]
println("D = ", D)
F = Int32[1 2; 3 4; 5 6]
println("F = ", F)
println("D * F = ", D * F)
# べき乗
println("A^2 = ", A^2)
# スカラー倍
println("2 * A = ", 2 * A)0 * A)
実行結果
A = Int32[5 9 1; 3 -1 -2; 1 0 8]
B = Int32[1 5 -1; -8 4 4; 2 4 -8]
C = Int32[-5 5 -9; 1 -6 6; -1 1 7]
E = Int32[1 0 0; 0 1 0; 0 0 1]
x = Int32[2, 1, -6]
A + B = Int32[6 14 0; -5 3 2; 3 4 0]
A - B = Int32[4 4 2; 11 -5 -6; -1 -4 16]
A * B = Int32[-65 65 23; 7 3 9; 17 37 -65]
A * x = Int32[13, 17, -46]
+ 結合法則 true
* 結合法則 true
+ 交換法則 true
* 交換法則 false
* 分配法則 (スカラー) true
* 分配法則 (行列) true
* 単位行列 (交換法則) true
D = Int32[2 1 0; 0 1 2]
F = Int32[1 2; 3 4; 5 6]
D * F = Int32[5 8; 13 16]
A^2 = Int32[53 36 -5; 10 28 -11; 13 9 65]
2 * A = [10 18 2; 6 -2 -4; 2 0 16]
一次変換
行列にベクトルを掛けるということは、そのベクトルを別のベクトルに変換するということです。この変換のことを一次変換と言います。(一次式だけからなるため)
単位行列にベクトルを掛けてもは元のベクトルのままです。
julia> E = Matrix{Float32}(I, 3, 3)
3×3 Matrix{Float32}:
1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
julia> x = [2.0; 5.0; 1.0]
3-element Vector{Float64}:
2.0
5.0
1.0
julia> E * x
3-element Vector{Float64}:
2.0
5.0
1.0
対角行列にベクトルを掛けると、ベクトルの方向は変わらず長さだけが変わります。
julia> D = Float32[2.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 3.0]
3×3 Matrix{Float32}:
2.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 3.0
julia> D * x
3-element Vector{Float64}:
4.0
5.0
3.0
次のように定義される回転行列を掛けるとベクトルの偏角が変わりベクトルを回転させることができます。
julia> theta = pi / 4.0
0.7853981633974483
julia> T = [cos(theta) -sin(theta); sin(theta) cos(theta)]
2×2 Matrix{Float64}:
0.707107 -0.707107
0.707107 0.707107
julia> v = [1.0; 1.0]
2-element Vector{Float64}:
1.0
1.0
julia> T * v
2-element Vector{Float64}:
1.1102230246251565e-16
1.414213562373095
対角行列と回転行列を組み合わせてベクトルに掛けると長さと角度を変えたベクトルに変換できます。
これらは二次元だけでなく n 次元で定義されています。
一次独立と一次従属
n 個のベクトルの集合 (v1, v2, … vn) の中の任意のベクトル vk が残りすべてのベクトルの線形結合で表すことができなければ「一次独立」と言い、そうでなければ「一次従属」と言います。
線形結合とは c1*v1 + c2*v2 + … + cn*vn のように定数 ck とベクトル vk を掛けたものの和です。
2次元で言えば、一次独立は方向の異なるベクトルどうし、一次従属は方向が同じなベクトルどうしになります・
つまり、v2 = c1 * v1 のような関係です。3次元で言えば、3つのベクトルのうち他のベクトルが作る平面上にそのベクトルが存在すれば一次従属です。
つまり、c1*v1 + c2*v2 = v3 になるような数 c1, c2 が存在すれば一次従属です。
一次独立かどうかの判別は「ランク (階数)」によって知ることができます。
つまり、n 個のベクトルの集合 (v1, v2, … vn) のランクが n ならこれらのベクトルは一次独立、n 未満なら一次従属なベクトルが存在することになります。
Julia ではランクは LinearAlgebra.rank() で得ることができます。
次の例で行列 A で定義される3つのベクトルは一次独立ですが、B のランクは3つのベクトルに対してランクが 2 なので一次従属です。
julia> using LinearAlgebra
julia> A = [1 0 0; 0 2 1; 2 1 0]
3×3 Matrix{Int64}:
1 0 0
0 2 1
2 1 0
julia> B = [2 0 0; 1 0 0; 2 1 0]
3×3 Matrix{Int64}:
2 0 0
1 0 0
2 1 0
julia> rank(A)
3
julia> rank(B)
2
連立一次方程式
連立一次方程式を解くには各式の係数と右辺のベクトルを1つにした係数行列を作って「掃き出し法」を使うと機械的に解くことができるのでよく使用されます。
しかし、Julia ではずばり連立方程式を解くための関数があるのでそれを利用します。
ところで、連立方程式は必ず解けるとは限りません。係数行列のすべての行が一次独立でないと解くことができません。
つまり、係数行列のランクが係数行列の次元より小さいと解くことができません。
# 連立方程式
using LinearAlgebra
# Ax = b を LAPACK.gesvx! を使って解く。
A = [0.5 0.4 1.2; -0.8 -03 1.5; 9.1 4.4 -2.0]
b = [0.1; 5.2; -1.0] # 注意:[0.1 5.2 -1.0] だと Matrix になってしまう。Vector でないとエラーになる。
tx = LAPACK.gesvx!(A, b)
x = tx[1] # タプルの最初の要素が答え
println(x)
# 検算
y = A * x
println(y)
実行結果
[0.8275653436592451; -1.7869373184898356; 0.33416021297192644]
[0.09999999999999998; 5.2; -0.9999999999999998]
連立方程式の解を求める関数は、行列が特殊な場合などについていろいろ用意されています。
行列式
行列式は正方行列を基にある手順で計算したスカラー値です。
2次と3次の正方行列から計算される行列式の値は簡単に計算できますが、それ以上の正方行列の場合は手計算では面倒です。
Julia ではもちろん行列式の計算も関数を呼び出すことにより簡単に計算できます。
行列式を計算する関数は det() です。
julia> using LinearAlgebra
julia> A = Int32[0 2 2; 1 -1 0; 3 1 2]
3×3 Matrix{Int32}:
0 2 2
1 -1 0
3 1 2
julia> det(A)
4.0
行列式は「余因子展開」で使用されたり、2次元空間での閉じた図形を一次変換したときの面積の倍率に対応しています。
余因子展開
余因子とは
余因子とは行列の m 行 n 列の要素を取り去った小行列の行列式にある決められた方法で符号をつけたものです。
決められた方法とは m + n が偶数なら正、奇数なら負とするものです。
次の例では 3×3 正方行列 A の 3 行 3 列 を除いた 行 1:2, 列 1:2 を取り出した行列 B の余因子行列を計算しています。余因子の符号は 3 + 3 なのでそのままになります。
julia> A = Int32[1 2 3; 4 5 6; 7 8 9]
3×3 Matrix{Int32}:
1 2 3
4 5 6
7 8 9
julia> B = A[1:2, 1:2]
2×2 Matrix{Int32}:
1 2
4 5
julia> det(B)
-3.0
余因子展開
余因子展開は行列式の計算をより低次の行列式に帰着させたものです。3次正方行列で2行2列に着目した場合で言えば、次のようになります。a21 などは行列の要素の値、A21 などは余因子です。
det(A) = a21*A21+a22*A22+a23*A23
逆行列
逆行列は元の行列と掛けると単位行列になるという行列です。Julia では関数 inv() で簡単に逆行列を計算できます。つまり、E を単位行列とすれば次のようになります。
A * inv(A) = E
実際に逆行列を計算してみます。
julia> A = [0.0 2.1 5.2; 3.4 8.2 -1.1; 7.9 -2.5 6.9]
3×3 Matrix{Float64}:
0.0 2.1 5.2
3.4 8.2 -1.1
7.9 -2.5 6.9
julia> B = inv(A)
3×3 Matrix{Float64}:
-0.120003 0.0612835 0.100207
0.071672 0.0915797 -0.0394141
0.163363 -0.0369841 0.0159172
julia> A * B
3×3 Matrix{Float64}:
1.0 2.77556e-17 -1.38778e-17
2.77556e-17 1.0 3.46945e-18
2.22045e-16 -5.55112e-17 1.0
行列の掛け算では交換法則は成り立ちませんが、逆行列との掛け算では交換法則が成り立ちます。
julia> B * A
3×3 Matrix{Float64}:
1.0 0.0 0.0
0.0 1.0 -5.55112e-17
-2.77556e-17 6.93889e-18 1.0
この行列は正確には単位行列ではありませんが、ほぼ単位行列になっています。これはコンピュータによる浮動小数点演算の誤差によるものです。ほぼ等しいかどうかは isapprox() 関数で知ることができます。
julia> A = [1.0 0.0 0.0; 0.0 1.0 -5.55112e-17; -2.77556e-17 6.93889e-18 1.0]
3×3 Matrix{Float64}:
1.0 0.0 0.0
0.0 1.0 -5.55112e-17
-2.77556e-17 6.93889e-18 1.0
julia> E = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
3×3 Matrix{Float64}:
1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
julia> isapprox(A, E)
true
(参考) A*x = B を満たすベクトル x は B \ A で求めることができます。
固有値
固有値とはベクトルを線形変換したとき、元のベクトルのスカラー倍になっているとき、そのスカラー値を固有値と言います。また、そのベクトルを固有ベクトルと言います。
Julia では固有値と固有ベクトルは eigen() 関数で計算できます。
julia> A = [2.0 2.4 1.0; -0.6 -1.7 0.0; 1.8 3.8 2.4]
3×3 Matrix{Float64}:
2.0 2.4 1.0
-0.6 -1.7 0.0
1.8 3.8 2.4
julia> eivals, eivec = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
-1.4137209983476255
0.8506178385312904
3.2631031598163345
vectors:
3×3 Matrix{Float64}:
0.353439 -0.845957 0.540171
-0.740758 0.199001 -0.0653024
0.571277 0.494727 0.839018
これの意味ですが、values が固有値 (3つある)、vectors が固有ベクトル (3つある) です。
一つ目の固有ベクトルは [0.353439; -0.740758; 0.571277] です。2つ目、3つ目も同様に3×3 行列の列になります。
実際に元の行列に固有ベクトルを掛けると、固有値 * 固有ベクトルになるか検算してみます。
julia> x = [0.353439; -0.740758; 0.571277]
3-element Vector{Float64}:
0.353439
-0.740758
0.571277
julia> A * x
3-element Vector{Float64}:
-0.4996642
1.0472252000000002
-0.8076253999999996
julia> -1.4137209983476255 * 0.353439
-0.49966413593498643
julia> -0.740758 * -1.4137209983476255
1.0472251392939904
julia> 0.571277 * -1.4137209983476255
-0.8076262907730365
対角化
対角行列とは対角成分以外 0 で、対角成分が正方行列の固有値であるような行列です。
n 次正方行列 A の n 個の一次独立な固有ベクトル x1, x2, … xn を並べた行列を P とすると行列 A は次のようにして対角化できます。
inv(P) * A * P
対角行列ではべき乗が簡単にできます。つまり、対角行列の n 乗は対角成分の n 乗になります。
これを Julia で確認してみます。
julia> A
3×3 Matrix{Float64}:
2.0 2.4 1.0
-0.6 -1.7 0.0
1.8 3.8 2.4
julia> lambda, P = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
-1.4137209983476255
0.8506178385312904
3.2631031598163345
vectors:
3×3 Matrix{Float64}:
0.353439 -0.845957 0.540171
-0.740758 0.199001 -0.0653024
0.571277 0.494727 0.839018
julia> D = inv(P) * A * P
3×3 Matrix{Float64}:
-1.41372 6.31439e-16 -2.22045e-16
-1.66533e-16 0.850618 2.77556e-16
0.0 -4.44089e-16 3.2631
julia> D ^ 2
3×3 Matrix{Float64}:
1.99861 -3.55565e-16 -4.10645e-16
9.37755e-17 0.723551 1.14179e-15
7.39557e-32 -1.82686e-15 10.6478