読者です 読者をやめる 読者になる 読者になる

行列演算

Matlabにそっくりな文法なので、いろいろ行列演算をしてみました。「#」は、Matlabの「%」やJavaの「//」に相当するコメント記号で、そこから行末までが無視されます。

julia> A = randi(6, 3, 3)   # 1〜6の乱数で3x3の行列を作る
3x3 Int64 Array:
 5  6  5
 1  3  3
 6  6  6

julia> diag(A)   # 対角成分だけ取り出す
[5, 3, 6]

julia> norm(A)   # ノルムの計算
14.504793991339222

julia> rank(A)   # ランクの計算
3

julia> inv(A)   # 逆行列の計算
3x3 Float64 Array:
  0.0  -0.5   0.25    
  1.0   0.0  -0.833333
 -1.0   0.5   0.75    

julia> A * inv(A)   # 逆行列をかけてみる
3x3 Float64 Array:
 1.0  0.0   0.0        
 0.0  1.0  -4.44089e-16
 0.0  0.0   1.0        

julia> A \ A   # バックスラッシュ(円記号)でも同じことができる
3x3 Float64 Array:
  1.0   0.0  0.0
  0.0   1.0  0.0
 -0.0  -0.0  1.0

julia> typeof(A)   # Aの型を調べると、要素がInt64の2次元配列だとわかる
Array{Int64,2}

julia> mean(A,1)   # 1次元目で平均を計算
1x3 Float64 Array:
 4.0  5.0  4.66667

julia> mean(A,2)   # 2次元目で平均を計算
3x1 Float64 Array:
 5.33333
 2.33333
 6.0    

julia> A[:,1]   # 1列目だけを抽出
3x1 Int64 Array:
 5
 1
 6

julia> mean(A[:,1])   # 1列目だけの平均を計算したいのだけど(行列だから)うまくいかない
no method mean(Array{Int64,2},)
 in method_missing at /Users/marui/julia/j/base.j:58

julia> squeeze(A[:,1])   # squeezeしてやるとベクトルになった
[5, 1, 6]

julia> mean(squeeze(A[:,1]))   # ベクトルだと平均が計算できる
4.0

julia> help(svd)   # 特異値分解の使い方を調べる
(U, S, V) = svd(A)
  Compute the SVD of A

julia> (U,S,V) = svd(A)   # 特異値分解してみる
(
3x3 Float64 Array:
 -0.638685   0.040433  -0.768405
 -0.283767  -0.940608   0.186368
 -0.715233   0.337078   0.612226,

[14.5048, 1.52162, 0.543704],
3x3 Float64 Array:
 -0.535588  -0.618748  -0.574715
  0.843852  -0.365898  -0.39247 
  0.032553  -0.695176   0.718102)

julia> U * diagm(S) * V   # Sを対角成分に持つ行列を作り、Aに戻す
3x3 Float64 Array:
 5.0  6.0  5.0
 1.0  3.0  3.0
 6.0  6.0  6.0

最初から基本的な行列計算は色々できますね。特異値分解もできるので、主成分分析の実装もすぐにできそうです。