はじめに
RustにはNalgebraやndarrayなどの行列を扱うクレートがあります．
その2つ両方を試してみて，体感的にNalgebraのほうが使い心地がよかった(しオシャレな公式ページがあった)ので，そちらの入門記事を書いて見たいと思います．
Nalgebraとndarrayをざっくり比較
Nalgebraのいいところは，まずコンパイル時に行列のサイズを決めておけるところです．
ndarrayは実行時にサイズの決まる動的配列となってしまうので，ややスピードが遅くなると思われます(予想)．
また，サイズチェックなどがコンパイル時であるNalgebraのほうが，実行時エラーを防げます(サイズの異なる行列の積でコンパイルエラーが出ます)．
もちろんNalgebraでも動的配列も扱えます．
それから，網羅的には調べられなかったのですが，演算子のオーバーロードがNalgebraのほうが強い気がします(ndarrayだと+の順序を入れ替えないとコンパイルが通らなかったことがありました．私の調査不足かもしれないですが...)．
ndarrayのいいところはその名の通り多次元配列を扱えるというところです．また，行列の一部を表すスライスなども書きやすいようです．機械学習分野には向いていそうです．ちなみに，Nalgebraにはデフォルトで線形代数モジュールが入っているのに対し，ndarrayにはndarray-linalgという姉妹クレートがあり，線形代数を扱うことができます．
また，ndarrayのほうがNumPyのnpyファイルを扱えるようになるクレートや画像ファイルからの変換クレートなどもあり，何かと便利かもしれないです．
また，よくも悪くもNalgebraの行列にはCopyトレイトが実装されているので，コンパイルエラーは起こりにくいですが，無意識的に，行列のコピーをしてしまう可能性があります．ここは好き嫌いがあると思います(私は無いほうがよかったと感じている)．ndarrayはそんなことはありません．
より詳しいところについては，英語の記事ですが，このページなどで比較できます．
また，ndarrayの使い方に関してはndarray::doc::ndarray_for_numpy_usersが非常にわかりやすかったです．また，かゆいところに手が届く用に書いたつもりの拙著[Rust] ndarray入門も是非ご覧ください．
別のcgmathというクレートも3〜4次元程度しか扱わなければ，次元固定なので便利そうですね．クォータニオンも扱えるらしいですし．
コンパイル時サイズ決定の方法
サイズはnalgebra::U{n}({n}に0,1,...,127が入る)型で与えられます．これを型パラメータにいれることで，コンパイル時にサイズを決めます．またサイズを動的に決めたいときはnalgebra::Dynamicとすれば大丈夫です．
型がないサイズにはコンパイルできないことには注意が必要です．ただクレートtypenumを使うともっと大きな次元で利用できるそうです(ざっと見たところ1024までと10の累乗は使える)．詳しくはこちらに説明があります．このクレートでは，型同士の数字の四則演算ができるので，うまく使えば行列のサイズの型をconstなusizeのように扱えます．
行列の型
以下の説明では全て
extern crate nalgebra as na;
としておきます(公式もそうしているので)．
一般の行列については
na::MatrixMN<N, R, C>
という型を用います. Nは各成分のデータの型で，R，Cには上記のna::U{n}またはna::Dynamicが入り，Rが行数，Cが列数を表します．
各成分がf64の行列を作りたいときにはNにf64を入れておけば大丈夫です．
また等価な型も多くあり
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次元行ベクトル
などと定義されています(型パラメータの記号は説明のために適宜変えた)．
また，動的配列については
pub type DMatrix<N> = MatrixN<N, Dynamic>; // 行も列もDynamicなので正方行列とは限らない
動的ベクトルについては
pub type DVector<N> = Matrix<N, Dynamic, U1, VecStorage<N, Dynamic, U1>>; // 詳しい説明は省略
と定義されています．
初期化
初期化は
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(); // 零行列
などとします．ランダム初期化は
let m = na::Matrix3::new_random();
とします．ランダムのシード設定はよくわからなかったので教えてもらえると助かります．
動的行列については，公式サイトがよくまとまっていたのでコピペしておきます．
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);
その他の方法はここにあります．
主な演算
以下vはベクトル，m, m0, m1, m2は行列とします．
下記のようにします．静的な行列にしておけば，コンパイル時にサイズチェックを行います．
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; // スカラーとの計算
成分ごとの積・商については以下のようにします．
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と等価です．
各成分へのアクセス
let a = v[1];
let b = m[(1, 2)];
v[0] = 4;
m[(0, 1)] = -1;
などとします．get_uncheckedなどの，高速だがunsafeなアクセス方法もあります．
スライスの取得
スライスの取得方法には静的なものと動的なものがあります(できるだけ静的なスライスを用いると良いと思います)．ここに視覚的な説明があるので，併せて参考にしてください．
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という代入
例えば，
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));
となります．
サイズの変更
例えば
m.fixed_resize::<U2, U3>(57) // m.reshape(2, 3) ただし前の行列が2*3=6成分未満のときは足りない分を57で補った行列を返す
m.resize(2, 3, 57) // 上の動的ver
とできます．57という数字に意味はありません1．
その他
公式のクイックリファレンスから，よく使いそうなメソッドをいくつか紹介していきます．
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度だけしておくようにしましょう．
終わりに
余裕があったら，書き足していこうと思います！いいRustライフを！