🌱

Juliaで線形代数：ことはじめ

2022/01/29に公開

Julia

線形代数というと, ベクトル空間の公理系を出発点として議論を進めなければいけない気がする. しかし, 世の中の入門的な教科書とか講義はそのようにはなっておらず, 数字の羅列で遊ぶ程度のものである. というのも, 同じ次元のベクトル空間は同型であることが証明されているから, ベクトルは数字の羅列と考えて良いし, 線形写像も表現行列を考えればよい. 次元とか基底という概念を理解しておいて欲しいという気持ちはあるが, 数ベクトルや行列の扱いに慣れるための第一歩として, Juliaを使う教育があってもよいのではないだろうか.

LinearAlgebra.jlはJuliaに最初からインストールされている線形代数のパッケージである. Fortranで書かれた人類史上最強[要出典]の線形代数パッケージLAPACKを呼び出して利用しているので,　固有値の計算はもちろんのこと, ノルムや内積, クロス積なども簡単に計算できる. この記事ではベクトル・行列・行列式の性質や定理について, Juliaによる具体的な計算例を挙げながら解説する.

パッケージ
ベクトル・行列・行列式の使い方
シュヴァルツの不等式
三角不等式
次元定理
回転群
行列式の性質
Trotter分解
エルミート行列, 固有値問題, そして量子力学へ
参考文献

パッケージ

LinearAlgebra.jlはJuliaに最初からインストールされているので, ノートブック上でusing LinearAlgebraを宣言するだけでよい.

パッケージの読み込み

using LinearAlgebra

ベクトル・行列・行列式の使い方

簡単な説明に留めるので, 慣れたら公式のドキュメントやこちらの記事などを参照されたい. まず, 縦ベクトルはu = [1,2,3]やu = vcat(1,2,3)のように宣言する. 横ベクトルはu = [1 2 3]やu = hcat(1,2,3)のように宣言する(正確には1行n列の行列). ベクトルの和 $\pmb{u}+\pmb{v}$ や差 $\pmb{u}-\pmb{v}$ , ノルム $|\pmb{u}|$ , 内積 $\pmb{u}\cdot\pmb{v}$ やクロス積 $\pmb{u}\times\pmb{v}$ が利用できる. なお, クロス積は3次元ベクトルでしか使えないので注意. これは, 一般の次元を考えるときはクロス積の一部の性質を諦めることになるためである. 詳しくはグラスマン代数を参照されたい.

ベクトルについて

u = [1,2,3]
v = [2,2,2]

println("u          = ", u)
println("v          = ", v)
println("norm(u)    = ", norm(u))
println("norm(v)    = ", norm(v))
println("u + v      = ", u + v)
println("u - v      = ", u - v)
println("dot(u,v)   = ", dot(u,v))
println("cross(u,v) = ", cross(u,v))

出力

u          = [1, 2, 3]
v          = [2, 2, 2]
norm(u)    = 3.7416573867739413
norm(v)    = 3.4641016151377544
u + v      = [3, 4, 5]
u - v      = [-1, 0, 1]
dot(u,v)   = 12
cross(u,v) = [-2, 4, -2]

行列については, A = [1 2 1; 0 1 1; 1 -1 -2]のように宣言できる. プログラミング言語では;を改行として解釈することが多いので, 改行でも代用できる. transpose(A)またはA'で転置行列 $\boldsymbol{A}^\mathrm{T}$ を表す. ただし, 複素行列に対してはA'はエルミート転置 $\boldsymbol{A}^\dagger$ となるので, 単純な転置を取りたい場合はtranspose(A)を使用する. また, トレース $\mathrm{Tr}~\boldsymbol{A}$ はtr(A), 行列式 $\mathrm{det}~\boldsymbol{A}$ はdet(A), ランクは $\mathrm{rank}~\boldsymbol{A}$ はrank(A)である. その他にも, 行列ノルムnorm(A)や行列指数関数exp(A)も使用できる.

行列について

A = [1  2  1
     0  1  1
     1 -1 -2]

println("A       = ", A)
println("A'      = ", A')
println("tr(A)   = ", tr(A))
println("det(A)  = ", det(A))
println("rank(A) = ", rank(A))

出力

A       = [1 2 1; 0 1 1; 1 -1 -2]
A'      = [1 0 1; 2 1 -1; 1 1 -2]
tr(A)   = 0
det(A)  = 0.0
rank(A) = 2

シュヴァルツの不等式

詳細はWikipediaまたは斎藤正彦『基礎数学1 線型代数入門』(東京大学出版会, 1966) p.6 を参照されたい. 最初なので, 2つのベクトル

\pmb{u} = \begin{pmatrix} 1 \\ 2 \\ 3 \\ \end{pmatrix} , \pmb{v} = \begin{pmatrix} 2 \\ 2 \\ 2 \\ \end{pmatrix}

を具体例にして丁寧に説明する. これらの内積は

\pmb{u}\cdot\pmb{v} = 1\cdot2 + 2\cdot2 + 3\cdot2 = 12

である. また, それぞれのノルムは

\begin{align} |\pmb{u}| &= \sqrt{1\cdot1 + 2\cdot2 + 3\cdot3} = \sqrt{14} = 3.74... \\ |\pmb{v}| &= \sqrt{2\cdot2 + 2\cdot2 + 2\cdot2} = \sqrt{12} = 3.46... \end{align}

であるから, これらの積は

|\pmb{u}||\pmb{v}| = 12.96...

である. 従って, シュヴァルツの不等式

\pmb{u}\cdot\pmb{v} \leq |\pmb{v}||\pmb{u}|

が確かに成り立っている. Juliaで計算した例：

シュヴァルツの不等式

u = [1,2,3]
v = [2,2,2]

println("|u・v| = ", dot(u,v))
println("|u||v| = ", norm(u)*norm(v))

出力

|u・v| = 12
|u||v| = 12.961481396815719

線形従属のときに等式が成り立つことがわかる.

線形従属のとき

u = [3,3,3]
v = [2,2,2]

println("|u・v| = ", dot(u,v))
println("|u||v| = ", norm(u)*norm(v))

出力

|u・v| = 18
|u||v| = 18.0

三角不等式

詳細はWikipediaまたは斎藤正彦『基礎数学1 線型代数入門』(東京大学出版会, 1966) p.6 を参照されたい.

三角不等式

u = [1,2,3]
v = [2,2,2]

println("|u＋v|   = ", norm(u+v))
println("|u|＋|v| = ", norm(u)+norm(v))

出力

|u＋v|   = 7.0710678118654755
|u|＋|v| = 7.205759001911696

シュヴァルツの不等式と同様に, 線形従属のときに等式が成り立つことがわかる.

線形従属のとき

u = [3,3,3]
v = [2,2,2]

println("|u＋v|   = ", norm(u+v))
println("|u|＋|v| = ", norm(u)+norm(v))

出力

|u＋v|   = 8.660254037844387
|u|＋|v| = 8.660254037844386

次元定理

次元定理についてはWikipediaやこちらの記事を参照されたい. まず, 斎藤正彦『基礎数学4 線型代数演習』(東京大学出版会, 1985) p.3 §1の問題より

1 次の三つのベクトルは線型独立か.
1)

\begin{pmatrix} 1 \\ 0 \\ 1 \\ \end{pmatrix} , \begin{pmatrix} 2 \\ 1 \\ -1 \\ \end{pmatrix} , \begin{pmatrix} 1 \\ 1 \\ -2 \\ \end{pmatrix}

を考えよう. 線形独立の定義はWikipediaまたは斎藤正彦『基礎数学1 線型代数入門』(東京大学出版会, 1966) p.5 を参照されたい. 既に見たように, シュヴァルツの不等式や三角不等式を使えば, 2つのベクトルが線形独立かどうかは簡単に判定できる. 3つ以上のベクトルについても, 行列式やランクを使えば判定できる. 斎藤正彦『基礎数学1 線型代数入門』(東京大学出版会, 1966) 第1章§5, [5.1]の1)より, 3つのベクトルを並べた行列の行列式が0なので線形従属であることがわかる. ランクが2<3であることからも同じ結論が得られる. さらに, $\mathrm{rank}~A + \mathrm{dim}(\mathrm{Ker}~A)=3$ であるから, 次元定理が成り立っていることも確認できる.

次元定理の例

A = [1  2  1;
     0  1  1;
     1 -1 -2]

println("det(A)      = ", det(A))
println("rank(A)     = ", rank(A))
println("dim(Ker(A)) = ", size(nullspace(A))[2])

出力

det(A)      = 0.0
rank(A)     = 2
dim(Ker(A)) = 1

回転群

定義はWikipediaを参照されたい. ここでは $\theta=\pi/2$ とした $2\times2$ の回転行列を使って回転群を作る.

回転行列

R(θ) = [cos(θ) -sin(θ);
         sin(θ) cos(θ)]

A = R(π/2)

出力

     6.12323e-17  -1.0
     1.0           6.12323e-17

対角要素がきれいに0にならないので, 整数型の行列に変換した.

型変換

A = Int.(floor.((R(π/2))))

出力

     0  -1
     1   0

この行列は自身に4回作用させると元に戻ることが分かった. 従って, 位数4の回転群であることが予想される.（90°ずつ回転させているので当然ではある.）また, A⁴は単位行列なので, 行列の積について考えれば単位元である.

println("A¹ = ", A^1)
println("A² = ", A^2)
println("A³ = ", A^3)
println("A⁴ = ", A^4)
println("A⁵ = ", A^5)

出力

A¹ = [0 -1; 1 0]
A² = [-1 0; 0 -1]
A³ = [0 1; -1 0]
A⁴ = [1 0; 0 1]
A⁵ = [0 -1; 1 0]

次にように積表を生成し, 位数4の群であることが網羅的に確かめられる（単位元はA⁴）. 辞書(Dict)をうまく使うとよい.

群であることの確認

G = Dict()

G[A^1] = "A¹"
G[A^2] = "A²"
G[A^3] = "A³"
G[A^4] = "A⁴"

for j in 1:4
    for i in 1:4
        println(G[A^i], G[A^j], " = ", G[A^i*A^j])
    end
end

出力

A¹A¹ = A²
A²A¹ = A³
A³A¹ = A⁴
A⁴A¹ = A¹
A¹A² = A³
A²A² = A⁴
A³A² = A¹
A⁴A² = A²
A¹A³ = A⁴
A²A³ = A¹
A³A³ = A²
A⁴A³ = A³
A¹A⁴ = A¹
A²A⁴ = A²
A³A⁴ = A³
A⁴A⁴ = A⁴

また, 全ての元の行列式が1であるから, SO(2)の部分群である.

SO(2)の部分群であることの確認

println("det(A¹) = ", det(A^1))
println("det(A²) = ", det(A^2))
println("det(A³) = ", det(A^3))
println("det(A⁴) = ", det(A^4))
println("det(A⁵) = ", det(A^5))

出力

det(A¹) = 1.0
det(A²) = 1.0
det(A³) = 1.0
det(A⁴) = 1.0
det(A⁵) = 1.0

行列式の性質

詳細はWikipediaまたは斎藤正彦『基礎数学1 線型代数入門』(東京大学出版会, 1966) 定理2.2, 定理2.3, 系2.4, 定理2.7 を参照されたい.

1つの列を $c$ 倍すると行列式も $c$ 倍になる.（多重線形性, 定理2.2）
列を入れ換えると符号が逆になる（交代性, 定理2.3）
同じ列を含む場合, 必ず0. (系2.4)
2つの $n$ 次行列の積の行列式は, それぞれの行列式の積に等しい. (定理2.7)

を順に確認する.

行列式の性質

A = [1  2  2;
     0  2  1;
     1 -1 -2]

B = [2  2  2;
     0  2  1;
     2 -1 -2]

C = [2  1  2;
     2  0  1;
    -1  1 -2]

D = [1  1  2;
     0  0  1;
     1  1 -2]

println("det(A)   = ", det(A))
println("det(B)   = ", det(B)  , "\t← 定理2.2")
println("det(C)   = ", det(C)  , "\t\t← 定理2.3")
println("det(D)   = ", det(D)  , "\t\t← 系2.4")
println("det(A*B) = ", det(A*C), "\t← 定理2.7")

出力

det(A)   = -5.0
det(B)   = -10.0	← 定理2.2
det(C)   = 5.0		← 定理2.3
det(D)   = 0.0		← 系2.4
det(A*B) = -25.0	← 定理2.7

Trotter分解

Trotter分解を演算子の指数関数へ拡張したSuzuki-Trotter分解は量子モンテカルロ法の基礎となる重要な定理の一つである. Torotter分解の詳細はWikipedia を参照されたい. 通常の指数関数では

\exp(a+b) = \exp(a)\exp(b)

が成り立つのに対し, 行列指数関数では

\exp(A+B) = \lim_{n\rightarrow\infty} \left[ \exp(A/n)\exp(B/n) \right]^n

が成り立つ. ここでは, $n$ を大きくするにつれて左辺に収束する様子を確認する. 見ての通り収束は速くないので, 実用上は改良版が使われている.

Trotter分解

using Printf
MatrixPrint(n, M) = println( "n = $n\t\t" * @sprintf("[%.6f\t", M[1,1]) * @sprintf("%.6f;\t", M[2,1]) * @sprintf("%.6f\t", M[1,2]) * @sprintf("%.6f]", M[2,2]) )

A = [1 1;
     1 1]

B = [1 0;
     1 1]

for n in [1,2,3,4,5,10,20,30,40,50,100,200,300,400,500]
    MatrixPrint(n, (exp(A/n)*exp(B/n))^n)
end

MatrixPrint("∞", exp(A+B))

出力

n = 1		[20.085537	20.085537;	8.683628	11.401909]
n = 2		[18.416237	20.167777;	9.686841	13.572816]
n = 3		[17.701504	20.195096;	9.915277	14.396411]
n = 4		[17.318629	20.205912;	9.999229	14.818822]
n = 5		[17.081891	20.211157;	10.038831	15.074124]
n = 10		[16.595136	20.218397;	10.092389	15.585897]
n = 20		[16.346316	20.220252;	10.105916	15.841020]
n = 30		[16.262726	20.220597;	10.108427	15.925779]
n = 40		[16.220824	20.220718;	10.109306	15.968091]
n = 50		[16.195649	20.220775;	10.109713	15.993455]
n = 100		[16.145230	20.220849;	10.110256	16.044127]
n = 200		[16.119986	20.220868;	10.110392	16.069434]
n = 300		[16.111567	20.220871;	10.110417	16.077865]
n = 400		[16.107356	20.220873;	10.110426	16.082080]
n = 500		[16.104829	20.220873;	10.110430	16.084608]
n = ∞		[16.094720	20.220874;	10.110437	16.094720]

エルミート行列, 固有値問題, そして量子力学へ

エルミート行列の性質は, その固有値や固有ベクトルを考えてこそ面白みがある. 固有値問題については, 次の記事では解説するので, エルミート行列も次の記事で解説することにしよう.

参考文献

https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/
斎藤正彦『基礎数学1 線型代数入門』(東京大学出版会, 1966)

この記事の元になったノートブックはGistにアップロードしてある.

パッケージ

ベクトル・行列・行列式の使い方

シュヴァルツの不等式

三角不等式

次元定理

回転群

行列式の性質

Trotter分解

エルミート行列, 固有値問題, そして量子力学へ

参考文献

Discussion