数値計算においては、ある行列$A$とベクトル$v$の積$Av$を使った演算が色々ある。
例えば、連立方程式$Ax=b$の解を求めるための共役勾配法(CG法)における逐次演算では行列ベクトル積を何回も演算させる。
しかし、行列$A$があまりにも大きかった場合、メモリに$A$が入らずに計算ができない!というような場合がある。一つの方法としては、$A$が疎行列であれば、行列要素のうち0ではない部分だけを格納することで、メモリを節約することができる。Juliaではその場合にはSparseArrays.jlというパッケージがある。
一方、非ゼロ要素が沢山ある場合でも、行列自体をメモリに保持するのではなく、あるベクトル$v$に対して$Av$のみを返すようにすれば、計算が可能となる。
そのようなケースに対応するパッケージが、LinearMaps.jlというパッケージである。
https://github.com/Jutho/LinearMaps.jl
これを使えば、線形な演算子$A$が与えられているときに$Av$を計算するというようなより抽象的な演算も行う事ができる。
バージョン
Julia 1.0
線形演算の定義
まず、線形演算を定義してみよう。ベクトル$x$が与えられた時に、ベクトル$y$が返ってくる関数$y=f(x)$を作ってみる。
これは、
function set_diff(v)
function calc_diff!(y::AbstractVector, x::AbstractVector)
n = length(x)
length(y) == n || throw(DimensionMismatch())
for i=1:n
dx = 1
jp = i+dx
jp += ifelse(jp > n,-n,0) #+1方向
dx = -1
jm = i+dx
jm += ifelse(jm < 1,n,0) #-1方向
y[i] = v*(x[jp]+x[jm])
end
return y
end
(y,x) -> calc_diff!(y,x)
end
として、
A = set_diff(-1.0)
とすればよい。これで、A(y,x)というような関数を作ることができた。
この関数を呼ぶと、
y = D x \\
y_i = \sum_{i=1}^N (D_{i,i+1} x_{i+1}+D_{i,i-1} x_{i-1})
というような行列Dをxに演算した結果のベクトルyが出てくることになる。
functionの中にfunctionを入れることで、パラメータvが変化したDを色々作れるようになっている。
LinearMaps
LinearMaps.jlを使ってみよう。
行列として、
using LinearMaps
n = 100
D = LinearMap(A, n; ismutating=true,issymmetric=true)
とする。ここで、issymmetric=trueは行列が対称行列であるということを指示しており、ismutating=trueはA(x)ではなくA(y,x)という関数であることを意味している。これをtrueにしておくと、古いyとxを使って新しいyを作る演算が定義できるようになる。nは行列のサイズである。
この行列の固有値を求めたければ、Arpack.jlを使えばよくて、
import Arpack
ene,v = Arpack.eigs(D; nev=20, which=:SR)
println(ene)
とすれば、20個固有値を得ることができる。なお、issymmetric=trueをつけないとArpackが対称行列用のルーチンを使わなくなるので、固有値が複素数になっても大丈夫なように演算することとなる。あらかじめ対称行列だとわかっているのであれば、issymmetric=trueをつけた方が計算が軽い。
CG法で連立方程式を解きたい場合には、IterativeSolvers.jlを使えばよい。ArpackとIterativeSolversはLinearMapsに対応している。
LinearMapsを使えば、4次元時空中のディラック演算子を作用させるなど、行列で書くのはしんどい演算を簡単に定義することができて、パッケージを使って様々なことができるようになる。
また、GPUを使った計算や並列計算なども見通しよく記述できるはずである。