|
setup diary | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||
プログラム言語で数値計算をしようとするときに，行列を扱う必要がある場合がある．そのやり方が言語毎に特徴があるので，比較してみた．例として，行列を定義して，固有値を固有ベクトルを求めて，対角化を計算する場合について，やり方を見てみよう．
まず，私が数値計算をするときによく使っているRでは，インストールするだけで，行列を使うことができる．逆行列を求めるコマンドがsolveである点と，行列の掛け算が%*%であることに注意すれば，簡単に使うことができる．
a<-t(matrix(1:4,2,2)) r<-eigen(a) r$values r$vectors solve(r$vectors)%*%a%*%r$vectors
次に，近年人気が出てきたjuliaであるが，これもインストールした時点で行列自体は使えて，さらにLinearAlgebraを組み込むことによって，固有値などの計算をすることができる．
a=[1 2; 3 4] using LinearAlgebra r=eigen(a) r.values r.vectors inv(r.vectors)*a*r.vectors
数値計算の分野でもシェアを広げているpythonでは，numpyを使う必要がある．pipなどでこれをインストールできる．私の場合は，debianでaptitudeからpython3-numpyを入れた．numpyをimportして，linalgを使うと計算ができる．行列の掛け算がnp.dot()で書かないといけないのが面倒だったのだが，新しいversionでは@を使えるようになったようだ．
import numpy as np a=np.array([[1,2],[3,4]]) r=np.linalg.eig(a) r[0] r[1] #np.dot(np.dot(np.linalg.inv(r[1]),a),r[1]) np.linalg.inv(r[1])@a@r[1]
最後にrubyであるが，数値計算は比較的苦手である．行列計算はnarrayを使えば良いのだが，固有値の計算などは困難である．そこで，numo-linalgというものを使ってみた．このインストールは少し手間取ったが，以下のようにしたらdebianにインストールすることができた．
sudo apt install liblapacke libopenblas0 sudo apt install git ruby gcc ruby-dev rake make sudo gem install specific_install sudo gem specific_install https://github.com/ruby-numo/numo-narray.git sudo gem install numo-linalg -- --with-backend=openblas --with-openblas-lib=/usr/lib/x86_64-linux-gnu
すると，numo/linalgをrequireしたら計算できるのだが，いちいちNumoと書くのが面倒なので，それをincludeした．narrayの場合には，NMatrixを定義すると，行列の積は*で書けるのだが，numoの場合には，numo-narrayのdotを使って書かないといけないのが面倒に感じた．
require "numo/linalg" include Numo a=NArray[[1,2],[3,4]] r=Linalg.eig(a) r[0] r[2] Linalg.inv(r[2]).dot(a).dot(r[2])
これらの4つの言語について，インストールのしやすさと，使いやすさを主観でまとめると，下の表のようになった．Rとjuliaはそれなりに使いやすいかな．pythonは行列の積が書きやすくなったことにより，昔よりは書きやすくなった．rubyはインストールが面倒だし，methodをうまく定義して，もっと使いやすくできるポテンシャルを秘めているが，まだまだのように感じる．ニーズが無いのかな．
| install | use | |
|---|---|---|
| R | ◎ | ◎ |
| julia | ◎ | ◎ |
| python | ○ | △→○ |
| ruby | △ | ○ |
数値計算をする場合には，Rを中心にやっていくことになりそうだ．juliaも使うようにすると，Rが苦手な処理を書きやすくなるかも知れない．
| R | julia | ruby | python3 | |
| 余 | 5%%2 | 5.0%2.0 | 5.0%2.0 | 5.0%2.0 |
| 商 | 5%/%2 | 5.0÷2.0 | 5/2 or 5.0.div(2.0) | 5.0//2.0 |
| 分数 | 5//2 | 5r/2 | from fractions import Fraction;Fraction(5,2) | |
| xor | bitwXor(3, 5) | 3⊻5 | 3^5 | 3^5 |
| 累乗 | 2^0.5 | 2^0.5 | 2**0.5 | 2**0.5 |
| 平方根 | sqrt(2) | √2 | Math.sqrt(2) | import math;math.sqrt(2) |
| 複素数 | (1+2i)^0.5 | (1+2im)^0.5 | (1+2i)**0.5 | (1+2j)**0.5 |
| 複素数 | sqrt(1+2i) | √(1+2im) | Math.sqrt(1+2i) | import numpy;numpy.sqrt(1+2j) |
| 実部 | Re(1+2i) | real(1+2im) | (1+2i).real | (1+2j).real |
| 虚部 | Im(1+2i) | imag(1+2im) | (1+2i).imag | (1+2j).imag |
| 共役 | Conj(1+2i) | conj(1+2im) | (1+2i).conj | (1+2j).conjugate() |
| 偏角 | Arg(1+2i) | angle(1+2im) | (1+2i).arg | import cmath;cmath.phase(1+2j) |
| 絶対値 | abs(-1) | abs(-1) | -1.abs | abs(-1) |
| 円周率 | pi | π | Math::PI | import math;math.pi |
| Euler数 | exp(1) | ℯ | Math::E | import math;math.e |
| R | julia | ruby | python3 | |
| 繰返 | paste(rep("abc",3),collapse="") | "abc"^3 | "abc"*3 | "abc"*3 |
| 連結 | paste("abc","def",sep="") | "abc"*"def" | "abc"+"def" | "abc"+"def" |
| 結合 | paste(1:5,collapse=", ") | join(1:5,", ") | [*1..5]*", " | ", ".join(map(str,range(1,6))) |
| 取出 | substr("πℯ^2",2,2) | collect("πℯ^2")[2] | "πℯ^2"[1] | "πℯ^2"[1] |
| 部分 | substr("πℯ^2",2,3) | String(collect("πℯ^2")[2:3]) | "πℯ^2"[1..2] "πℯ^2"[1,2] | "πℯ^2"[1:3] |
| 文字数 | nchar("πℯ^2") | length("πℯ^2") | "πℯ^2".length "πℯ^2".size | len("πℯ^2") |
| 大文字 | toupper("This") | uppercase("This") | "This".upcase | "This".upper() |
| 小文字 | tolower("This") | lowercase("This") | "This".downcase | "This".lower() |
| 入れ替え | chartr("A-Za-z","a-zA-Z","This") | "This".swapcase | "This".swapcase() | |
| 埋込 | s="world" "Hello $(s)!" | s="world" "Hello #{s}!" | s="world" f"Hello {s}!" "Hello {}!".format(s) | |
| 書式 | sprintf("%d %s",3,"times") | using Printf;@sprintf("%d %s",3,"times") | sprintf("%d %s",3,"times") "%d %s"%[3,"times"] | "%d %s"%(3,"times") |
| 文字列 | as.character(123) | string(123) | 123.to_s | str(123) |
| 整数 | as.integer("1") | parse(Int,"1") | "1".to_i | int("1") |
| 浮動小数点 | as.numeric("1") | parse(Float64,"1") | "1".to_f | float("1") |
| 複素数 | as.complex("1+2i") | parse(Complex{Int},"1+2im") | "1+2i".to_c | complex("1+2j") |
| コード変換 | utf8ToInt("A") | Int('A') | "A".ord | ord("A") |
| 文字変換 | intToUtf8(960) | Char(960) | 960.chr(Encoding::UTF_8) | chr(960) |
| 検索 | regexpr("[\\d\\.]+","pi=3.14",perl=TRUE) | match(r"[\d\.]+","pi=3.14") | "pi=3.14"=~/[\d\.]+/ | import re;re.search("[\d\.]+","pi=3.14") |
| 置換 | sub("(\\d)\\.",".\\1","pi=3.14") | replace("pi=3.14", r"(\d)\." => s".\1") | "pi=3.14".sub(/(\d)\./,'.\1') | import re;re.sub("(\d)\.",".\\1","pi=3.14") |
| 分割 | strsplit("1,2, 3",",\\s*")[[1]] | split("1,2, 3", r",\s*") | "1,2, 3".split(/,\s*/) | import re;re.split(",\s*","1,2, 3") |
| R | julia | ruby | python3 | |
| if | a<-3 if(a==1){cat(1) }else if(a==2){cat(2) }else{cat(3)} | a=3 if a==1 print(1) elseif a==2 print(2) else print(3) end | a=3 if a==1 then p 1 elsif a==2 then p 2 else p 3 end | a=3 if a==1:print(1) elif a==2:print(2) else:print(3) |
| 三項演算子 | a<-2 cat(ifelse(a>0,"+","-")) | a=2 print(a>0 ? "+" : "-") | a=2 print(a>0?"+":"-") | a=2 print('+' if a>0 else '-') |
| switch | i<-"a" cat(switch(i,a=1,b=2)) | i="a" p case i when "a" then 1 when "b" then 2 end | ||
| for | for(i in 1:10){cat(i);cat("\n")} | for i=1:10 println(i) end | for i in 1..10 do p i end | for i in range(1,11):print(i) |
| while | i<-0 while(i<10){cat(i<-i+1);cat("\n")} | let i=0 while i<10 println(i+=1) end end | i=0 while i<10 do p i+=1 end | i=0 while i<10:i+=1;print(i) |
| 無限ループ | repeat{cat(1);cat("\n")} | while true println(1) end | loop{p 1} | while True:print(1) |
| ループを抜ける | break | break | break | break |
| 次のループ | next | continue | continue | next |
| 例外処理 | tryCatch(cat("1"/0), error=function(e)cat(-1)) | try print("1"/0) catch print(-1) end | begin p "1"/0 rescue p -1 end | try:print("1"/0) except:print(-1) |
| 終了処理 | tryCatch({cat(1);cat("\n")}, finally=cat(2)) | try println(1) finally println(2) end | begin p 1 ensure p 2 end | try:print(1) finally:print(2) |
| 複文 | a<-{cat(1);2} | a=(print(1);2) | a=(p 1;2) |