【Julia】基礎: Matrix

0. はじめに

本稿では Julia 言語の行列型のプログラムについて説明します．
参考という意味でスター * を節名につけましょう．

環 $R$ 上の $n \times n$ 行列 ($n \in \mathbb{N} $) 全体 $\mathrm{M}(n,R) $ は行列和と行列積に関して環をなします．

1. 行列計算を始める

行列計算を始めてみましょう．

1.1 undef

まず初期化されていない $2 \times 2 $ の行列を構成します．

n = 2
display(Matrix{Float64}(undef,n,n))

出力例:

2×2 Matrix{Float64}:
 2.23609e-314  2.23609e-314
 2.23609e-314  2.23609e-314

1.2 初期化*

初期化する方法には nothing や missing を使う方法があります．¹

X1 = Array{Union{Nothing, Int}}(nothing, 2, 3)
X2 = Array{Union{Missing, Int}}(missing, 2, 3)
display(X1)
display(X2)
X2[1,2] = 1
display(X2)

2×3 Matrix{Union{Nothing, Int64}}:
 nothing  nothing  nothing
 nothing  nothing  nothing
2×3 Matrix{Union{Missing, Int64}}:
 missing  missing  missing
 missing  missing  missing
2×3 Matrix{Union{Missing, Int64}}:
 missing  1         missing
 missing   missing  missing

1.2 ゼロ行列

Julia でゼロ行列を生成するには zeros を使います．それが Z1 です．

Z2 と Z3 は Julia の標準の行列の記法であり，Z3 は MATLAB の公式ドキュメントにも記載があります．²

Z4 は推奨されないと思われますが同じ結果を与えます．

# n = 2
Z1 = zeros(n,n)
Z2 = [
    0 0
    0 0
]
Z3 = [0 0; 0 0]
Z4 = 0 * Matrix{Float64}(undef,n,n)

println(Z1 == Z2 == Z3 == Z4)

コンマとセミコロンを併用するとエラーが返されます．

# A1 = [1, 0; 0, 1]
# ERROR: syntax: unexpected semicolon in array expression around ...(absolute path)

1.3 単位行列

単位行列の生成は現在 LinearAlgebra パッケージの仕様が推奨される事が多いです．³

• 整数型単位行列
• 浮動小数点型 〜
• LinearAlgebra の I で Float64 型として生成
• Base.one で Z1 から生成 ⁴

下3つは等しいです．

using LinearAlgebra
I2_int = [
    1 0
    0 1
]
I2_float = [
    1.0 0.0
    0.0 1.0
]
I2_alg = Matrix{Float64}(I,n,n)
I2_one = one(Z1)
display(I2_int)
display(I2_float == I2_alg == I2_one)

2×2 Matrix{Int64}:
 1  0
 0  1
true

1.4 すべての成分が1の行列

Base.ones ですべての成分が環 $R$ の単位元 $1_{R}$ であるような行列を生成できます．$n \ne 1$ ならばそれは行列環の単位元 $1_{\mathrm{M}(n,R)}$ とは異なります．$n = 1$ ならばそれらは等しくなります．

# n = 2
display(ones(n,n))
println(ones(1,1) == one(ones(1,1)))

2×2 Matrix{Float64}:
 1.0  1.0
 1.0  1.0
true

2. 行列演算

行列環の多項演算を見ていきましょう．

2.1 行列和

2.1.1 加法

行列環は環，つまり加法 $+$ に関して加群です．

#%%addition
A1 = [1 0; 0 1]
A2 = [2 0; 0 2]
A3 = A1 + A2
A4 = A2 + A1# inverse order
display(A3)
println(A3 == A4)# commutative

2×2 Matrix{Int64}:
 3  0
 0  3
true

要素全てに同じ数を足すには以下のようにします．
ドット記法 .+ ではなく + を使うとエラーが出ます．

M = [1 0; 0 1]
r = 1
# M + r
# ERROR: MethodError: no method matching...
Mplusr = M .+ r
display(Mplusr)
display(Mplusr == M + r * ones(Int,n,n))

2×2 Matrix{Int64}:
 2  1
 1  2
true

2.1.2 減法

差は - で計算します．

# A1 = [1 0; 0 1]
# A2 = [2 0; 0 2]
@show A1 - A2, A2 - A1
# (A1 - A2, A2 - A1) = ([-1 0; 0 -1], [1 0; 0 1])

2.2 行列積

2.2.1 行列積

行列積は $n \ne 1$ の場合一般には非可換です．
ただし $R$ が可換環で $n = 1$ ならば行列積も可換になります．

#%%multiplication
A1 = [1 1; 0 1]
A2 = [1 0; 1 1]
A3 = A1 * A2
A4 = A2 * A1# noncommutative
@show A1, A2, A3, A4
begin end

(A1, A2, A3, A4) = ([1 1; 0 1], [1 0; 1 1], [2 1; 1 1], [1 1; 1 2])

2.2.2 除法*

除法は / で計算します．
特に I に対する除法は積に関する逆元を返します．

# using LinearAlgebra
A = [1 1; 1 0]
B = I / A
display(B)
println(A * B == B * A == I)

2×2 Matrix{Float64}:
 0.0   1.0
 1.0  -1.0
true

Backslach \ で連立方程式を解く事ができます．

#%% Ax = b
A = [1 0; 0 -1]
b = [0, 1]
x = A \ b
@show x

x = [0.0, -1.0]

この場合 x または b は A の固有ベクトルになっています．
A は物理数学で Pauli 行列と呼ばれるものの3番めです．

2.3 冪

---

冪はハット ^ を使います．

---

MATLAB を基本にしているためです．⁵
Python の ** は使えません．

#%%power
A = [1 1; 1 0]
N = 5
for j = 1:N
    @show (A^j)# Fibonacci numbers
end

A ^ j = [1 1; 1 0]
A ^ j = [2 1; 1 1]
A ^ j = [3 2; 2 1]
A ^ j = [5 3; 3 2]
A ^ j = [8 5; 5 3]

ゼロ乗は単位行列を返します．

display(A^0)
#=
2×2 Matrix{Int64}:
 1  0
 0  1
=#

ドット記法 .^ を使うと要素ごとのべき乗になります．⁶

2.4 逆行列

---

逆行列とは行列環の積に関する逆元です．マグマ (集合と2項演算) の可逆元をとることを1項演算と解釈することはしばしばあります．

---

逆行列の生成には inv() を使います．⁷

負冪も使えます．

A = [1 1; 1 0]
@show inv(A), A^-1

(inv(A), A ^ -1) = ([0.0 1.0; 1.0 -1.0], [0.0 1.0; 1.0 -1.0])

(参考) 一応時間を測ってみましたが大差はありませんでした．Julia は出力に 1 ms - 100 ms かかるのが普通なので出力も考えませんでした．負冪が速いことも多々ありました．
メモリ使用量は 10000×10000 行列で 約 1.5 GiB (ギビバイト) でした

N1 = 10000
@time inv(rand(N1,N1))
@time (rand(N1,N1))^-1

 22.878215 seconds (930.38 k allocations: 1.540 GiB, 0.19% gc time, 1.69% compilation time)
 25.530230 seconds (521 allocations: 1.495 GiB, 1.16% gc time, 0.01% compilation time)

エラー処理については ⁸ が詳しいです．

2.5 対合

環には対合の構造を考えることができます．⁹

2.5.1 転置

転置行列は関数 transpose() で得られます．
これは実行列環の対合です．

#%%transpose()
A0100 = [0 1; 0 0]
display(transpose(A0100))

2×2 transpose(::Matrix{Int64}) with eltype Int64:
 0  0
 1  0

2.5.2 エルミート共役

エルミート共役は single quotation mark ' で得られます．
これは複素行列環の対合です．

#%%Hermitian conjugate
A0im00 = [0 im; 0 0]
display(A0im00')

2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 0+0im  0+0im
 0-1im  0+0im

2.6 スカラー倍

行列環は加法に関して加群である．
(単位的) 環 $R$ による「(左)作用」は以下で定まる．

rA :=\left( rA_{jk}\right) _{j,k} \quad (A\in \mathrm{M}\left( n,R\right) ,\,r\in R).

$(\mathrm{M}\left( n,R\right);R)$ は環上の加群をなし，この演算はスカラー倍と呼ばれる．
スカラー倍は Julia では * で計算できる．

# n = 2
X = Matrix{Float64}(I,n,n)
a = 2
display(a * X)

2×2 Matrix{Float64}:
 2.0  0.0
 0.0  2.0

3. 対角化

---

MATLAB 同様 Julia では対角化を高速に計算することができます．

---

3.1. 2準位系の理論計算

2準位系の固有エネルギーおよびその摂動論を考えてみましょう．

H=\begin{pmatrix}
E_{1} & \epsilon  \\
\epsilon  & E_{2}
\end{pmatrix}

を Hamiltonian とします．

$\mathbf{1}_{2}$ を2次単位行列とし固有方程式

\det \left( H-\lambda \mathbf{1}_{2}\right) 
= \left( E_{1}-\lambda \right) \left( E_{2}-\lambda \right) -\epsilon ^{2}

を解きます．

右辺は $ \lambda ^{2}-\left( E_{1}+E_{2}\right) \lambda +E_{1}E_{2} -\epsilon ^{2} $ だから2次方程式の解の公式から

\lambda =\dfrac{E_{1}+E_{2}\pm \sqrt{\left( E_{1}+E_{2}\right) ^{2}-4(E_{1}E_{2} -\epsilon ^{2})}}{2} .

$\left( E_{1}+E_{2}\right) ^{2}-4(E_{1}E_{2} -\epsilon ^{2}) = \left( E_{1}-E_{2}\right) ^{2}+4\epsilon ^{2}$ なので $E_{\mathrm{mean}} := \dfrac{E_{1}+E_{2}}{2}$, $\Delta E :=\left| E_{2}-E_{1}\right| $ とおくと

\lambda =E_{\mathrm{mean}}\pm \dfrac{\Delta E}{2}\sqrt{1 + \left(2 \frac{\epsilon}{\Delta E}\right)^{2}} .

$\frac{\epsilon}{\Delta E} \ll 1 $ とすると摂動論が帰結するはずで，$\sqrt{1 + \left(2 \frac{\epsilon}{\Delta E}\right)^{2}} \simeq 1 + \frac{1}{2} \left(2 \frac{\epsilon}{\Delta E}\right)^{2} $ なので $E_{1} \le E_{2} $ ならば

\lambda =
E_{2}+\left( \dfrac{\epsilon }{\Delta E}\right) ^{2} ,\,
E_{1}-\left( \dfrac{\epsilon }{\Delta E}\right) ^{2} .

---

物理的な解釈として 1s 軌道の 2原子合成が考えられます．

$(1,2)$ 成分と $(2,1)$ 成分が同じ $\epsilon$, (特に同符号です) Hamiltonian は自己共役作用素の条件を満たしています．

相互作用系ではエネルギー準位が「反発」し基底状態は $E_{1}$ より小さいので，2原子がそこへ落ち込んだ分子軌道が実現するでしょう．

つまり2つの原子の間に「引力」が働くでしょう．

---

一方もし $(2,1)$ 成分が逆符号で $ - \epsilon$ だったら「斥力」が働くわけですから分子軌道は実現しません．そもそも この 行列 (または作用素) は自己共役でないので量子力学の公理としての Hamiltonian の資格を備えていません．

固有ベクトルの理論計算は省略します．

3.2. 2準位系の数値計算

少し長めに理論計算してきましたが Julia で数値計算すれば (病的でない数値であれば) 8行で終わってしまいます．

3.2.1. eigen()

eigen() で固有値と固有ベクトルが得られます．

E1 = 0
E2 = 1
ϵ  = 0.01
H = [
    E1 ϵ
    ϵ  E2
]
eigen(H)

Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
 -9.999000199950014e-5
  1.0000999900019996
vectors:
2×2 Matrix{Float64}:
 -0.99995    0.0099985
  0.0099985  0.99995

理論計算と一致していますね．

3.2.2. eigvals()

eigvals() で固有値だけが得られます．

display(eigvals(H))

2-element Vector{Float64}:
 -9.999000199950014e-5
  1.0000999900019996

3.2.3. eigvecs()

また eigvecs() で固有ベクトルだけが得られます．

display(eigvecs(H))

2×2 Matrix{Float64}:
 -0.99995    0.0099985
  0.0099985  0.99995

4. おわりに

Julia の行列計算について説明してきました．
MATLAB の文法を引き継いでいる Julia では Python よりも数倍から数百倍高速な線型代数計算が可能です．
文法も簡単なのでぜひ覚えて使っていきましょう．

1. Arrays · The Julia Language ↩

2. 行列および配列 - MATLAB & Simulink - MathWorks 日本 ↩

3. 昔は eye を使っていたようです． ↩

4. https://docs.julialang.org/en/v1/base/numbers/#Base.one ↩

5. べき乗と指数 - MATLAB & Simulink - MathWorks 日本 ↩

6. 要素単位のべき乗 - MATLAB power .^ - MathWorks 日本 ↩

7. Linear Algebra · The Julia Language ↩

8. https://stanford.edu/class/engr108/lectures/julia_inverses_slides.pdf ↩

9. *-algebra - Wikipedia ↩