Rustで線形代数 (ndarray-linalgの使い方)

ndarray-linalg 0.9.0のリリース記念に簡単な使い方をまとめます

GitHub: https://github.com/termoshtt/ndarray-linalg
docs.rs: https://docs.rs/ndarray-linalg/0.9.0/ndarray_linalg/

Cargo.toml

基本的にはrust-ndarrayと組み合わせて使います

Cargo.toml

[dependencies]
ndarray = "0.11"
ndarray-linalg = { version = "0.9", features = ["openblas"] }

featureにはLAPACKのバックエンドを入れます。一つだけ指定しないとエラーになります。現在対応してるのは以下の通り：

• OpenBLAS (openblas)
• Netlib (netlib)
• Accelerate (accelerate)
• Intel MKL (intel-mkl)

ndarray-linalgの使い方

lib.rs

#[macro_use]
extern crate ndarray;
#[macro_use]
extern crate ndarray_linalg;

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

以下名前空間を省略してます

複素数 c32, c64

type c32 = Complex<f32>;
type c64 = Complex<f64>;

Complexはnum-complexのものです

Generate

let a: Array3<f64> = ndarray_linalg::random((3, 2, 3));
let a: Array1<c64> = ndarray_linalg::random(10);
let h: Array2<f64> = ndarray_linalg::random_hermite(5);
let h: Array2<f64> = ndarray_linalg::random_hpd(5);

乱数で多次元配列を初期化します。主に初期化・テスト用です。random_hpdはHermite Positive Definiteの略です。
戻り値の型で判定するので推定できない場合はアノテーション付ける必要があります。

Norm

let a = ndarray_linalg::random((3, 3, 4));
a.norm_l1();  // L1-norm
a.norm_l2();  // L2-norm
a.norm_max(); // 最大値ノルム

全要素の絶対和(L1)二乗和(L2)最大値のいずれかでノルムを取ります。任意の次元で使えます。
行列の場合は作用素ノルムも使えます：

let a = ndarray_linalg::random((3, 3));
a.opnorm_one();
a.opnorm_inf();
a.opnorm_fro();

QR分解

let a: Array2<f64> = ndarray_linalg::random((5,4));
let (q, r) = a.qr().unwrap();

これはaを非破壊で結果を新しい配列(Array2<f64>)に入れます。正方行列でなくても使えます。

let a: Array2<f64> = ndarray_linalg::random((5,4));
let (q, r) = a.qr_into().unwrap();

qr_into()にするとaをmoveします。配列の確保が減ります。

let a: Array2<f64> = ndarray_linalg::random((10, 10));
let (q, r) = a.qr_square().unwrap();

正方行列を仮定してる場合です。ちょっとアルゴリズムが簡単になります。aは消費しません。qr_square_into()にすると消費します。消費する/しないは引数がself(消費する)か&self(消費しない)を見て判断します。

let mut a: Array2<f64> = ndarray_linalg::random((10, 10));
let (q, r) = a.qr_square_inplace().unwrap();

最後にinplaceで変換するタイプです。rは新しく生成します。qにはaの&mutが入っています。

このようにa.method()の形で基本的に呼び出せるようにしてあって、元の行列を消費する・しないをメソッド名で決定します。

SVD

let a: Array2<f64> = ndarray_linalg::random((10, 10));
let (u, s, vt) = a.svd(True, True).unwrap();

基本的にはQR分解と一緒ですが、特異値分解 $A = U\Sigma V^T$で$U$と$V^T$を計算しない選択肢があるのでそれをboolであげます。なのでu, vtはOption<Array2<f64>>とかになります。これもsvd_into(...), svd_inplace(...)があります。

Solve(h)

let a: Array2<f64> = ndarray_linalg::random((3,3));
let b: Array1<f64> = ndarray_linalg::random(3);
let x = a.solve_into(b).unwrap();

線形方程式 $Ax = b$を解きます。LU分解で解きますが複数$b$がある場合$A$の分解を記憶しておく方がいいの、その場合は以下のように明示的に分解します：

let a: Array2<f64> = random((3, 3));
let f = a.factorize_into().unwrap(); // LU factorize A (A is consumed)
for _ in 0..10 {
    let b: Array1<f64> = random(3);
    let x = f.solve_into(b).unwrap(); // Solve A * x = b using factorized L, U
}

Hermite（あるいは実対称）行列の場合は同じインターフェースがsolvehモジュールに用意されています。
https://docs.rs/ndarray-linalg/0.9.0/ndarray_linalg/solveh/index.html

最後に

他にもいくつか機能があります！（飽きた(´・ω・｀)
続きはWEBで！ https://docs.rs/ndarray-linalg/0.9.0/ndarray_linalg/index.html