[Rust] Nalgebra入門(ndarrayとの比較付き)

# はじめに

Rustには[Nalgebra](https://www.nalgebra.org/)や[ndarray](https://github.com/rust-ndarray/ndarray)などの行列を扱うクレートがあります．
その2つ両方を試してみて，体感的にNalgebraのほうが使い心地がよかった(しオシャレな公式ページがあった)ので，そちらの入門記事を書いて見たいと思います．

# Nalgebraとndarrayをざっくり比較

Nalgebraのいいところは，まず**コンパイル時に行列のサイズを決めておける**ところです．
ndarrayは実行時にサイズの決まる動的配列となってしまうので，ややスピードが遅くなると思われます(予想)．
また，サイズチェックなどがコンパイル時であるNalgebraのほうが，実行時エラーを防げます(サイズの異なる行列の積でコンパイルエラーが出ます)．
もちろんNalgebraでも動的配列も扱えます．

それから，網羅的には調べられなかったのですが，演算子のオーバーロードがNalgebraのほうが強い気がします(ndarrayだと+の順序を入れ替えないとコンパイルが通らなかったことがありました．私の調査不足かもしれないですが...)．

ndarrayのいいところはその名の通り**多次元配列を扱える**というところです．また，行列の一部を表すスライスなども書きやすいようです．機械学習分野には向いていそうです．ちなみに，Nalgebraにはデフォルトで線形代数モジュールが入っているのに対し，ndarrayにはndarray-linalgという姉妹クレートがあり，線形代数を扱うことができます．

また，ndarrayのほうが[NumPyのnpyファイルを扱えるようになるクレート](https://docs.rs/ndarray-npy/0.7.1/ndarray_npy/)や[画像ファイルからの変換クレート](https://docs.rs/ndarray-image/0.3.0/ndarray_image/)などもあり，何かと便利かもしれないです．

また，よくも悪くもNalgebraの行列にはCopyトレイトが実装されているので，コンパイルエラーは起こりにくいですが，無意識的に，行列のコピーをしてしまう可能性があります．ここは好き嫌いがあると思います(私は無いほうがよかったと感じている)．ndarrayはそんなことはありません．

より詳しいところについては，英語の記事ですが，[このページ](https://www.reddit.com/r/rust/comments/btn1cz/ndarray_vs_nalgebra/)などで比較できます．

また，ndarrayの使い方に関しては[ndarray::doc::ndarray_for_numpy_users](https://docs.rs/ndarray/0.14.0/ndarray/doc/ndarray_for_numpy_users/index.html)が非常にわかりやすかったです．また，かゆいところに手が届く用に書いたつもりの拙著[[Rust] ndarray入門](https://qiita.com/ttttkkkkk31525/items/b75574bbff52e6cd213a)も是非ご覧ください．

別の[cgmath](https://docs.rs/crate/cgmath/0.18.0)というクレートも3〜4次元程度しか扱わなければ，次元固定なので便利そうですね．クォータニオンも扱えるらしいですし．


# コンパイル時サイズ決定の方法
サイズは`nalgebra::U{n}`({n}に0,1,...,127が入る)型で与えられます．これを型パラメータにいれることで，コンパイル時にサイズを決めます．またサイズを動的に決めたいときは`nalgebra::Dynamic`とすれば大丈夫です．

型がないサイズにはコンパイルできないことには注意が必要です．ただクレート[typenum](https://docs.rs/crate/typenum/)を使うともっと大きな次元で利用できるそうです(ざっと見たところ1024までと10の累乗は使える)．詳しくは[こちら](https://nalgebra.org/vectors_and_matrices/)に説明があります．このクレートでは，型同士の数字の四則演算ができるので，うまく使えば行列のサイズの型を`const`な`usize`のように扱えます．

# 行列の型

以下の説明では全て

```rust
extern crate nalgebra as na;
```

としておきます(公式もそうしているので)．
一般の行列については

```rust

na::MatrixMN<N, R, C>
```

という型を用います. Nは各成分のデータの型で，R，Cには上記の`na::U{n}`または`na::Dynamic`が入り，Rが行数，Cが列数を表します．
各成分が`f64`の行列を作りたいときにはNに`f64`を入れておけば大丈夫です．

また等価な型も多くあり

```rust
use na::{U1, U2, U3, Dynamic};

pub type MatrixN<N, D> = MatrixMN<N, D, D>;     // 正方行列
pub type VectorN<N, R> = MatrixMN<N, R, U1>;    // 列ベクトル
pub type RowVectorN<N, C> = MatrixMN<N, U1, C>; // 行ベクトル

pub type Matrix2x3<N> = MatrixMN<N, U2, U3>;    // 2x3行列
pub type Matrix2<N> = MatrixN<N, U2>;           // 2次正方行列
pub type Vector3<N> = VectorN<N, U3>;           // 3次元列ベクトル
pub type RowVector3<N> = RowVectorN<N, U3>;     // 3次元行ベクトル
```
などと定義されています(型パラメータの記号は説明のために適宜変えた)．

また，動的配列については

```rust
pub type DMatrix<N> = MatrixN<N, Dynamic>; // 行も列もDynamicなので正方行列とは限らない
```
動的ベクトルについては

```rust
pub type DVector<N> = Matrix<N, Dynamic, U1, VecStorage<N, Dynamic, U1>>; // 詳しい説明は省略
```
と定義されています．

# 初期化
初期化は

```rust
let v = na::Vector3::new(1, 2, 3);

let m = na::Matrix3x4::new(11, 12, 13, 14,
                           21, 22, 23, 24,
                           31, 32, 33, 34);
let u = na::Matrix3x4::zeros(); // 零行列
```
などとします．ランダム初期化は

```rust
let m = na::Matrix3::new_random();
```
とします．ランダムのシード設定はよくわからなかったので教えてもらえると助かります．

動的行列については，公式サイトがよくまとまっていたのでコピペしておきます．

```rust
let dm = DMatrix::from_row_slice(4, 3, &[
    1.0, 0.0, 0.0,
    0.0, 1.0, 0.0,
    0.0, 0.0, 1.0,
    0.0, 0.0, 0.0
]);
let dm1 = DMatrix::from_vec(4, 3, vec![1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]);
let dm2 = DMatrix::from_diagonal_element(4, 3, 1.0);
let dm3 = DMatrix::identity(4, 3);
let dm4 = DMatrix::from_fn(4, 3, |r, c| if r == c { 1.0 } else { 0.0 });
let dm5 = DMatrix::from_iterator(4, 3, [
    // Components listed column-by-column.
    1.0, 0.0, 0.0, 0.0,
    0.0, 1.0, 0.0, 0.0,
    0.0, 0.0, 1.0, 0.0
].iter().cloned());

assert_eq!(dm, dm1); assert_eq!(dm, dm2);
assert_eq!(dm, dm3); assert_eq!(dm, dm4);
assert_eq!(dm, dm5);
```

その他の方法は[ここ](https://www.nalgebra.org/vectors_and_matrices/#matrix-construction)にあります．

# 主な演算
以下`v`はベクトル，`m`, `m0`, `m1`, `m2`は行列とします．

下記のようにします．静的な行列にしておけば，コンパイル時にサイズチェックを行います．

```rust
let _ =  m0 +  m1; // ok
let _ =  m0 + &m1; // ok
let _ = &m0 +  m1; // ok
let _ = &m0 + &m1; // ok
// その他の演算子も上と同様に可能
let _ = &m0 - &m1; 
let _ = &m0 * &m1; // 行列積(成分ごとの積ではない)

let _ = m / 2 * 3 - 4 * &m1; // スカラーとの計算
```
成分ごとの積・商については以下のようにします．

```rust
m0.component_mul(&m1)        // 成分ごとの積(アダマール積)
m0.component_mul_assign(&m1) // 破壊的に変更するアダマール積
m0.component_div(&m1)        // 成分ごとの商
m0.component_div_assign(&m1) // 破壊的に変更する成分ごとの商
```

また，`map`，`apply`(破壊的)，`zip_map`，`zip_apply`等のメソッドで各成分に関数を適用することも可能です．NumPyのような値の比較も`m.map(|a| a > 0)`や`m0.zip_map(&m1, |x, y| x > y)`などとします．

スカラーを足すときは`m.add_scalar(1)`や`m.add_scalar_mut(1)`(破壊的)などとします．`sub_scalar`はないので，-1を足すというようにします．

そして，ベクトルの内積には`dot`というメソッドがあり，
`dot`は`self.transpose() * rhs`と等価です．

複素行列用の
`dot_c`は`self.adjoint() * rhs`と等価です．

# 各成分へのアクセス

```rust
let a = v[1];
let b = m[(1, 2)];
v[0] = 4;
m[(0, 1)] = -1;
```
などとします．`get_unchecked`などの，高速だが`unsafe`なアクセス方法もあります．


# スライスの取得

スライスの取得方法には静的なものと動的なものがあります(できるだけ静的なスライスを用いると良いと思います)．[ここ](https://www.nalgebra.org/vectors_and_matrices/#matrix-slicing)に視覚的な説明があるので，併せて参考にしてください．

```rust
m.row(1) // mの1行目 NumPy風に書けばm[1]
m.column(0) // mの0列目 NumPy風に書けばm[:, 0]

// 静的なスライス
m.fixed_rows::<U2>(i)                    // NumPy風に書けばm[i:i+2]の2行
m.fixed_columns::<U2>(i)                    // NumPy風に書けばm[:, i:i+2]の2列
m.fixed_slice::<U2, U3>(i, j)            // NumPy風に書けばm[i:i+2, j:j+3]
m.fixed_rows_mut::<U2>(i).copy_from(&m2) // m[i:i+2] = m2という代入

// 動的なスライス
m.rows(i, 2)                             // NumPy風に書けばm[i:i+2]の2行
m.columns(i, 2)                             // NumPy風に書けばm[:, i:i+2]の2列
m.slice((i, j), (2, 3))                  // NumPy風に書けばm[i:i+2, j:j+3]
m.rows_mut(i, 2).copy_from(&m2)          // m[i:i+2] = m2という代入
```

例えば，

```rust
let mut x = na::Vector5::new(0, 1, 2, 3, 4);
let v = na::Vector3::new(0, -1, -2);

x.fixed_rows_mut::<na::U3>(2).copy_from(&v);
assert_eq!(x, na::Vector5::new(0, 1, 0, -1, -2));
```
となります．

# サイズの変更
例えば

```rust
m.fixed_resize::<U2, U3>(57) // m.reshape(2, 3) ただし前の行列が2*3=6成分未満のときは足りない分を57で補った行列を返す
m.resize(2, 3, 57) // 上の動的ver
```
とできます．57という数字に意味はありません[^1]．
[^1]: 57は素数ではない．

# その他
[公式のクイックリファレンス](https://www.nalgebra.org/quick_reference/)から，よく使いそうなメソッドをいくつか紹介していきます．

L2ノルム(ベクトルの長さ)は`v.norm()`
その2乗は`v.norm_squared()`で求められます(2乗のほうが平方根を求めない分速い).

転置行列は`m.transpose()`で得られます．また，In-placeで得たいときには`m.transpose_mut()`で`m`を変更できます(正方行列か動的な行列のみ)．`m.transpose_to(output)`でoutputに転置行列の結果が入ります．
`transpose`を`adjoint`に読み替えることで，複素行列の共役転置(エルミート共役)も可能です．

逆行列は`m.try_inverse()`で`Option<_>`で得られます．
また，`m`を破壊的に変更して逆行列を求める`m.try_inverse_mut()`もあります(正則でない，つまり逆行列が存在しないときは`false`を返す)．

なお，大きな行列で$A^{-1}b$ ($A$は行列で$b$はベクトル)を求めるときは，逆行列を求めるのではなく，$Ax = b$という方程式を($A$を一度LU分解して)解くようにしましょう．`A.lu().solve(&b)`[^2]や`A.lu().solve_mut(&mut b)`(破壊的)で解くことができます．そのほうが計算量的にかなり有利です．計算量として$7n^3/6$くらい得をします($n$は行列のサイズ)．
繰り返し同じ$A$で計算を行うときには，LU分解を1度だけしておくようにしましょう．
[^2]: 説明の都合上，行列の変数名を大文字にしましたが，変数名は小文字から始めるのをおすすめします．

# 終わりに
余裕があったら，書き足していこうと思います！いいRustライフを！