ジュリアで線形代数パッケージを使用する方法
概要
Juliaは、MATLABレベルの線形代数をサポートしている。むしろMATLABよりも進化した、直感的で美しい構文を見ると、Juliaが作られた時点でよく設計されていたと感じられる1。
コード
julia> A = [
1 0 3
0 5 1
3 1 9
]
3×3 Matrix{Int64}:
1 0 3
0 5 1
3 1 9
見ての通り、行列を定義する段階で既に直感的で便利だ。ここで、普通に必要とされるいくつかの関数について学ぶ。小見出しに関連記事をリンクし、別の説明は省略する。
トレース tr()
julia> tr(A)
15
行列式 det()
julia> det(A)
-1.000000000000003
逆行列 inv()
julia> inv(A)
3×3 Matrix{Float64}:
-44.0 -3.0 15.0
-3.0 6.10623e-16 1.0
15.0 1.0 -5.0
julia> round.(Int64, inv(A))
3×3 Matrix{Int64}:
-44 -3 15
-3 0 1
15 1 -5
対角行列と対角成分 diag(), diagm()
julia> diag(A)
3-element Vector{Int64}:
1
5
9
julia> diagm([1,5,9])
3×3 Matrix{Int64}:
1 0 0
0 5 0
0 0 9
ノルム norm()
julia> norm(A, 1)
23.0
固有値 eigvals()
julia> eigvals(A)
3-element Vector{Float64}:
-0.020282065792505244
4.846013411157458
10.174268654635046
julia> eigvecs(A)
3×3 Matrix{Float64}:
-0.944804 0.117887 0.305692
-0.0640048 -0.981459 0.180669
0.321322 0.151132 0.934832
julia> eigmax(A)
10.174268654635046
行列分解 factorize()
julia> factorize(A)
BunchKaufman{Float64, Matrix{Float64}}
D factor:
3×3 Tridiagonal{Float64, Vector{Float64}}:
-0.0227273 0.0 ⋅
0.0 4.88889 0.0
⋅ 0.0 9.0
U factor:
3×3 UnitUpperTriangular{Float64, Matrix{Float64}}:
1.0 -0.0681818 0.333333
⋅ 1.0 0.111111
⋅ ⋅ 1.0
permutation:
3-element Vector{Int64}:
1
2
3
julia> svd(A)
SVD{Float64, Float64, Matrix{Float64}}
U factor:
3×3 Matrix{Float64}:
-0.305692 0.117887 -0.944804
-0.180669 -0.981459 -0.0640048
-0.934832 0.151132 0.321322
singular values:
3-element Vector{Float64}:
10.174268654635044
4.846013411157461
0.02028206579250516
Vt factor:
3×3 Matrix{Float64}:
-0.305692 -0.180669 -0.934832
0.117887 -0.981459 0.151132
0.944804 0.0640048 -0.321322
行列代数カテゴリの行列分解を参照せよ。行列の形に応じて適切な分解法を自動的に選んで分解してくれる。もちろん、条件を満たすなら、具体的な分解関数を直接使っても良い。
行列の操作
julia> B = [
1 0 1
1 1 0
2 1 1
]
3×3 Matrix{Int64}:
1 0 1
1 1 0
2 1 1
julia> A + B
3×3 Matrix{Int64}:
2 0 4
1 6 1
5 2 10
julia> A - B
3×3 Matrix{Int64}:
0 0 2
-1 4 1
1 0 8
julia> A * B
3×3 Matrix{Int64}:
7 3 4
7 6 1
22 10 12
julia> A .* B
3×3 Matrix{Int64}:
1 0 3
0 5 0
6 1 9
julia> B / A
3×3 Matrix{Float64}:
-29.0 -2.0 10.0
-47.0 -3.0 16.0
-76.0 -5.0 26.0
julia> B * inv(A)
3×3 Matrix{Float64}:
-29.0 -2.0 10.0
-47.0 -3.0 16.0
-76.0 -5.0 26.0
julia> A / B
ERROR: SingularException(3)
我々が考える常識的な操作は全て通用する。割り算は、当然乗算の逆元である逆行列を掛けるのと同じで、Bのように逆行列が存在しない場合は、シンギュラー例外をレイズする。
ブロック行列 []
他の言語と比べて、ブロック行列を非常に便利に作ることができる。
julia> [A B]
3×6 Matrix{Int64}:
1 0 3 1 0 1
0 5 1 1 1 0
3 1 9 2 1 1
julia> [A;B]
6×3 Matrix{Int64}:
1 0 3
0 5 1
3 1 9
1 0 1
1 1 0
2 1 1
julia> [A,B]
2-element Vector{Matrix{Int64}}:
[1 0 3; 0 5 1; 3 1 9]
[1 0 1; 1 1 0; 2 1 1]
二つの行列の間にスペースを置くと横に積み上げ、セミコロンを置くと縦に積み上げる。カンマは行列を積み上げるわけではなく、一般的に配列で使用していた構文そのままで、行列の配列になる。
全体のコード
複素行列や内積に関連する内容は省略したが、全体のコードには含まれている。
using LinearAlgebra
A = [
1 0 3
0 5 1
3 1 9
]
tr(A)
det(A)
inv(A)
round.(Int64, inv(A))
diag(A)
diagm([1,5,9])
norm(A, 1)
eigvals(A)
eigvecs(A)
eigmax(A)
factorize(A)
svd(A)
B = [
1 0 1
1 1 0
2 1 1
]
det(B)
rank(B)
eigvals(B)
Symmetric(B) # |> issymmetric
transpose(B)
B'
C = [
im im 1
2 im 0
im 1 2
]
C'
B'B
x = [1,2,3]
y = [0,1,2]
x'y
A + B
A - B
A * B
A .* B
B / A
B * inv(A)
[A B]
[A;B]
[A,B]
x' * y
y * x'
環境
- OS: Windows
- julia: v1.7.0