⇧ 前 目 次

Juliaで行列の計算

目次

なぜJuliaか

行列の計算は煩雑でかなりの頻度で途中で間違えてしまう。また、添字付きの変数の変形も間違えやすい。そこでコンピュータにやらせたい。

WolframAlphaも優秀だが入力はなかなか面倒。入力欄に何かを入れてボタンを押すものより、プログラムを組む方が柔軟に発展させられるとの予想から、Maxma, Python, R, Julia をざっとみて、たまたま Julia にひっかかったので、ちょっとメモを作っておく。新しいものにはなにかあるかもしれない。

対話モード

juliaコマンドを打つことで、対話モードになる。REPL (read-eval-print loop)と呼ばれる状態。

adachi@ebook:julia$ julia
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.4.1
 _/ |\__'_|_|_|\__'_|  |  Ubuntu ⛬  julia/1.4.1+dfsg-1
|__/                   |

julia> A=[1 1 2; 0 2 1; 0 0 3]
3×3 Array{Int64,2}:
 1  1  2
 0  2  1
 0  0  3

A=[1 1 2; 0 2 1; 0 0 3]と行列を入力すると、3x3に並べて書いてくれる

Aの値は記憶されて、Aと打てば値を3x3で出す。2とAの積は*を略して2Aと打てば各要素が2倍になり、A*Aは行列としての掛け算が実行される。

julia> A
3×3 Array{Int64,2}:
 1  1  2
 0  2  1
 0  0  3

julia> 2A
3×3 Array{Int64,2}:
 2  2  4
 0  4  2
 0  0  6

julia> A*A
3×3 Array{Int64,2}:
 1  3  9
 0  4  5
 0  0  9

julia> exit()

これはこれで便利だ。コマンドには履歴も使える。

しかし、プログラムをファイルに書いてスクリプト言語のようにも使う方がしっくりくるかな。

プログラムファイルで実行

スクリプト言語のようにファイルに記述して実行できる。

matrix001.jl

#!/usr/bin/julia
#行列の定義
A = [
1 1 2
0 2 1
0 0 3
]
#行列の計算と表示
print("A   = ")
println(A)
print("2A  = ")
println(2A)
print("A*A = ")
println(A*A)

juliaコマンドの引数にファイルを指定すれば、実行できる。#!の記述があるので./matrix001.jl でも実行できる(こちらの方法では実行フラグが立っている必要はある)。

adachi@ebook:julia$ julia matrix001.jl 
A  =  [1 1 2; 0 2 1; 0 0 3]
2A  = [2 2 4; 0 4 2; 0 0 6]
A*A = [1 3 9; 0 4 5; 0 0 9]

対話モード(REPL)と異なり、表示が1行になってしまう。

プログラムで3x3表示

プログラムでも行列の形に書き出すために、2つの方法を見つけた。

matrix002.jl

#!/usr/bin/julia
#行列の定義
A = [
   1 1 2
   0 2 1
   0 0 3
]
print("A   = ")
println(A)
#方法1
println("using display")
display(A)
println()
#方法2
println("using print_array")
Base.print_array(stdout, A);
println("")

displayはREPLと同様の表示。print_arrayは基本的にファイルへの出力を想定しているらしい。

adachi@ebook:julia$ ./matrix002.jl 
A   = [1 1 2; 0 2 1; 0 0 3]
using display
3×3 Array{Int64,2}
 1  1  2
 0  2  1
 0  0  3
using print_array
 1  1  2
 0  2  1
 0  0  3

println("")またはprintln()がないと次の出力が繋がってしまう。

ベクトルと行列

ベクトルは配列、行列は2次元の配列で表現しているらしい。ベクトルは縦ベクトル、1列の行列として扱われる。

コンマ区切りはベクトル

v=[1,2,3]
println("[1,2,3]→")
display(v)
println()

縦ベクトル

[1,2,3]→
3-element Array{Int64,1}:
 1
 2
 3

スペース区切りは1行中の項の区切り

M1l=[1 2 3]
println("[1 2 3]→")
display(M1l)
println()

1行3列の行列

[1 2 3]→
1×3 Array{Int64,2}:
 1  2  3

セミコロン区切りは行の区切り

M1c=[1;2;3]
println("[1;2;3]→")
display(M1c)
println()

書き方からは3行1列の行列だが、ベクトルと同じ扱いにも見える

[1;2;3]→
3-element Array{Int64,1}:
 1
 2
 3

改行区切りも行の区切り

M1l2=
[
 1 
 2
 3
]
println("3行に分けて記述→")
display(M1l2)
println()

;の代わりに改行でも3行1列の行列になる。

3行に分けて記述→
3-element Array{Int64,1}:
 1
 2
 3

ベクトルから行列を作る

縦ベクトルを横に並べて行列にできる。固有ベクトルから変換行列を作るときなどわかりやすくなる。

3次のベクトル3つから3×3の行列を作る

v1=[1,0,0]
v2=[1,1,0]
v3=[3,2,2]
P = [v1 v2 v3]
println("P=")
display(P)
println()

並びが行列Aと似ているが、ちょっと違う。Aの固有ベクトルを並べたもの。

P=
3×3 Array{Int64,2}:
 1  1  3
 0  1  2
 0  0  2

inv()関数による逆行列

逆行列には、inv()という関数が使える。ただし、using LinearAlgebra の宣言が必要。

Pの逆行列

using LinearAlgebra

R=inv(P)
println("P^-1=")
display(R)
println()

floatになってしまうが問題はない。

P^-1=
3×3 Array{Float64,2}:
 1.0  -1.0  -0.5
 0.0   1.0  -1.0
 0.0   0.0   0.5

\による割り算ともう一つの逆行列

実は、inv()よりも、こちらの計算を勧めているのを見かけた。理由は計算時間か結果の精度かよくわからないが、今回の範囲では計算結果に違いはない。わかりやすさと、書きやすさでその都度選択する。

Pの逆行列をP-1、単位行列をEとすれば、

PP-1=E

ですから、両辺に左からP-1をかけるというのを、左からPで割り算をすると形式的に考えて、 

P-1=P\E

と書くことができる。

これがjuliaの他では見慣れない割り算の演算子 \ の由来と考えられる。

\を使ってPの逆行列

E = Matrix{Int64}(I, 3, 3)
display(E)
println()
R=P\E
println("P^-1=")
display(R)
println()

単位行列の作り方がちょっと美しくないが、正しく計算できている。

3×3 Array{Int64,2}:
 1  0  0
 0  1  0
 0  0  1
P^-1=
3×3 Array{Float64,2}:
 1.0  -1.0  -0.5
 0.0   1.0  -1.0
 0.0   0.0   0.5

指定された対角行列の生成

diagm()を使えばベクトルから対角行列が作れる

diagm()使用例

G = diagm([5,4,3])
println("diagm")
display(G)
println()

diagm
3×3 Array{Int64,2}:
 5  0  0
 0  4  0
 0  0  3

要素を全部1にすれば単位行列も作れる。

逆行列を使って対角化

元の行列の対角化をしてみる。

matrix004.jl

#!/usr/bin/julia
using LinearAlgebra
A = [
   1 1 2
   0 2 1
   0 0 3
]
print("A   = ")
display(A)
println("")

v1=[1,0,0]
v2=[1,1,0]
v3=[3,2,2]
P = [v1 v2 v3]
println("P=")
display(P)
println()

R=inv(P)
println("P^-1=")
display(R)
println()

D=R*A*P
println("D=")
display(D)
println()

実行結果

A   = 3×3 Array{Int64,2}:
 1  1  2
 0  2  1
 0  0  3
P=
3×3 Array{Int64,2}:
 1  1  3
 0  1  2
 0  0  2
P^-1=
3×3 Array{Float64,2}:
 1.0  -1.0  -0.5
 0.0   1.0  -1.0
 0.0   0.0   0.5
D=
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  2.0  0.0
 0.0  0.0  3.0

対角化だけならDiagonal()関数

Diagonal()を使えば

A = [
   1 1 2
   0 2 1
   0 0 3
]
print("A   = ")
display(A)
println("")

D=Diagonal(A)
println("D=")
display(D)
println()

Diagonal(A)の計算結果は普通の行列でないが、計算は合う。

A   = 3×3 Array{Int64,2}:
 1  1  2
 0  2  1
 0  0  3
D=
3×3 Diagonal{Int64,Array{Int64,1}}:
 1  ⋅  ⋅
 ⋅  2  ⋅
 ⋅  ⋅  3

固有値と固有ベクトル

ここで使った例では固有値はすぐにわかるが、固有値と固有ベクトルを求める関数も用意されている。

この関数の戻り値は2つで、固有値の配列(e)と固有ベクトル(ベクトルは縦)を横に並べた行列(u)である。

eigen()を使う

using LinearAlgebra
A = [
   1 1 2
   0 2 1
   0 0 3
]
print("A   = ")
display(A)
println("")

e,u = eigen(A) 
println("固有値")
display(e)
println("")
println("固有ベクトル")
display(u)
println("")

R = inv(u)

println("P^-1=")
display(R)
println()
D = R*A*u
println("対角化")
display(D)
println()

固有ベクトルは長さが1に規格化されている。対角化まで求めている。

adachi@ebook:julia$ ./matrix006.jl 
A   = 3×3 Array{Int64,2}:
 1  1  2
 0  2  1
 0  0  3
固有値
3-element Array{Float64,1}:
 1.0
 2.0
 3.0
固有ベクトル
3×3 Array{Float64,2}:
 1.0  0.707107  0.727607
 0.0  0.707107  0.485071
 0.0  0.0       0.485071
P^-1=
3×3 Array{Float64,2}:
 1.0  -1.0      -0.5
 0.0   1.41421  -1.41421
 0.0   0.0       2.06155
対角化
3×3 Array{Float64,2}:
 1.0  0.0   0.0
 0.0  2.0  -2.22045e-16
 0.0  0.0   3.0

若干の誤差が見られるが、精度は64ビット(8バイト)で倍精度である。println で表示させると生の値になる。displayは表示では調整して出力しているらしい。

u=
[1.0 0.7071067811865475 0.7276068751089989;
 0.0 0.7071067811865475 0.48507125007266594;
 0.0 0.0                0.48507125007266594]
D=
[1.0   0.0   0.0;
 0.0  2.0 -2.220446049250313e-16;
 0.0   0.0   2.9999999999999996]