Rustによる数値計算: 線形代数編 #Rust

Rust Advent Calendar 12/4の記事です

以前にRustで1次元のNewton法を実装する記事を書きましたが、実際の科学技術計算では多次元のベクトルの扱いが必要になります。
PythonではNumPyで提供されているnumpy.ndarrayが中心的な役割を果たしますが、
今回はRustによる多次元配列の実装であるrust-ndarrayを使用します。
しかしこのcrateにはベクトルの演算と行列xベクトルの演算までしか実装されておらずQR分解や特異値分解、
あるいは高速Fourier変換のような演算は実装されていません。
このままだと自分でアルゴリズムを書くのに支障があるため、基本的な部分を実装してきました。

rust-ndarray

rust-ndarrayはnumpyのものと同じで一つの連続したメモリとstridesと呼ばれる各インデックス毎のとび幅からなっています。
この辺はnumpyのstrideのドキュメントを参照してください。
このシンプルな構成が好みでrust-ndarrayを選択しました。
簡単にrust-ndarrayの使い方を確認しておきましょう。

extern crate ndarray;
use ndarray::prelude::*;

let a = Array::<f64, _>::zeros(3)
// [0.0, 0.0, 0.0]
let b = Array::range(0., 9., 1.).into_shape((3, 3)).unwrap();
// [[0, 1, 2]
//  [3, 4, 5]
//  [6, 7, 8]]

乱数を生成するときはndarray-randを使います：

extern crate rand;
extern crate ndarray_rand;
use rand::distributions::*;
use ndarray_rand::RandomExt;

let r_dist = Range::new(0., 1.);
let a = Array::<f64, _>::random((10, 10), r_dist); // 10x10の乱数でうまった行列

ndarray-linalg

線形代数演算をrust-ndarray上に実装したものがndarray-linalgです。
OpenBLAS(旧:GotoBLAS)をrustにバインドしたlapack crateをrust-ndarrayに用にバインドしたものです。
類似のプロジェクトとしてlinxalがあります。
ndarray-linalgではQR分解・SVD等の演算をtraitとして実装してあります。
現状まだ実数の行列にしか対応していません(v0.2で対応の予定)。以下の演算が実装されています：

trait Matrix
- singular-value decomposition
- LU decomposition
- QR decomposition
- 作用素norm (L1/L^inf)
- Frobeiuns norm
trait SquareMatrix
- inverse of matrix
- trace of matrix
- [WIP] eigenvalue
trait TriangularMatrix
- solve linear problem with upper triangular matrix
- solve linear problem with lower triangular matrix
trait HermiteMatrix
- eigenvalue analysis
- symmetric square root
- Cholesky factorization

すべてtraitで実装されているため、rust-ndarrayのinstanceに対して勝手に演算が定義されています：

extern crate ndarray;
extern crate ndarray_linalg;

use ndarray::prelude::*;
use ndarray_linalg::prelude::*;

fn main() {
    let a = arr2(&[[3.0, 1.0, 1.0], [1.0, 3.0, 1.0], [1.0, 1.0, 3.0]]);
    let (e, vecs) = a.clone().eigh().unwrap(); // 固有値計算(エルミート)
    println!("eigenvalues = \n{:?}", e);
    println!("V = \n{:?}", vecs);
    let av = a.dot(&vecs);
    println!("AV = \n{:?}", av);
}

各種演算は失敗するかもしれないのでResult<_, LinalgError>をほとんどの演算が返します。

（おまけ）crateを公開する

cargo package
cargo publish

によって自作のcrateをcrates.ioに公開できます(事前にcrates.ioにログインしておく必要があります)。
ここに公開させるとしばらくしてからdocs.rsに反映されます。
ここではドキュメントがcargo docで生成され、自動的にホストしてくれます。

最後に

現在までに他にもrust-ndarrayの拡張として以下のcrateを作ってるのでissue/PR募集中です。

ndarray-odeint : ODEのソルバ集(の予定)
ndarray-fftw : FFTWのバインド
ndarray-numtest : テスト用コードの集積