Rust と nalgebra で MLP を実装した話

## はじめに

はじめまして，CADDi でエンジニアをしている和田です．

この記事は [CADDi Advent Calendar](https://qiita.com/advent-calendar/2020/caddi) 14 日目の記事です．昨日は，[@ryokotmng](https://twitter.com/ryokotmng)さんによる[Rustでパラメーター化テスト](https://caddi.tech/archives/1849)でした！

CADDi では図面画像を解析するアルゴリズムの開発を Rust を用いて行っています．この際， 各種線形代数演算を実現するために[nalgebra](https://nalgebra.org/)を利用しています．この記事では，学習の一貫として nalgebra を利用して MLP(Multi Layer Perceptron) を実装し，MNIST に対して分類を行った例を紹介させていただきます．

## nalgebra について
nalgebra は Rust で記述された線形代数ライブラリです．現在でも精力的に開発が続けれられています．nalgebra は非常に多機能なライブラリです．従って，MLP の実装に用いた機能に絞って以下に紹介します．

### 行列の初期化
nalgebra には用途に応じた複数の行列を表現する型とエイリアスを提供します．今回，任意の次元が指定可能な `MatrixMN` を使用します．これは，各レイヤの入出力の次元を柔軟に設定可能とするためです．ここで，`N` は要素の型を，`R` は行数を，`C` は列数をそれぞれ指定します．行数および列数の指定は nalgebra が提供する型レベル整数を使用します．

初期化についても同様に，用途に応じた複数の方法を提供します．従って，今回の実装で用いらる４つの方法を以下に紹介します．なお，`U1`，`U2`および`U3`は nalgebra が提供する型レベル整数です．

```rust
// 全ての要素を０で初期化する
let m0 = MatrixMN::<f64, U2, U3>::zeros();

// 指定されたスライスで行優先に要素を初期化
let m1 = MatrixMN::<f64, U2, U3>::from_row_slices(&[
    1.0, 2.0, 3.0,
    4.0, 5.0, 6.0
]);

// 指定された行ベクトル（行数が１の行列）で行を初期化
let m2 = MatrixMN::<f64, U2, U3>::from_rows(&[
    MatrixMN::<f64, U1, U3>::from_row_slices(&[1.0, 2.0, 3.0]),
    MatrixMN::<f64, U1, U3>::from_row_slices(&[4.0, 5.0, 6.0]),
]);

// 全ての要素を任意の分布で初期化する．ここでは，標準正規分布を用います．
use rand;
use rand_distr::StandardNormal;
let m3 = MatrixMN::<f64, U2, U3>::from_distribution(
    &StandardNormal,
    &mut rand::thread_rng(),
)
```

### 型レベル整数による次元の指定
上述の例のように，型レベル整数を用いることで任意次元の行列を利用することができます．しかしながら，現状では以下の様な場合に対応することは出来ません．

- nalgebra が提供している最大値(`U128`)より大きな値を使いたい
- 型レベル整数に対する四則演算の結果を次元として扱いたい

そこで，型レベル整数同士の四則演算によってこの問題を解決します．この機能も nalgebra によって(正確には nalgebra が利用している typenum によって)提供されています．今回の MLP 実装では，積を表す `DimProd` を MNIST の画像を表現するために，和を表現する `DimSum` を Affine レイヤーの重みを表現するために用います．

以下に `DimProd` と `DimSum` を使用した行列の宣言の例を示します．

```rust
// 5(2 + 3)行，6(2 * 3) 列の行列を定義
type Matrix5x6 = MatrixMN::<f64, DimSum<U2, U3>, DimProd<U2, U3>>;
```

### 行列を引数として受け取る関数の定義
これまでの解説で，解析の対象となる任意の次元を持つ行列を宣言することが出来ました．そうすると，これらに対して関数によって処理を適用したくなります．行列の次元 `R` と `C` とをジェネリクスを用いて受け取る関数の例を以下に示します．

```rust
use na::{DimName, DefaulyAllocator, allocator::Allocator}

fn identity<R, C>(x: MatrixMN<f64, R, C>) -> MatrixMN<f64, R, C>
where
    R: DimName,
    C: DimName,
    DefaultAllocator: Allocator<f64, R, C>
{
    x
}
```

上記の例において，`R` と `C` は次元を表す型レベル整数でなければなりません．従って，次元を表す型レベル整数を表現する `DimName` をトレイト境界追加します．また，要素の型が `f64` で `R` 行 `C` 列の行列が確保可能であることをコンパイラに伝える必要があります．そのため，`DefaultAllocator` に対して `Allocator` をトレイト境界を追加します．

次に戻り値の行列の次元が，実引数の行列の次元に依存する場合を考えます．以下に示す例は，渡された行列の行数を`+1`，列数を`*2`した次元の零行列を返す関数です．

```rust
use na::{DimName, DimSum, DimProd, DimAdd, DimMul, DefaulyAllocator, allocator::Allocator}

fn expand_zeros<R, C>(x: MatrixMN<f64, R, C>) -> MatrixMN<f64, DimSum<R, U1>, DimProd<C, U2>>
where
    R: DimName + DimAdd<U1>,
    C: DimName + DimMul<U2>,
    DimSum<R, U1>: DimName,
    DimProd<C, U2>: DimName,
    DefaultAllocator: Allocator<f64, R, C> + Allocator<f64, DimSum<R, U1>, DimProd<C, U2>>
{
    MatrixMN<f64, DimSum<R, U1>, DimProd<C, U2>>::zeros()
}
```

この例では，前述した `DimSum` と `DimProd` を利用して戻り値の行列の次元を計算しています．これらの計算結果も同様に行列の次元を表します．従って，型レベル整数である必要があります．そのため，`DimName` がそれらのトレイト境界に追加されます．

さらに， `DimAdd` と `DimMul` とがそれぞれ `R` と `C` のトレイト境界に追加されていることが確認出来ます．これらは，それぞれ指定された型レベル整数と和および積が計算可能であることを表しています．

また，`DefaultAllocator` にも新しい制約が課せられていることが解ります．これは，戻り値の行列が確保可能であること表していいます．

## MLP の実装

### 特徴量の定義
まず初めに，MLP が取り扱う特徴量を表す型を定義します．ここでは，バッチサイズを `N` ，特徴量の次元を `C` とします．

```rust
type Feature<N, C> = MatrixMN<f64, N, C>;
```

### レイヤーの定義
今回の実装では，複数のレイヤーを直列に接続することで MLP を構成します．従って，各レイヤーに共通した機能を抽象化することで，実装の簡略化が期待できます．そこで，これらの機能を以下に示す Layer トレイトにまとめました．

```rust
trait Layer {
    // レイヤーの入力型
    type Input;
    // レイヤーの出力型
    type Output;

    // 順伝搬処理
    fn forward(&mut self, input: Self::Input) -> Self::Output;
    // 逆伝搬処理
    fn backward(&mut self, output: Self::Output) -> Self::Input;
    // パラメータの更新
    fn update(&mut self, _lr: f64, _momentum: f64) {}
}
```

以降，実際にレイヤーを実装していきます

#### Relu
以降において，$i$を特徴量のインデックス，$j$を次元のインデックス，$x_{i,j}$をレイヤーへの入力，$y_{i,j}$をレイヤーの出力，$L$を損失関数の結果とします．すると，Relu レイヤーの順伝搬$y_{i,j}$と逆伝搬$\left(\frac{\partial L}{\partial x}\right)_{i,j}$は以下の様に定義されます．

順伝搬

```math
y_{i,j} =  \begin{cases}
0  & (x_{i,j} \leq 0)\\
x_{i,j} & (x_{i,j} \gt 0)
\end{cases}
```

逆伝搬

```math
\left(\frac{\partial L}{\partial x}\right)_{i,j} = \begin{cases}
0  & (x_{i,j} \leq 0)\\
\left(\frac{\partial L}{\partial y}\right)_{i,j} & (x_{i,j} \gt 0)
\end{cases}
```

なお，Relu レイヤーには更新するパラメータが存在しないため，`update`　はデフォルト実装を用います．実装は以下の様になります．

```rust
struct Relu<N: DimName, C: DimName>
where
    DefaultAllocator: Allocator<f64, N, C>,
{
    input: Option<Feature<N, C>>,
}

impl<N: DimName, C: DimName> Relu<N, C>
where
    DefaultAllocator: Allocator<f64, N, C>,
{
    fn new() -> Self {
        Self { input: None }
    }
}

impl<N: DimName, C: DimName> Layer for Relu<N, C>
where
    DefaultAllocator: Allocator<f64, N, C>,
{
    type Input = Feature<N, C>;
    type Output = Feature<N, C>;

    fn forward(&mut self, input: Self::Input) -> Self::Output {
        let output = input.map(|v| v.max(0.0));
        self.input = Some(input);
        output
    }
    fn backward(&mut self, dldy: Self::Output) -> Self::Input {
        self.input
            .as_ref()
            .unwrap()
            .zip_map(&dldy, |x, d| if 0.0 < x { d } else { 0.0 })
    }
}
```

#### Affine

Softmax レイヤーの順伝搬 $y_{i,j}$ と逆伝搬$\left(\frac{\partial L}{\partial x}\right)\_{i,j} $ は以下の様に定義されます．ここで，$C_{in}$は入力図元を， $C_{out}$ は出力次元を，$w_{i,j}$はAffineレイヤーのパラメータを表します．

順伝搬処理

```math
y_{i,j} = \sum_{k}^{C_{in}} x_{i,k} w_{k,j} + w_{C_{in}+1,j}
```

逆伝搬処理

```math
\left(\frac{\partial L}{\partial x}\right)_{i,j} = \sum_{k}^{C_{out}} w_{j, k} \left(\frac{\partial L}{\partial y}\right)_{i,k }
```

次に，パラメータの勾配$\left(\frac{\partial L}{\partial w}\right)_{i,j}$を以下に示します．

```math
\left(\frac{\partial L}{\partial w}\right) =  \begin{cases}
    \sum_{k} \left(\frac{\partial L}{\partial y}\right)_{k,j} & (i = C_{in} + 1) \\
    \sum_{k} x_{k,j} \left(\frac{\partial L}{\partial y}\right)_{k,j} & (\mathrm{otherwise})
\end{cases}
```

また，`update` には，Momentumによる方法を，パラメータの初期化方法としてHeによる初期化を用いました．Affineレイヤーの実装は以下の様になります．

```rust
fn to_homogeneous<N: DimName, C: DimName + DimAdd<U1>>(
    input: Feature<N, C>,
) -> Feature<N, DimSum<C, U1>>
where
    DefaultAllocator: Allocator<f64, N, C> + Reallocator<f64, N, C, N, DimSum<C, U1>>,
{
    input.insert_column(C::dim(), 1.0)
}

fn he_initialize<C0: DimName + DimAdd<U1>, C1: DimName>() -> MatrixMN<f64, DimSum<C0, U1>, C1>
where
    DimSum<C0, U1>: DimName,
    DefaultAllocator: Allocator<f64, DimSum<C0, U1>, C1>,
{
    let c = (2.0 / C0::dim() as f64).sqrt();
    c * MatrixMN::<f64, DimSum<C0, U1>, C1>::from_distribution(
        &StandardNormal,
        &mut rand::thread_rng(),
    )
}

struct Affine<N: DimName, C0: DimName + DimAdd<U1>, C1: DimName>
where
    DefaultAllocator: Allocator<f64, N, DimSum<C0, U1>> + Allocator<f64, DimSum<C0, U1>, C1>,
{
    input: Option<Feature<N, DimSum<C0, U1>>>,
    d_w: MatrixMN<f64, DimSum<C0, U1>, C1>,
    v: MatrixMN<f64, DimSum<C0, U1>, C1>,
    w: MatrixMN<f64, DimSum<C0, U1>, C1>,
}

impl<N: DimName, C0: DimName + DimAdd<U1>, C1: DimName> Affine<N, C0, C1>
where
    DimSum<C0, U1>: DimName,
    DefaultAllocator: Allocator<f64, N, DimSum<C0, U1>> + Allocator<f64, DimSum<C0, U1>, C1>,
{
    fn new() -> Self {
        Self {
            input: None,
            d_w: MatrixMN::<f64, DimSum<C0, U1>, C1>::zeros(),
            v: MatrixMN::<f64, DimSum<C0, U1>, C1>::zeros(),
            w: he_initialize::<C0, C1>(),
        }
    }
}

impl<N: DimName, C0: DimName + DimAdd<U1>, C1: DimName> Layer for Affine<N, C0, C1>
where
    DimSum<C0, U1>: DimName,
    DefaultAllocator: Allocator<f64, N, C0>
        + Allocator<f64, C1, C0>
        + Allocator<f64, N, C1>
        + Allocator<f64, DimSum<C0, U1>, N>
        + Allocator<f64, DimSum<C0, U1>, C1>
        + Reallocator<f64, N, C0, N, DimSum<C0, U1>>,
{
    type Input = Feature<N, C0>;
    type Output = Feature<N, C1>;

    fn forward(&mut self, input: Self::Input) -> Self::Output {
        let input = to_homogeneous(input);
        let output = &input * &self.w;
        self.input = Some(input);
        output
    }
    fn backward(&mut self, dldy: Self::Output) -> Self::Input {
        self.d_w = self.input.as_ref().unwrap().transpose() * &dldy;
        dldy * self.w.fixed_slice::<C0, C1>(0, 0).transpose()
    }

    fn update(&mut self, lr: f64, momentum: f64) {
        self.v = momentum * &self.v - lr * &self.d_w;
        self.w += &self.v;
    }
}
```

#### Softmax
Softmax レイヤーの順伝搬$y_{i,j}$と逆伝搬$\left(\frac{\partial L}{\partial x}\right)_{i,j}$は以下の様に定義されます．

順伝搬処理

```math
y_{i,j} = \frac{\exp\left(x_{i,j} - \max_{k}x_{i,k}\right)}{\sum_{l}\exp\left(x_{i,l} - \max_{k} x_{i,k}\right)}
```

逆伝搬処理

```math
\left(\frac{\partial L}{\partial x}\right)_{i,j} = y_{i,j} \left(\left(\frac{\partial L}{\partial y}\right)_{i,j}  - \sum_{k} y_{i,k}  \left(\frac{\partial L}{\partial y}\right)_{i,k} \right)
```

Softmax レイヤーも更新するパラメータが存在しないため，`update`　はデフォルト実装を用います．Softmax レイヤーの実装は以下の様になります．

```rust
struct Softmax<N: DimName, C: DimName>
where
    DefaultAllocator: Allocator<f64, N, C>,
{
    input: Option<Feature<N, C>>,
}

impl<N: DimName, C: DimName> Softmax<N, C>
where
    DefaultAllocator: Allocator<f64, N, C> + Allocator<f64, U1, C> + Allocator<f64, C, U1>,
{
    fn new() -> Self {
        Self { input: None }
    }
    fn _forward(input: &Feature<N, C>) -> Feature<N, C> {
        Feature::<N, C>::from_rows(
            &input
                .row_iter()
                .map(|row| row.map(|v| (v - row.max()).exp()))
                .map(|row| &row / row.sum())
                .collect::<Vec<_>>(),
        )
    }
}

impl<N: DimName, C: DimName> Layer for Softmax<N, C>
where
    DefaultAllocator: Allocator<f64, N, C> + Allocator<f64, U1, C> + Allocator<f64, C, U1>,
{
    type Input = Feature<N, C>;
    type Output = Feature<N, C>;

    fn forward(&mut self, input: Feature<N, C>) -> Feature<N, C> {
        let output = Self::_forward(&input);
        self.input = Some(input);
        output
    }
    fn backward(&mut self, dldy: Feature<N, C>) -> Feature<N, C> {
        let output = Self::_forward(self.input.as_ref().unwrap());
        Feature::<N, C>::from_rows(
            &dldy
                .row_iter()
                .zip(output.row_iter())
                .map(|(d, o)| o.component_mul(&d.add_scalar(-(d * o.transpose())[0])))
                .collect::<Vec<_>>(),
        )
    }
}

```

### ネットワークの定義
ネットワークは以下に示す`Network`トレイトを実装することで実現します．また，推論は`forward`によって，学習は`forward`，`backward`, `update`の繰り返しによって実現します．

```rust
trait Network {
    // ネットワークの入力型
    type Input;
    // ネットワークの出力型
    type Output;

    // 順伝搬処理
    fn forward(&mut self, input: Self::Input) -> Self::Output;
    // 逆伝搬処理
    fn backward(&mut self, output: Self::Output) -> Self::Input;
    // 各レイヤーのパラメータの更新
    fn update(&mut self, lr: f64, momentum: f64);
}
```

### マクロによるネットワークの実装
今回の実装では，マクロによってレイヤーを接続しネットワークを構築します．これは，型の異なる複数のレイヤーを上手く取り扱う方法が見つからなかったためです．以下に各種マクロの定義を示します．



```rust
#[macro_export]
macro_rules! impl_forward {
    ($self:expr, $input:expr, $head:ident,) => {
        $self.$head.forward($input)
    };
    ($self:expr, $input:expr, $head:ident, $($tail:ident,)*) => {
        impl_forward!($self, $self.$head.forward($input), $($tail,)*)
    };
}

#[macro_export]
macro_rules! impl_backward {
    ($self:expr, $input:expr, ) => {$input};
    ($self:expr, $input:expr, $head:ident, $($tail:ident,)*) => {
        $self.$head.backward(
            impl_backward!($self, $input, $($tail,)*)
        )
    };
}

#[macro_export]
macro_rules! impl_update {
    ($self:expr, $lr:expr, $momentum:expr, ) => {};
    ($self:expr, $lr:expr, $momentum:expr, $head:ident, $($tail:ident,)*) => {
        $self.$head.update($lr, $momentum);
        impl_update!($self, $lr, $momentum, $($tail,)*)
    };
}

#[macro_export]
macro_rules! first_layer {
    (($head:ty, $($_:ty,)*)) => {
        $head
    };
}

#[macro_export]
macro_rules! last_layer {
    (($head:ty,)) => {
        $head
    };
    (($head:ty, $($tail:ty,)*)) => {
        last_layer!(($($tail,)*))
    };
}

#[macro_export]
macro_rules! define_network {
    ($network_name:ident, $(($name:ident, $layer:ty)),*) => {
        struct $network_name{$($name:$layer,)*}
        impl $network_name {
            pub fn new() -> Self {
                Self{$($name:<$layer>::new(),)*}
            }
        }
        impl Network for $network_name {
            type Input = <first_layer!(($($layer,)*)) as Layer>::Input;
            type Output = <last_layer!(($($layer,)*)) as Layer>::Output;

            fn forward(&mut self, input: Self::Input) -> Self::Output{
                impl_forward!(self, input, $($name,)*)
            }
            fn backward(&mut self, output: Self::Output) -> Self::Input{
                impl_backward!(self, output,$($name,)*)
            }
            fn update(&mut self, lr: f64, momentum: f64) {
                impl_update!(self, lr, momentum, $($name,)*);
            }
        }
    };
}
```

実際にネットワークのを定義するには `define_network` マクロを使用します．以下に，AffineレイヤーとSoftmaxレイヤーからなる単純なネットワークを定義する例を示します．

```rust
define_network!(
    SimpleNetwork,  // ネットワークの型名
    (fc0, Affine<U1, U8, U2>), // fc0という名称でAffineレイヤーを定義
    (softmax0, Softmax<U1, U2>) // softmax0という名称でSoftmaxレイヤーを定義
);
```

そして，このマクロは以下の様に展開されます．

```rust
struct SimpleNetwork {
    fc0: Affine<U1, U8, U2>,
    softmax0: Softmax<U1, U2>,
}
impl SimpleNetwork {
    pub fn new() -> Self {
        Self {
            fc0: <Affine<U1, U8, U2>>::new(),
            softmax0: <Softmax<U1, U2>>::new(),
        }
    }
}
impl Network for SimpleNetwork {
    type Input = <Affine<U1, U8, U2> as Layer>::Input;
    type Output = <Softmax<U1, U2> as Layer>::Output;
    fn forward(&mut self, input: Self::Input) -> Self::Output {
        self.softmax0.forward(self.fc0.forward(input))
    }
    fn backward(&mut self, output: Self::Output) -> Self::Input {
        self.fc0.backward(self.softmax0.backward(output))
    }
    fn update(&mut self, lr: f64, momentum: f64) {
        self.fc0.update(lr, momentum);
        self.softmax0.update(lr, momentum);
    }
}
```

以上の結果より，`SimpleNetwork`は`Network`トレイトを実装していることが確認出来ます．また，`forward`では各レイヤーの`forward`が，`backward`では各レイヤーの`backward`が，`update`では各レイヤーの`update`が呼ばれていることが確認できます．

## 学習と評価の実装

これまでに実装したレイヤーを用いて，３層の全結合層からなるネットワーク`MNISTNetwork`を定義します．なお，全結合層の次元は80としました．`define_network`マクロを用いた定義を以下に示します．

```rust
type N = U64; // バッチサイズ
type F = DimProd<U28, U28>; // 入力次元
type M = U80; // 中間層の次元
type C = U10; // 出力次元（識別クラス数）

define_network!(
    MNISTNetwork,
    (fc0, Affine<N, F, M>),
    (act0, Relu<N, M>),
    (fc1, Affine<N, M, M>),
    (act1, Relu<N, M>),
    (fc2, Affine<N, M, C>),
    (softmax0, Softmax<N, C>)
);
```

次に，ネットワークの学習に用いる関数`load_dataset`，`cross_entropy`，`accuracy` を定義します．`load_dataset`は引数として渡された画像とラベルのスライスをネットワークに入力可能なフォーマットに整形します．`load_dataset`の実装を以下に示します．

```rust
fn load_dataset<N, F, C>(
    img: &[f32],
    label: &[u8],
) -> impl Iterator<Item = (Feature<N, F>, Feature<N, C>)>
where
    N: DimName, // バッチサイズ
    F: DimName, // 画像の次元
    C: DimaName, // ラベルの次元
    DefaultAllocator: Allocator<f64, N, F> + Allocator<f64, N, C>,
{
    fn create_mini_batches<N: DimName, C: DimName, T>(feature: &[T]) -> Vec<Feature<N, C>>
    where
        N: DimName, // バッチサイズ
        C: DimName, // 特徴量の次元
        T: AsPrimitive<f64>, // 特徴量要素の型
        DefaultAllocator: Allocator<f64, N, C>,
    {
        feature
            .into_iter()
            .map(|&v| v.as_())
            .collect::<Vec<_>>()
            .chunks_exact(N::dim() * C::dim()) // N * C毎にサンプルをまとめる
            .map(|chunk| Feature::<N, C>::from_row_slice(chunk))
            .collect::<Vec<_>>()
    }
    let x = create_mini_batches::<N, F, f32>(img);
    let y = create_mini_batches::<N, C, u8>(label);

    x.into_iter().zip(y.into_iter()) 
}
```

`cross_entropy`は名前の通り，ネットワークの出力と対応するラベルからクロスエントロピーを計算します．`cross_entropy`の実装を以下に示します．

```rust
fn cross_entropy<N: DimName, C: DimName>(
    x: &Feature<N, C>,
    y: &Feature<N, C>,
) -> (f64, Feature<N, C>)
where
    DefaultAllocator: Allocator<f64, N, C>,
{
    let n = N::dim() as f64;
    let epsillon = 1e-12;
    let forward = -x.zip_map(y, |_x, _y| _y * (_x + epsillon).ln()).sum() / n;
    let backward = x.zip_map(y, |_x, _y| -_y / (_x + epsillon)) / n;
    (forward, backward)
}
```

`accuracy`も名前の通り，ネットワークの出力と対応するラベルから精度を計算します．`accuracy`の実装を以下に示します．

```rust
fn accuracy<N: DimName, C: DimName>(y_hat: &Feature<N, C>, y: &Feature<N, C>) -> f64
where
    DefaultAllocator: Allocator<f64, N, C> + Allocator<f64, C, U1>,
{
    let acc = y_hat
        .row_iter()
        .zip(y.row_iter())
        .map(|(_y_hat, _y)| _y_hat.transpose().imax() == _y.transpose().imax())
        .fold(0, |acc, m| acc + m as u32);
    (acc as f64) / (N::dim() as f64)
}
```

それでは，実際にネットワークの学習を行っていきます．まず初めに，データセットのダウンロードと読み込みを行います．これには，[mnist](https://crates.io/crates/mnist)というクレートが提供する`MnistBuilder`を利用します．`MnistBuilder`はBuilderパターンを用いて，データ数やフォーマットを指定します．今回は，以下の条件でデータのロードを行います．

* トレーニング画像50000枚
* バリデーション画像10000枚
* テスト画像10000枚
* ラベルをone-hotエンコーディング
* 画像を正規化

従って，`MnistBuilder`は以下の様に設定されます．

```rust
    let NormalizedMnist {
        trn_img, // トレーニング画像
        trn_lbl, // トレーニングラベル
        val_img, // バリデーション画像
        val_lbl, // バリデーションラベル
        tst_img, // テスト画像
        tst_lbl, // テストラベル
    } = MnistBuilder::new()
        .label_format_one_hot()
        .training_set_length(50000)
        .validation_set_length(10000)
        .test_set_length(10000)
        .download_and_extract()
        .finalize()
        .normalize();
```

次に，ネットワークの学習を行います．ネットワークの学習は，エポック毎にトレーニング画像とトレーニングラベルによるパラメータの更新と，バリデーション画像とバリデーションラベルによる汎化性能の検証を行います．具体的なコードを以下に示します．なお，損失関数には`cross_entropy`を利用します．

```rust
    let mut network = MNISTNetwork::new();
    let epochs = 12;
    let lr = 0.01;
    let momentum = 0.9;

    for _ in 0..epochs {
        let train_loss = {
            let (acc_train_loss, n) = load_dataset::<N, F, C>(&trn_img, &trn_lbl).fold(
                (0.0, 0),
                |(acc, n), (train_x, train_y)| {
                    // 順伝搬結果とラベルから損失とその勾配を計算
                    let (loss, delta) = cross_entropy(&network.forward(train_x), &train_y);
                    // 逆伝搬
                    network.backward(delta);
                    // パラメータ更新
                    network.update(lr, momentum);
                    (acc + loss, n + 1) // ロスとバッチ数の積算
                },
            );
            acc_train_loss / (n as f64)
        };

        let val_loss = {
            let (acc_val_loss, n) = load_dataset::<N, F, C>(&val_img, &val_lbl).fold(
                (0.0, 0),
                |(acc, n), (val_x, val_y)| {
                    // 順伝搬結果とラベルから損失を計算（パラメータの更新は行わないので勾配は不要）
                    let (loss, _) = cross_entropy(&network.forward(val_x), &val_y);
                    (acc + loss, n + 1)
                },
            );
            acc_val_loss / (n as f64)
        };
        println!("train_loss :{:?}, val_loss:{:?}", train_loss, val_loss);
    }
```

最後にテスト画像とテストラベルを用いて，ネットワークの性能を検証します．今回，性能指標として`accuracy`を用いました．

```rust
    let accuracy = {
        let (acc_accuracy, n) = load_dataset::<N, F, C>(&tst_img, &tst_lbl)
            .fold((0.0, 0), |(acc, n), (tst_x, tst_y)| {
                (acc + accuracy(&network.forward(tst_x), &tst_y), n + 1)
            });
        acc_accuracy / (n as f64)
    };
    println!("accuracy: {:?}", accuracy);
```

## 分類器の学習と評価
これまでに実装したコードを用いて，分類器の学習と評価を行います．プロジェクト構成を以下に示します．

```bash
$ tree
.
├── Cargo.toml
└── src
    └── main.rs

1 directory, 2 files
```

`main.rs`は上述のコードを一つのファイルにまとめたものです．[こちら](https://gist.github.com/ar90n/eeee7976b2c416e7100b29189283d69d#file-main-rs)から取得可能です．`Cargo.toml`は以下の通りです．

```toml
[package]
name = "neural_network_rust_nalgebra"
version = "0.1.0"
authors = ["Masahiro Wada <masahiro_wada@caddi.jp>"]
edition = "2018"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html


[dependencies]
nalgebra = '*'
approx = '*'
rand = "*"
rand_distr = "*"
num = "*"

[dependencies.mnist]
version = "*"

features = ["download"]
```

実行結果を以下に示します．

```bash
$ cargo run --release
    Finished release [optimized] target(s) in 0.34s
     Running `target/release/neural_network_rust_nalgebra`
Attempting to download and extract train-images-idx3-ubyte.gz...
  File "data/train-images-idx3-ubyte.gz" already exists, skipping downloading.
  Extracted file "data/train-images-idx3-ubyte" already exists, skipping extraction.
Attempting to download and extract train-labels-idx1-ubyte.gz...
  File "data/train-labels-idx1-ubyte.gz" already exists, skipping downloading.
  Extracted file "data/train-labels-idx1-ubyte" already exists, skipping extraction.
Attempting to download and extract t10k-images-idx3-ubyte.gz...
  File "data/t10k-images-idx3-ubyte.gz" already exists, skipping downloading.
  Extracted file "data/t10k-images-idx3-ubyte" already exists, skipping extraction.
Attempting to download and extract t10k-labels-idx1-ubyte.gz...
  File "data/t10k-labels-idx1-ubyte.gz" already exists, skipping downloading.
  Extracted file "data/t10k-labels-idx1-ubyte" already exists, skipping extraction.
train_loss :0.3990251154045798, val_loss:0.20073474741457015
train_loss :0.18740507891471156, val_loss:0.15270709888764708
train_loss :0.1384824042466224, val_loss:0.13499420871190054
train_loss :0.11040477926959635, val_loss:0.12465709407120212
train_loss :0.09132263514236212, val_loss:0.11778158499482444
train_loss :0.0768990102236621, val_loss:0.11111220500240461
train_loss :0.06550682241861264, val_loss:0.10474639523700195
train_loss :0.05632068071386103, val_loss:0.09988697537468595
train_loss :0.048288227996312404, val_loss:0.09742129201130649
train_loss :0.041381218803573774, val_loss:0.09715223463915591
train_loss :0.03555275317739821, val_loss:0.09735193014844665
train_loss :0.030610562793088675, val_loss:0.0989283094945219
accuracy: 0.9735576923076923
```

以上の結果より，約97.4%の精度が得られていることが確認出来ました．

## まとめ

- Rust と nalgebra を用いて MLP を実装しました
- MNIST にて 約97.4%の精度が得られました

## 最後に

キャディでは，フロントエンド，バックエンド，アルゴリズム，SREと幅広くエンジニアを募集しております．興味をお持ちいただけた方は，[こちら](https://share.hsforms.com/1XUTUewm3SP2xrN6MKJNg_Q3x4tb)からご連絡ください！

明日は，@m-yskさんによる「Rust入門者がrust-analyzerにコントリビュートするまで」です!!


## 参考文献

1. 斎藤康毅(2016)，ゼロから作る Deep Learning: Python で学ぶディープラーニングの理論と実装，オライリー・ジャパン，