線形代数#

using LinearAlgebra

ベクトル#

Juliaは列ベクトルが基本．

\[ \mathbf{a}=[a_1, a_2, \ldots, a_n]^\top \]
a = [1, 2, 3]

3-element Vector{Int64}:
 1
 2
 3

行列#

\[\begin{split} \begin{aligned} \mathbf{A} =\left[\begin{array}{ccc} a_{11} & \cdots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{1n} & \cdots & a_{nn} \end{array}\right] \end{aligned} \in \mathbb{R}^n \end{split}\]

行列\(\mathbf{A}\)の\((i, j)\)成分を\(\mathbf{A}_{ij}=a_{ij}\)とする．

A = []

Any[]

要素ごとの演算#

行列の行・列ごとの正規化#

シミュレーションにおいてニューロン間の重み行列を行あるいは列ごとに正規化 (weight normalization)する場合がある．これは各ニューロンへの入力の大きさを同じにする働きや重みの発散を防ぐ役割がある．以下では行ごとの和を1にする．

W = rand(3,3)

3×3 Matrix{Float64}:
 0.208063  0.593138  0.381483
 0.319687  0.351225  0.394025
 0.142155  0.120187  0.499655

Wnormed = W ./ sum(W, dims=1)

3×3 Matrix{Float64}:
 0.310586  0.557172  0.299164
 0.477212  0.329928  0.309
 0.212202  0.1129    0.391837

println(sum(Wnormed, dims=1))

[1.0 0.9999999999999999 1.0]

行列の結合 (concatenate)#

行列の結合はMATLABに近い形式で行うことができる．まず，2つの行列A, Bを用意する．

A = [1 2; 3 4]

2×2 Matrix{Int64}:
 1  2
 3  4

B = [4 5 6; 7 8 9]

2×3 Matrix{Int64}:
 4  5  6
 7  8  9

水平結合 (Horizontal concatenation)#

hcatを使うやり方と，[ ]を使うやり方がある．

H1 = hcat(A,B)

2×5 Matrix{Int64}:
 1  2  4  5  6
 3  4  7  8  9

H2 = [A B]

2×5 Matrix{Int64}:
 1  2  4  5  6
 3  4  7  8  9

なお，MATLABのように次のようにすると正しく結合はされない．

H3 = [A, B]

2-element Vector{Matrix{Int64}}:
 [1 2; 3 4]
 [4 5 6; 7 8 9]

垂直結合 (Vertical concatenation)#

V1 = vcat(A, B')

5×2 Matrix{Int64}:
 1  2
 3  4
 4  7
 5  8
 6  9

V2 = [A; B']

5×2 Matrix{Int64}:
 1  2
 3  4
 4  7
 5  8
 6  9

[V2 [A;B']]

5×4 Matrix{Int64}:
 1  2  1  2
 3  4  3  4
 4  7  4  7
 5  8  5  8
 6  9  6  9

配列に新しい軸を追加#

要はnumpyでのA[None, :]やA[np.newaxis, :]のようなことがしたい場合．やや面倒だが，reshapeを使うか，[CartesianIndex()]を用いる．

v = rand(3)

3-element Vector{Float64}:
 0.08078390066519736
 0.21880164577595063
 0.38638066402308535

newaxis = [CartesianIndex()]
v1 = v[newaxis, :]

1×3 Matrix{Float64}:
 0.0807839  0.218802  0.386381

単位行列#

I

UniformScaling{Bool}
true*I

I(3)

3×3 Diagonal{Bool, Vector{Bool}}:
 1  ⋅  ⋅
 ⋅  1  ⋅
 ⋅  ⋅  1

対角行列#

線形行列方程式#

\[ \mathbf{A}\mathbf{x}=\mathbf{b} \]

\(\mathbf{A}\)が正則の場合，逆行列が存在し，

\[ \mathbf{x}=\mathbf{A}^{-1}\mathbf{b} \]
A = rand(2,2)
b = rand(2)

2-element Vector{Float64}:
 0.6642499593758039
 0.8357523767362349

x = inv(A) * b

2-element Vector{Float64}:
 -0.5413347924076956
  1.368260994277419

Juliaではバックスラッシュ演算子 \を用いることで明示的に逆行列を計算せずに解を求めることができる．

x = A \ b

2-element Vector{Float64}:
 -0.5413347924076958
  1.3682609942774193

Roth’s column lemma#

Roth’s column lemmaは，例えば，\(A, B, C\)が与えられていて，\(X\)を未知とするときの方程式 \(AXB = C\)を考えると，この方程式は

\[ (B^\top \otimes A)\text{vec}(X) = \text{vec}(AXB)=\text{vec}(C) \]

の形に書き下すことができる，というものである．\(\text{vec}(\cdot)\)はvec作用素（行列を列ベクトル化する作用素）である．vec(X) = vcat(X...)で実現できる．Roth’s column lemmaを用いれば，\(AXB = C\)の解は

\[ X = \text{vec}^{-1}\left((B^\top \otimes A)^{-1}\text{vec}(C)\right) \]

として得られる．ただし，\(\text{vec}(\cdot)^{-1}\)は列ベクトルを行列に戻す作用素(inverse of the vectorization operator)である．reshape()で実現できる．2つの作用素をまとめると，

\[\begin{split} \begin{align} \text{vec} &: R^{m\times n}\to R^{mn}\\ \text{vec}^{−1} &: R^{mn}\to R^{m×n} \end{align} \end{split}\]

であり，\(\text{vec}^{−1}\left(\text{vec}(X)\right)=X\ (\text{for all}\ X\in R^{m\times n})，\text{vec}\left(\text{vec}^{−1}(x)\right)=x\ (\text{for all}\ x \in R^{mn})\)となる．

using LinearAlgebra, Kronecker, Random

m = 4
A = randn(m, m)
B = randn(m, m)
C = convert(Array{Float64}, reshape(1:16, (m, m)))

4×4 Matrix{Float64}:
 1.0  5.0   9.0  13.0
 2.0  6.0  10.0  14.0
 3.0  7.0  11.0  15.0
 4.0  8.0  12.0  16.0

X = reshape((B' ⊗ A) \ vec(C), (m, m))

4×4 Matrix{Float64}:
   148.347  -74.8909  -392.707   133.712
   113.584  -44.518   -278.883    99.2009
 -1055.49   507.847   2751.61   -945.167
  -544.143  269.801   1432.13   -489.247

A * X * B

4×4 Matrix{Float64}:
 1.0  5.0   9.0  13.0
 2.0  6.0  10.0  14.0
 3.0  7.0  11.0  15.0
 4.0  8.0  12.0  16.0

配列の1次元化#

配列を一次元化(flatten)する方法．まずは3次元配列を作成する．

B = rand(2, 2, 2)

2×2×2 Array{Float64, 3}:
[:, :, 1] =
 0.681655  0.487036
 0.513269  0.277588

[:, :, 2] =
 0.262146  0.991841
 0.966683  0.713706

用意されているflattenを素直に用いると次のようになる．

import Base.Iterators: flatten
collect(flatten(B))

8-element Vector{Float64}:
 0.6816546797447299
 0.5132688758361276
 0.4870356238821668
 0.27758760252750925
 0.26214565051023875
 0.966683275932893
 0.9918412964721871
 0.713705675056993

ただし，単にB[:]とするだけでもよい．

B[:]

8-element Vector{Float64}:
 0.6816546797447299
 0.5132688758361276
 0.4870356238821668
 0.27758760252750925
 0.26214565051023875
 0.966683275932893
 0.9918412964721871
 0.713705675056993

reshapeにおける残りの次元の指定#

numpyにおいては(2, 3, 5)次元の配列に対し，reshape(-1, 5)を行うと(6, 5)次元の配列となった．これと同様なことは，Juliaでは:を使うことで実装できる．

a = rand(2,3,5)
b = reshape(a, (:, 5))

6×5 Matrix{Float64}:
 0.918377  0.403922  0.621976   0.53253   0.857224
 0.18758   0.804664  0.552209   0.434259  0.834326
 0.738768  0.261772  0.621643   0.310129  0.993374
 0.501559  0.313444  0.564713   0.935804  0.726968
 0.726934  0.272811  0.123096   0.104811  0.605749
 0.526304  0.912738  0.0609884  0.291267  0.487083

Array{Array{Float64, x},1}をArray{Float64, x+1}に変換#

numpyではarray([matrix for i in range()])などを用いると，1次元配列のリストを2次元配列に変換できた．Juliaでも同様にする場合はhcat(...)やcat(...)を用いる．...はsplat operatorと呼ばれる．

A1 = [i*rand(3) for i=1:5]

println("Type : ", typeof(A1))
println("Size : ", size(A1))

Type : Vector{Vector{Float64}}
Size : (5,)

A2 = hcat(A1...)'

println("Type : ", typeof(A2))
println("Size : ", size(A2))

Type : Adjoint{Float64, Matrix{Float64}}
Size : (5, 3)

以下は多次元配列の場合．cat(...)で配列を結合し，permitedimsで転置する．

B1 = [i*rand(3, 4, 5) for i=1:6]

println("Type : ", typeof(B1))
println("Size : ", size(B1))

Type : Vector{Array{Float64, 3}}
Size : (6,)

B2 = permutedims(cat(B1..., dims=4), [4, 1, 2, 3])

println("Type : ", typeof(B2))
println("Size : ", size(B2))

Type : Array{Float64, 4}
Size : (6, 3, 4, 5)