【Julia】基礎: Matrix

<!-- 【Julia】基礎: Matrix -->

# 0. はじめに

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

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

# 1. 行列計算を始める

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

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

```julia
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` を使う方法があります．[^1]

[^1]: [Arrays · The Julia Language](https://docs.julialang.org/en/v1/base/arrays/#Base.Matrix)

<!-- ```julia:initialize -->
```julia
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 の公式ドキュメントにも記載があります．[^2]

[^2]: [行列および配列 - MATLAB & Simulink - MathWorks 日本](https://jp.mathworks.com/help/matlab/learn_matlab/matrices-and-arrays.html)

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

<!-- ```julia:zero matrix -->
```julia
# 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` パッケージの仕様が推奨される事が多いです．[^eye]

[^eye]: 昔は `eye` を使っていたようです．

* 整数型単位行列
* 浮動小数点型 〜
* `LinearAlgebra` の `I` で `Float64` 型として生成
* `Base.one` で `Z1` から生成 [^3]

[^3]: https://docs.julialang.org/en/v1/base/numbers/#Base.one

下3つは等しいです．

```julia
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$ ならばそれらは等しくなります．

```julia
# 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 加法

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

```julia
#%%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
```

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

```julia
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 減法

差は `-` で計算します．

```julia
# 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$ ならば行列積も可換になります．

```julia
#%%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` に対する除法は積に関する逆元を返します．

```julia
# 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 `\` で連立方程式を解く事ができます．

```julia
#%% 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 を基本にしているためです．[^power]
Python の `**` は使えません．

[^power]: [べき乗と指数 - MATLAB & Simulink - MathWorks 日本](https://jp.mathworks.com/help/matlab/math/powers-and-exponentials.html)

```julia
#%%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]
```

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

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

ドット記法 `.^` を使うと要素ごとのべき乗になります．[^power-element]

[^power-element]: [要素単位のべき乗 - MATLAB power .^ - MathWorks 日本](https://jp.mathworks.com/help/matlab/ref/power.html)

# 2.4 逆行列

---

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

---

逆行列の生成には `inv()` を使います．[^inv]

[^inv]: [Linear Algebra · The Julia Language](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/)

負冪も使えます．

```julia
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 (ギビバイト) でした

```julia
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)
```

エラー処理については [^stanford] が詳しいです．

[^stanford]: https://stanford.edu/class/engr108/lectures/julia_inverses_slides.pdf


# 2.5 対合

環には対合の構造を考えることができます．[^star-alg]

[^star-alg]: [*-algebra - Wikipedia](https://en.wikipedia.org/wiki/*-algebra)

# 2.5.1 転置

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

```julia
#%%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 `'` で得られます．
これは複素行列環の対合です．


```julia
#%%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$ による「(左)作用」は以下で定まる．

```math
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 では `*` で計算できる．

```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準位系の固有エネルギーおよびその摂動論を考えてみましょう．

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

を Hamiltonian とします．

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

```math
\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次方程式の解の公式から

```math
\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| $ とおくと


```math
\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} $ ならば

```math
\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()` で固有値と固有ベクトルが得られます．

```julia
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()` で固有値だけが得られます．

```julia
display(eigvals(H))
```

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

# 3.2.3. `eigvecs()`

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


```julia
display(eigvecs(H))
```

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

<!-- 対角化の数値計算の理論については割愛します． -->


# 4. おわりに

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

<!-- [Arrays · The Julia Language](https://docs.julialang.org/en/v1/base/arrays/#Base.Matrix) -->