💯

Juliaは行列をシンプルに書ける言語です!

2023/01/06に公開

例えば, ij列目の要素がi+10jであるような3\times4の行列Aを宣言するとき, 次のようにたったの1行で内包的に宣言することができます.

julia> A = [i + 10j for i=1:3, j=1:4]
3×4 Matrix{Int64}:
 11  21  31  41
 12  22  32  42
 13  23  33  43

非常に便利ですね. 他の例も見ていきましょう.

ベクトルの外延的・内包的宣言

まずベクトルの宣言を見ていきましょう. 外延的には次のように書けます.

julia> x = [1,2,3,4]
4-element Vector{Int64}:
 1
 2
 3
 4

内包的には次のように書けます.

julia> x = [i for i=1:4]
4-element Vector{Int64}:
 1
 2
 3
 4

行列の外延的・内包的宣言

外延的には次のように書けます. インデントの仕方は自由です. 改行の代わりにセミコロン;を使うこともできます.

julia> A = [5  4  3  2  1
            4  4  3  2  1
            3  3  3  2  1
            2  2  2  2  1
            1  1  1  1  1]
5×5 Matrix{Int64}:
 5  4  3  2  1
 4  4  3  2  1
 3  3  3  2  1
 2  2  2  2  1
 1  1  1  1  1

内包的には, ベクトルの場合と同様にforを使います. iが行, jが列の添え字です.

julia> A = [5 - max(i,j) + 1 for i=1:5, j=1:5]
5×5 Matrix{Int64}:
 5  4  3  2  1
 4  4  3  2  1
 3  3  3  2  1
 2  2  2  2  1
 1  1  1  1  1

このように, 1行で行列を宣言することができます.

フランク行列の固有値問題を解く

100\times100フランク行列の固有値問題Ax=\lambda xを解きます. REPL上またはノートブック上でusing LinearAlgebraを宣言すればeigen(A)の戻り値として固有値λと固有ベクトルxが求められます.

using LinearAlgebra
n = 100
A = [n - max(i,j) + 1 for i=1:n, j=1:n]
λ,x = eigen(A)

例えば最小固有値は

julia> λ[1]
0.2500610827206968

であり, Fortran/LAPACKで求めた値0.25006108272069305と一致します. 最後の2桁がズレている理由は, スレッド数, BLASのディストリビューション, LAPACKのバージョン, 収束閾値などによって差が出るためです.

おわりに

外延的な記法ではどの言語も同じようなものですが, 内包的な記法だとJuliaのシンプルさが際立ちますね. 2次元配列を宣言してメモリを確保し, 2重ループを回して各要素に代入する, なんてことをいちいちやる時代は終わりました. 非常に便利なので積極的に使っていきましょう!

Discussion