数値計算においては、ある行列$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を使った計算や並列計算なども見通しよく記述できるはずである。