Keyboard shortcuts

Press ← or → to navigate between chapters

Press S or / to search in the book

Press ? to show this help

Press Esc to hide this help

はじめに

近年、人工知能は私たちの生活に溶け込み、様々な利便性をもたらしてきました。対話型のAIなど、今ではAIは仕事やプライベートといった際に欠かせない存在となってきています。そんなAIにたよる日常生活が当たり前になりつつあるこの時代において、AIの信憑性や、人の学習能力の低下、情報の操作、そして情報漏洩といったAIの危険性が重視され始めています。こういった危険性が顕在化している背景に、AIの急速な大衆化があると考えます。AIの研究者とは違い、未知のものである一般の人々がインターネットやメディアにより、ここ数年であっという間にAIというものが身近なものとなり、簡単に利用できるようになりました。しかし、人々は突然現れた、日常で欠かせない得体の知れないAIというものを「ただ便利だから」、「みんなが使っているから」、ということで十分理解せずに利用しようとしているのではないでしょうか。人々の日常生活に溶け込む人工知能はもはや、現代社会において教養として学ぶ必要があります。なんとなく文を送ったら、返事が来たというような、裏側を見ない表面的なAIの理解だけでの使用は、浅はかで、無責任なAIの使用と言えるのではないでしょうか。AIはどのような仕組みで構成されているのか、裏側でどのような処理がされているのか、そういったAIの裏側を学ぶことで、私たちは真のAIとの共存が可能になるでしょう。このドキュメントでは、現代のAIの基礎である深層学習の仕組みを、新しいプログラミング言語Rustで一から実装するプロセスを通じて、AIの裏側を学びます。高校生である筆者らにとって、この開発は大胆な挑戦であり、多少荒々しい冒険となりますが、ぜひ最後までお付き合いください。

フレームワーク StuCrs のコンセプト

現在、深層学習のフレームワークはTensorFlowやPytorchというように複数存在します。それらのフレームワークには様々な特徴、コンセプトがありますが、そんな中で私たちは一からRustで深層学習のフレームワーク： StuCrs を開発、実装しました。

TensorFlowやPyTorchといったフレームワークはオープンソースであり、公開されています。試しにGitHubでTensorFlowの中のコードを見てみましょう。するとどうでしょうか、ファイルやフォルダーがたくさん存在し非常に複雑に構成されています。また、ファイルの中のコードもパフォーマンスを最適化するために、様々な関数が組み合わさって、読み手にとって解読しづらいものとなっています。既存の多くのフレームワークは様々なプロジェクトでフレームワークとして利用されるためにパフォーマンスや機能の豊富さ、バグの少なさなどが重視されます。そこで私たちは逆の立場として、シンプルな構造のフレームワークの構築を目指し、それを読者が理解して1からフレームワークを実装してもらうことで深層学習の原理を探究してもらう、そんな目標をもってこのような研究に至りました。

また、Rustという言語を選択した理由として、Rust言語のコンピューターに寄り添う仕様にあります。近年人気の高いpython言語はAIの分野で高い地位を得ていますが、それはその言語が人間に理解しやすい言語であり、設計がわかりやすいからです。しかし、そのかりやすさが故にコンピューターの本来の特性を十分理解せずに利用できてしまうという特徴もあります。一方Rust言語はコンピューターに寄り添った言語で、確かに人間にとって理解しづらい言語です。しかしながら、Rustを一つ一つ自らの手で実装することで、コンピューターとの真の対話をし、本当のプログラミング、そしてコンピューターによるAIの計算を真髄から学ぶことができるのではないでしょうか。私たちはコンピューター・AIの計算処理のアルゴリズムを実際に実装し、学ぶ教育ツールの一つとしてこのようなアプローチを提案します。

どのような人を対象としているか

このドキュメントを見ている方はすでにこの研究がどのようなものかある程度把握されていると思いますが、改めてどのような方のためか簡単に示します。

  • Rust言語を学ぶためにのRustの練習のサンプルとして利用したい方。
  • 深層学習のフレームワークを一から実装し、原理から学びたい方。
  • Rustでの機械学習、深層学習の分野の開発に興味がある方。

学ぶに当たってあらかじめ必要な知識

このドキュメントは一から実装するにあたって、学びたいと考えている様々なユーザーが理解できるよう分かりやすく書きましたが、十分理解して進めるためにあらかじめいくつかの知識が必要だと考えています。

  • Rust言語について基本的な仕様を理解していること。

    このドキュメントではRustの特徴と深層学習の知識をなるべく多く説明していますが、主に深層学習の原理を学ぶことをこのドキュメントの第一の目的としているため、多少高度なRustの知識の説明を省略していることがあります。Rustの仕様をあらかじめ理解してもらうのが良いと思います。ただし、このドキュメントで使用するRustのコードは検索すれば簡単に理解できる基礎的なものを使用しているので、特に高度なスキルは必要ありません。

  • 深層学習の概要をあらかじめ理解していること。

    このドキュメントは深層学習の原理を探究するものですが、一から実装するという高度な開発を行うため、事前に深層学習の全体像を把握しておいた方が理解が進みます。

  • 数学の知識。

    AIの分野ではとりわけ数学が多く応用されています。このドキュメントの場合、主に微分と行列の知識が必要です。例えば、sin,cosの微分はもちろん、tanhや偏微分、行列積などといった高度な数学をモデル化します。

ドキュメントの進め方

はじめにバイナリクレートを一つ用意してください。このクレートはプログラムをテストするためのファイルですので、名前は何でも構いません。（test_⋯.rsなど）はじめのうちはこのファイルのmain()の外に構造体やトレイトを定義し、プログラムをこのファイルのmain()で書き、そこで実行します。ある程度関数や構造体が増え、複雑になってきたら、ライブラリクレートをを新たに作成します。このクレート名はstucrsという感じにフレームワークの名前にするのがおすすめです。そこに今まで作成した構造体たちを移し、ライブラリクレートからテスト用ファイルにインポートする形を取ります。なお、移行するタイミングやその方法は記載してあります。

注意

このドキュメントはPCでの閲覧を推奨しています。グラフや数式が正しく表示されない可能性があります.

環境

本研究のリポジトリの中に含まれるDockerfile、compose.yamlファイルを用いてコンテナを立ち上げてください。(Dockerに関する使い方の説明は省略します。)
注意 DockerfileのFROM のところは自分の環境に適したCUDAのバージョンを指定してください。また、NVIDIA製のGPUを使用しない場合、ubuntu20.04をお使いください。

docker build -t cuda-im ./ #イメージ名はcuda-im
docker compose up -d # イメージ名をcuda-imと設定しているので、イメージ名を変更したい場合はyamlファイルの中を変更してください。

環境の詳細はGitHubのREADMEをご覧ください。

変数の作成

変数とはデータを格納する箱のようなものです。はじめにこの変数をVariableという名前で構造体として実装してみましょう。

fn main() {
    let mut a = Variable::new(1.0);
    println!("{}", a.data); //1.0
    a.data = 5.0;
    println!("{}", a.data); //5.0
}

struct Variable {
    data: f32,
}

impl Variable {
    fn new(data: f32) -> Self {
        Variable { data }
    }
}

まず Variable という構造体を実装します。フィールドとして小数をあつかえるf32型をdataとして保持します。また、コンストラクタを生成するための初期化の関数 new() を定義します。これにより、f32のデータを渡すことでデータを保持する変数、Variable型を生成することができます。main関数を実行すると、Variableのデータを見ることができます。

pythonを学んだ人からすると、関数の戻り値がしっかり示されていないと感じるかもしれません。pythonだと戻り値はreturnを用いますが、Rustの場合、戻り値となるものには後ろにセミコロンをつけないのが通常です。この場合、fn new()の戻り値はVariable{data}ですが、後ろにセミコロンがないので戻り値として扱われます。

関数の作成

先ほど変数を実装しましたが、変数はどのように作り出されるのでしょうか。それは 関数 です。変数を入力として関数に渡し、出力として新たな変数を生み出します。それでは関数を実装していきましょう。

Functionトレイトの実装

はじめにRust独自の仕様であるトレイトを理解する必要があります。トレイトとは様々な構造体がある特定のふるまいをすることを保証するものです。ここでいうふるまいとは、トレイトを継承した構造体がすべて同じある関数を保持しているということです。　　 今後関数は様々な構造体で使用するのでトレイトを使用して関数を定義します。では実際に関数をトレイトで実装してみましょう

fn main() {
 　let x = Variable::new(2.0);
    　println!("{}", x.data); //2.0
   　 let f = Square {};
    　let y = f.call(&x);
   　 println!("{}", y.data); //4.0
 }

trait Function {
    fn call(&self, input: &Variable) -> Variable;
    fn forward(&self, x: f32) -> f32;
}

struct Square {}

impl Function for Square {

    fn call(&self, input: &Variable) -> Variable {
        let x = input.data;
        let y = self.forward(x);
        let output = Variable::new(y);
        output
    }

    fn forward(&self, x: f32) -> f32 {
        x.powf(2.0)
    }
}

まずはFunctionというtraitを実装します。ここでは call と forward というメソッドを作ります。今後このtraitには多くの関数（Exp関数やsin関数など）を追加するためcallにはすべての関数に共通する「Variableからデータを取り出す」、「計算結果をVariable型にして返す」という２つの機能のみを追加します。具体的な計算はforwardにやらせます。ここで大事なのは input 、output はVariable型、x,yはf32型であるということです。

次に Squre という構造体を実装します. その後implキーワードという定義されたメソッドをimplブロック内で具体的に実装できる機能を用いることでSqure関数という２乗の計算を実装します。main関数を実行すると２の２乗の結果がでます。

Exp関数の実装

新しい関数を１つ実装します。今回はe(ネイピア数)のｘ乗という関数を実装します。

struct Exp {}

impl Function for Exp {

    fn call(&self, input: &Variable) -> Variable {
        let x = input.data;
        let y = self.forward(x);
        let output = Variable::new(y);
        output
    }
    fn forward(&self, x: f32) -> f32 {
        x.exp()
    }
}

Squre構造体と同じように実装します。変更点はなかの計算がeのｘ乗に変わっただけです。 ※今後、このExp構造体のように、Functionトレイトを実装した関数の構造体(SinやTanhなど)を多く実装していきます。その際、これらの構造体を Function構造体 と呼ぶことにします。

Function構造体を呼び出す関数

「Functionトレイトの実装」 のコードを見てみましょう。Function構造体を用いて計算する際、let f = Square{} と、 let y = f.call(&x) という二つのコードで実行しています。この処理はSquareというFunction構造体を作成し、作成した構造体に変数を渡して計算するという今後多く使われる基本的な処理です。これを一つの関数として定義しましょう。

fn square(&x:Variable) -> Variable {
    let f = Square{};
    let y = f.call(x);
    y
}

このように一つの関数としてまとめると、前のコードをよりわかりやすく書けます。

fn main() {
 　let x = Variable::new(2.0);
    　println!("{}", x.data); //2.0
    　let y = square(&x);   
   　 println!("{}", y.data); //4.0
 }

ここで大事なのは、構造体ははじめの一文字は 大文字で、関数の場合はすべて小文字で書くというルールを決めることです。そうすることで、このコードは 構造体、関数 どちらを指しているのかがすぐわかります。では同じようにExpにも実装してみましょう。

微分の理論

微分は機械学習の分野において重要なものです。私たちはこれから微分を使って機械学習の核心を実装していきます。それに先立ち、この章では微分の効率的な求め方について解説します。なお、ここでは合成関数の微分の知識を前提としています。

チェーンルール

チェーンルールとは複数の関数が組み合わさった「合成関数」を微分する際に、それぞれの関数の微分を掛け合わせることで、合成関数の微分（導関数）を求められるというものです。
式にすると

$$\frac{dy}{dx}=\frac{dy}{dy} \cdot\frac{dy}{db} \cdot\frac{db}{da} \cdot\frac{da}{dx}$$

となり、ｘに関するｙの微分は各関数の微分の積によって表わすことができます。つまり合成関数の微分は各関数の局所的な微分へ分解できるということです。

バックプロパゲーションの理論

$$\frac{dy}{dx}=\frac{dy}{dy} \cdot\frac{dy}{db} \cdot\frac{db}{da} \cdot\frac{da}{dx}$$

の式を変形すると

$$\frac{dy}{dx} = \left(\left(\frac{dy}{dy}\cdot\frac{dy}{db}\right)\frac{db}{da}\right)\frac{da}{dx}$$

この式は出力方向から入力方向へ微分を計算していくことを表しています。()で閉じたところはdy/dbになっています。 これは最初dy/dyだったのが、C’(b)のbackward（もしくは導関数）によってdy/dbを導くことができたということです。

これを同様にB’(a)、A’(x)でも行うと、dy/dxが求まります。この計算をグラフによって表すと

graph RL
 A(("$$dy/dy$$")) --> B["$$C'(b)$$"]
 B --> C(("$$dy/db$$"))
 C --> D["$$B'(a)$$"]
 D --> E(("$$dy/da$$"))
 E --> F["$$A'(x)$$"]
 F --> G(("$$dy/dx$$"))

グラフにすることで微分の流れがわかりやすく、このグラフを見るとｙの各変数に関する微分が伝播することでｘに関するｙの微分が求まっていることがわかります。これがバックプロパゲーションです。ここでの重要な点は伝播するデータはすべて「ｙに関する微分」ということです。
なぜこのバックプロパゲーションが微分の効率的な求め方なのかというと機械学習は基本的に大量のパラメータを入力として「損失関数」を求めていくものです。つまり私たちは損失関数の各パラメータに関する微分を求める必要があります。そのような場合バックプロパゲーションを用いれば一回の伝播で全てのパラメータに関する微分を求められるのです。

ここで通常の計算と微分を求める計算の比較をしてみると

graph LR
 A(("$$x$$")) --> B["$$A$$"]
 B --> C(("$$a$$"))
 C --> D["$$B$$"]
 D --> E(("$$b$$"))
 E --> F["$$C$$"]
 F --> G(("$$y$$"))

graph RL
 A(("$$dy/dy$$")) --> B["$$C'(b)$$"]
 B --> C(("$$dy/db$$"))
 C --> D["$$B'(a)$$"]
 D --> E(("$$dy/da$$"))
 E --> F["$$A'(x)$$"]
 F --> G(("$$dy/dx$$"))

この２つの図を見てわかるように順伝播と逆伝播は対応していることがわかります。そしてC‘(b)に注目するとこの微分の計算をするにはｂという値が必要になります。 同様なことが各計算にもいえます。つまり逆伝播をするためには順伝播で求めたデータが必要であるということです。 この微分の計算を最初のxまでたどれば、yをxで微分したことになります。 そのためバックプロパゲーションを実装するには順伝播を行い、その各変数の値を保持しなければなりません。

いろいろ長く話してきましたが、 大事なことは、 変数と関数のつながりとデータをしっかり保持し、この流れの源流にたどれるようにすることです。 はじめはすんなり理解できないと思いますが、 手で実装するにつれ、理解できるようになります。 詳しい説明は実際に実装する時に行います。では実際にバックプロパゲーションを実装していきましょう。

微分の実装(手動)

前の章では微分の理論について説明してきましたが、この章ではVariable構造体とFunctionトレイトを拡張して微分をもとめていきます。

Variable構造体とFunctionトレイトの拡張

Variable構造体は微分の値を保持するために、通常の値に加えてそれに対応した微分の値を持つように変更します。いままでのFunctionトレイトは通常の計算をする機能しか持っていませんでした。そのため微分の計算をするbackward機能と、逆伝播のために順伝播する際に入力された値を保持する機能を追加します。

#[derive(Debug)]

struct Variable {
    data: f32,
    grad: Option<f32>,

}
impl Variable {
    fn new(data: f32) -> Self {
        Variable { data, grad: None }
    }
}

trait Function {
    fn call(&mut self) -> Variable;
    fn forward(&self, x: f32) -> f32; // 引数f32
    fn backward(&self, gy: f32) -> f32; // 引数f32

   }

Variableはフィールドとして初期状態や勾配が不要な場合はNoneとし、逆伝播で計算された後にはSome(f32)で値を持つようにするためOption型でgradを保持させます。関数new()にgradとしてNoneを持たせgradを保持します。今後、構造体のコンストラクタを作成する際は、このnew()を使って作成します。pythonでいうinit__(self、)です。

Functionトレイトは call() も変更し、backward() を追加します。以前のcall()は &self にしていましたが、&mut self に変更します。なぜならcall()の中でFunction構造体のフィールドであるinputを変更するため、selfを可変で渡す必要があるからです。 また前のcall()にはinputを渡していましたが、今回はそれを削除します。その理由は次の4.2の変更でわかります。

Function構造体の拡張

続いて、具体的な関数の逆伝播を実装していきます。まずは2乗の計算をするSqure構造体です。ｙ＝Ｘ２の微分は２Ｘとなることから実装します。次にｙ＝eXの計算をするExp構造体です。この微分の値はe**Xとなりこれをもとに実装していきます。

今後も数学的な関数を実装していくうえで、微分の説明は省略することがあります。

/// 新しい設定
struct Square{
    input: Variable,
}
impl Function for Square {
    fn call(&mut self) -> Variable {
        let x = self.input.data; //inputのデータをフィールドから得る
        let y = self.forward(x);
        let output = Variable::new(y);
        output
    }

    fn forward(&self, x: f32) -> f32 {
        x.powf(2.0)
    }

    fn backward(&self, gy: f32) -> f32 {
        let x = self.input.data;
        2.0 * x * gy // gxを表す
    }

}

impl Square {
    fn new(input:&Variable) -> Self {
        Self { input: input.clone() }　//ここでinputをフィールドとして持つ
    }
}

fn square(&x:Variable) -> Variable {
    let f = Square::new(x);
    let y = f.call();
    y
}
// 前のcall()の設計の場合　この下のコードは以前の設計のものです。間違えないようにしてください。　

struct Square{
    input: Option<Variable>, //Option型にしなければならない
}
impl Function for Square {
    fn call(&mut self,input:&Variable) -> Variable {
        let x = input.data; //引数のinputからデータを得る
        self.input = Some(input.clone()); //ここでinputを保持
        let y = self.forward(x);
        let output = Variable::new(y);
        output
    }

    fn forward(&self, x: f32) -> f32 {
        x.powf(2.0)
    }

    fn backward(&self, gy: f32) -> f32 {
        let x = self.input.unwrap().data; //unwarp()を使わないといけない
        2.0 * x * gy // gxを表す
    }

}

impl Square {
    fn new() -> Self {
        Self { input: None }　//はじめは持っていないので、None
    }
}

fn square(&x:Variable) -> Variable {
    let f = Square::new();
    let y = f.call(x);
    y
}

call()のところはFunctionトレイトと同じように揃えます。call()の中の変更点としてinputのデータであるxの取得方法です。前はcall()にinputを渡していましたが、今回はフィールドから得ています。なぜこのようなことができるようになったかというと、2.3で実装した構造体を呼び出す関数のおかげです。この関数で構造体を呼び出す際に同時にinputを構造体にinputフィールドとしてはじめから持たせることで、call()に渡さなくても、inputにアクセスできるというわけです。

なぜこのように変更したかというと、Option型を多用しないため です。

以前の設計の場合のコードを見てみましょう。新しい関数インスタンスを作成するためのfn new() ->Selfの戻り値は Self { input: None } と、初期状態では入力が未設定であるため、inputフィールドをNoneで初期化する必要があります。するとここでinputフィールドはOption型で保持しなくてはなりません。Option型はエラーの原因となりやすく、またunwrap()など、取り出すのにコードが長くなったりと、あまり多用すべきものではありません。一方、新しい設計の場合、Option型を使わずに実装できています。

new()関数はすべてのFunction構造体 に定義されますが、Functionトレイトから外し、個々の構造体にそれぞれ実装しています。その理由は5.6の可変長への拡張のところで説明します。

逆伝播のbackwardメソッドを追加します。順伝播時(つまりcall()の処理中)に記憶しておいた入力x の値をself.input.data として取り出します。このメソッドでは出力側から伝わる微分が渡されます。合成関数の微分の公式より（3.2を参照）、「引数で渡された出力側から伝わる微分の値」（backwardの中のgy を表す）と「その関数のinputの値」（backwardの中のx を表す）を使って計算したその関数での微分の値(squareの場合、導関数は2xなので、xにinputの値を代入して出た値)を掛け算してその値をf32型として返していきます。

バックプロパゲーションの実装

実際に微分をしてみましょう。

fn main() {
    let mut x = Variable::new(0.5);
    println!("{:?}", x);

    let mut A = Square::new();
    let mut B = Exp::new();
    let mut C = Square::new();

    let mut a = A.call(&x);
    let mut b = B.call(&a);
    let mut y = C.call(&b);

    y.grad = Some(1.0);
    b.grad = Some(C.backward(y.grad.unwrap()));
    a.grad = Some(B.backward(b.grad.unwrap()));
    x.grad = Some(A.backward(a.grad.unwrap()));
    println!("{}", x.grad.unwrap());
}
graph LR
 A(("$$x$$")) --> B["$$A$$"]
 B --> C(("$$a$$"))
 C --> D["$$B$$"]
 D --> E(("$$b$$"))
 E --> F["$$C$$"]
 F --> G(("$$y$$"))


graph RL
 A(("$$dy/dy$$")) --> B["$$C'(b)$$"]
 B --> C(("$$dy/db$$"))
 C --> D["$$B'(a)$$"]
 D --> E(("$$dy/da$$"))
 E --> F["$$A'(x)$$"]
 F --> G(("$$dy/dx$$"))
 A ---|"outputのgradを渡す"| B
 B ---|"inputのデータを渡す"| C
 C ---|"次の関数へ"| D

先ほどの「バックプロパゲーションの理論」の図と同じものを用意しました。前半のコードは図の上の順伝播の計算、いわゆる普通のy = f(x)でxに値を代入してyを求める処理をしています。後半のコードは図の下の逆伝播（バックプロパゲーション）で、yをxで微分した値を求めています。図の上と下では赤い線で結ばれていますが、これは例えばXとdy/dxなら、Xはdy/dxをX.gradとして保持しているということです。変数のgradの更新が行われていますが、それはそれぞれの関数で微分の値を次々とxの方へ流しているイメージです。例えばy.gradはdy/dy を、b.grad はdy/db を表しています。そしてdy/dbは関数C のbackward (数学で言うなら導関数)にdy/dyを渡すことで求まります。これを次はdy/da,dy/dxというように繰り返すこのような計算を行うことでdy/dxが求まるのです。これをより一般的に説明するならば、

  1. ある変数を生み出した関数を呼び出し、
  2. その関数のinput変数のdata、output変数のgradを用いてその関数の導関数であるbackwardによってその関数の微分の値を求め、
  3. その値をinput変数のgradに変更する。
  4. そしてinput変数を生み出した新たな関数で同様のことを行う。

この処理を自動化することがこの章の一番の目的です。以上が微分の手作業による実装です。次はこの処理を実際に自動化していきます。

微分の実装（自動化）

前の段階で微分を実装することができました。しかし逆伝播の計算を自分自身で書く必要があります。今回は関数が3つの単純な関数でしたが、これが100、200個と、長い計算グラフを考えたときx.grad =・・・ という逆伝播のコードを数百個すべて手作業で書かなくてはなりません。この章では順伝播の関数に対する逆伝播が自動的に行われる仕組みを作ります。

微分自動化のための変更点

微分を自動化するためには、変数と関数の「つながり」について考えなければなりません。手動で微分したコードでは人が自分たちでgradの値を変更していたため、つながりを考えなくても正しく微分できました。しかし、自動化するにあたって、このつながりを正しくさかのぼる必要があります。正しく遡るにあたって、変数と関数がどのようにつながっているか間違えることのないよう正しく保存しておかなくてはなりません。

関数の目線から変数がどのようにみえるかというと、変数は「入力される変数」と「出力される変数」の２種類存在します。

graph LR
 A(("$$x$$")) --> B["$$f(x)$$"]
 B --> C(("$$y$$"))
 A ---|"input"| B
 B ---|"output"| C

続いて、逆に変数の目線から関数がどう見えるのか考えてみましょう。ここで注目すべき点は「関数の作成」の章で言ったように変数は関数によって作り出されるということです。関数は変数を入力として関数に渡し、出力として新たな変数を生み出します。言い換えると変数にとって関数は「creator（生みの親）」です。

graph LR
 A(("$$x$$")) --> B["$$f(x)$$"]
 B --> C(("$$y$$"))
 B ---|"creator"| C

ではその関数と変数の「つながり」を私たちのコードに取り入れましょう。

struct Variable {
    data: f32,
    grad: Option<f32>,
    creator: Option<Rc<RefCell<dyn Functions>>>,
    name: Option<String>,
}

impl Variable {
    fn new(data: f32) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Variable {
            data,
            grad: None,
            creator: None,
            name: None,
        }))
    }
    fn set_creator(&mut self, func: &Rc<RefCell<dyn Functions>>) {
        self.creator = Some(Rc::clone(func));
    }
}

このコードを説明する前に Rc<Refcell>型 について説明します。 通常Rustの所有権の考え方では、データは一度に一つの変数しか所有できません。今までは一変数の関数のみを扱ってきたので、所有権の扱いは簡単でした。しかし、これからいろいろな関数(多変数関数)を実装していく中で、複雑な関数も実装するためRustの所有権を管理するのが大変になってきます。なので共同保有という考え方をVariable構造体に導入します。

graph LR
 A[A] --> a((a))
 a --> B[B]
 B --> b((b))
 a --> C[C]
 C --> c((c))

上の計算グラフの場合を考えてみましょう。この関数では関数Aを用いて変数aを作り出し、関数Bが変数aを参照し変数bを、関数Cが変数aを参照し変数cを作り出しています。ここで重要なのは関数Bと関数Cがどちらとも変数aという１つの変数を参照していることです。

これを解決するためにRustは Clone() というトレイトを実装しています。このトレイトは新しいメモリを確保し、完全なコピーをします。しかし私たちが作るフレームワークは複雑な関数を何度も使用するため、毎回コピーすると処理が重く、メモリも莫大に必要になってきます。

Rcは「参照カウント」型で複数の所有者を可能にしますが内部のデータへの不変な参照しか提供しません。RefCellは borrow_() 、 borrow_mut() によって、実行時における可変性（Interior Mutability）を可能にします。つまりRc<Refcell>型は所有権の共有と内部のデータを可変に操作できるというRustでは難しい特徴を両立できる型なのです。まとめると、Rc型は所有権の共同保有を、Refcell型は共同保有されたものを可変に扱うことを可能にしているのです。

Variable構造体にRc<RefCell> を導入すると、Rc<Refcell<Variable>> 構造体となります。これはもとのVariable構造体を可変な共同保有ができるようにしたものです。しかし、Variableの内部のデータにアクセスするにはborrow()を多用しなければなりません。またこの構造体を型で示すとき、毎回、Rc<Refcell<Variable>>と書かなくてはならず、面倒です。そこでこのRc<Refcell<Variable>>構造体を一つの構造体として実装してみましょう。ここではRc型を用いているのでRcVariable型とします。

では可変な参照の共有ができるRc<Refcell>型を用いてRcVariable構造体 を実装してみましょう。

ここで用いるRc、RefCellはRustの中でも扱いが難しい概念です。特に borrow() などは扱い方を知らないと簡単にエラーが起きます。なので事前に調べておくことをお勧めします。これらの参考資料はGitHubのreadmeから見ることもできます。

struct Variable {
    data: f32,
    grad: Option<f32>,
    creator: Option<Rc<RefCell<dyn Function>>>,
    name: Option<String>,
}

impl Variable {
　　fn new_rc(data: f32) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Variable {
            data: data,
            grad: None,
            creator: None,
            name: None,
            id: id_generator(),
        }))
    }

    fn set_creator(&mut self, func: &Rc<RefCell<dyn Function>>) {
        self.creator = Some(Rc::clone(func));
    }
}

#[derive(Debug, Clone)]
pub struct RcVariable(pub Rc<RefCell<Variable>>);

impl RcVariable {
    pub fn new(data: f32) -> Self {
        RcVariable(Variable::new_rc(data.to_owned()))
    }
   
    pub fn data(&self) -> f32 {
        self.0.borrow().data.clone()
    }

    pub fn grad(&self) -> Option<RcVariable> {
        self.0.borrow().grad.clone()
    }
    
}

trait Function{
    fn call(&mut self, input: &RcVariable) -> RcVariable;
    fn forward(&self, x: f32) -> f32; // 引数f32
    fn backward(&self, gy: f32) -> f32; // 引数f32
    fn get_input(&self) -> RcVariable;
    fn get_output(&self) -> RcVariable;
}

  • pub fn new_rc(data: f32) 初期化したRcVariableを呼び出す関数。前のVariableを生成するnew()とは別物なので注意。

  • self.0.borrow_mut()
    可変借用 を取得し、内部のVariableに対してbackwardメソッドを呼び出すことで勾配情報を更新している。

  • fn data(&self)
    RcVariable に含まれるVariableの値(dataフィールド)の取得

  • fn grad(&self) RcVariable に含まれるVariableの微分の値(gradフィールド)の取得

ここで変更点は主に二つ、一つ目はVariableをRc化して共同所有 できるようにすること、そして、もう一つは共同所有ができるようになったVariableをRcVariableとして新たな構造体として定義する ことです。

はじめにVariableの変更点について説明します。

フィールドとしてOption<Rc<RefCell<Functions>>>型を creator として保持し、Option<String>型を name として保持します。また、初期化の関数 new() にもcreatorとnameを追加します。次にcreatorを保持するための関数 set_creater を追加します。この関数はFunctionの共有された所有権をコピーし、Variableのcreatorに格納します。これにより、Variableは、自分を生み出したFunctionをたどることができ、そのFunctionの情報にアクセスできるようになります。

RcVariableはRc<Refcell<Variable>>をタプル構造体として定義します。なので、実際には(Rc<Refcell<Variable>>、)となっていてタプルとしてVariableを保持します。タプルの要素にアクセスする際は taple.0 など、インデックスの値をつけるのでした。 だから self.0.borrow_mut() のとき.0を使う のです。このRcVariable構造体のイメージは饅頭です。中の餡がVariable でそのVariable のdataやgrad にアクセスするのに、self.0 を用います。 なのでRcVariableはあくまでVariableを使いやすくするためのもので、データの保持など、本質的な構造体は中身のVariableなのです。

RcVariableとして一つの構造体を定義することで、borrow() を用いる関数などをこの構造体に実装し、まとめることができます。また、構造体として定義したおかげで、RcVariable特有の関数を簡単に実装でき、よりコードの可読性が増します。

これによって私たちは共同保有の考えを用いて、可変な参照の共有ができるRcVariable構造体を構築することができました。

Functionトレイトの変更点について説明します。

関数としてget_inputとget_outputを追加します。この関数は２つともFunction構造体のインプットとアウトプットのRcVariable型を返すものです。この関数はinputとoutputとFunction構造体のつながりを保つという意味で重要なものとなります。 ここでなんでフィールドにアクセスする関数をわざわざFunctionトレイトに実装しないといけないのか疑問を持たれるかもしれません。これに関しては後の5.3のbackward()メソッドを実装する際に説明します。 ちなみにFunctionトレイトで関数をいくつか実装しましたがここで型を指定する際、RcVariableと何度も書きましたが、RcVariableを導入していないと、全部Rc<RefCell<Variable>>と書かなくてなりません。

各関数の変更点

5.1によってVariableの自動微分への対応はできましたが、まだ各関数には対応していません。なので次はSqure構造体を例にして実装していきましょう。

下のコードは変更した関数のみを表示しています。forward()やbackward()などは省略しています。

struct Square{
    input: Option<RcVariable>,
    output: Option<Weak<RefCell<Variable>>>,
}

impl Function for Square {
       fn call(&mut self, input: &RcVariable) -> RcVariable {
        let x = input.borrow().data;
        let y = self.forward(x);

        let output = Variable::new(y);

        self.input = Rc::clone(input); //inputを保存
       
　　　　 self.output = Rc::downgrade(&output); //outputを保存、弱参照で
       
        let self_f: Rc<RefCell<dyn Function>> =        Rc::new(RefCell::new(self.clone()));
        output.borrow_mut().set_creator(&self_f); //outputに自分を覚えさせる

        output
    }
    
    fn get_input(&self) -> RcVariable {
        Rc::clone(&self.input)
    }

    fn get_output(&self) -> RcVariable {
        Rc::clone(self.output.upgrade().as_ref().unwrap())
    }
}
impl Square {
    fn new() -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Self {
            input: Variable::new(0.0),
            output: Rc::downgrade(&Variable::new(0.0)),
        }))
    }
}

ここでは以前の関数構造体から大きく変更しているので、一つ一つ段階を踏んで説明していきます。なおここではSquare構造体のみ説明していきます。Exp構造体はこれをもとに変更してください。

  • フィールドの変更
    Squareを構造体として定義していますが、そのフィールドを変更していきます。はじめにinputのVariableをRcVariableに変更します。

  • Square::new()の変更点
    Rc<RefCell<Self>>>を返すことで、Square関数自体が共有可能かつ内部可変なオブジェクトとして扱われるようにします。 逆伝播で値を取得するために使用するためVariableへの強い参照（所有権を共有）つまりinputの値を長時間保持しなくてはならないため。self.input = Rc ::clone(input); でcloneを使っています。ここでFunctionがinputとのつながりを保持します。

self.output = Rc::downgrade(&output); はなぜcloneしないのかというと もしFunctionがoutputを強い参照で持ってしまうと、VariableのcreatorフィールドがFunctionを参照し、Functionがoutputを参照するという循環参照 が発生します。

graph LR
 f["$$f(x)$$"] ---|"output"| y(("$$y$$"))
 f --> y
 f ---|"input"| y

この図のように関数と変数が互いに参照を持っています。参照カウントの場合、参照カウントが0になるとオブジェクトは削除されるのですが、この場合、互いに一つのカウントを持っているため、決して消えることがありません。そのため参照カウントを増やさないdowngradeを使っています。

output.borrow_mut().set_creator(...) では出力Variableのcreatorフィールドに、この計算を行ったFunctionインスタンスへの参照を与えることで逆伝播の経路をつくります。

まとめると、Function構造体のこれらの処理は、

  • Function構造体がinputを覚える。
  • outputを生成し、それも覚える。ただし、弱参照で。
  • outputに自分自身が creator であることを覚えさせる。

となっています。これらの処理は順伝播を行う際に呼び出されるFunction構造体の call() によって自動で行われます。 　　

TODO: つながりを確認するコードを追加予定。

これらの設定により、順伝播の計算をするときに同時に、inputとFunctionとoutputのつながりを自動で保存します。これこそが動的に「つながり」を作る核心の部分です。

これを参考にして自分自身の手でExp関数を変更してみましょう。わからないときはGithubリポジトリを参照してください。

backwardメソッドの追加

前のところで変数と関数のつながりを自動で保存するシステムを構築しました。次はこの作られたつながりを自動で辿ることです。

4.3で手動で微分をした図とコードを改めて見てみましょう。

graph LR
 A(("$$x$$")) --> B["$$A$$"]
 B --> C(("$$a$$"))
 C --> D["$$B$$"]
 D --> E(("$$b$$"))
 E --> F["$$C$$"]
 F --> G(("$$y$$"))

graph RL
 A(("$$dy/dy$$")) --> B["$$C'(b)$$"]
 B --> C(("$$dy/db$$"))
 C --> D["$$B'(a)$$"]
 D --> E(("$$dy/da$$"))
 E --> F["$$A'(x)$$"]
 F --> G(("$$dy/dx$$"))

y.grad = Some(1.0);
b.grad = Some(C.backward(y.grad.unwrap()));
a.grad = Some(B.backward(b.grad.unwrap()));
x.grad = Some(A.backward(a.grad.unwrap()));
println!("{}", x.grad.unwrap());

ここは手動で微分したコードの部分です。はじめにy.gradを初期化しています。それはdy/dy=1.0だからです。次にbのgradであるb.gradの値を変更しています。bのgradはいわばdy/dbです。 これをもとめる数式は(dy/dy)・(dy/db)です。dy/dyは先ほどのy.gradの1.0、そして、dy/dbはCのbackwardにinputであるbのdataを渡すことで求まります。試しにy=2*bを考えてみましょう。この時y=C(b)とかけ、C’(b)は2bです。この導関数にbの値を入れることで微分の値が出ます。

  • creatorである関数を取得
  • 関数の入力のRcVariable(またはVariable)を取得
  • 関数のbackwardメソッドを呼ぶ

という処理を繰り返していることがわかります。その都度コードを書かなくてはなりません。 そこで繰り返しの処理を自動化できるようにVariable構造体にbackwardというメソッドを追加します。

impl Variable {
 fn backward(&self) {
        let mut funcs: Vec<Rc<RefCell<dyn Functions>>> =
            vec![Rc::clone(self.creator.as_ref().unwrap())];

    let mut last_variable = true;
        while let Some(f_rc) = funcs.pop() {
            let f_borrowed = f_rc.borrow();
            
            let x= f_borrowed.get_input();
            let y = f_borrowed.get_output();
            let y_grad: RcVariable;

            if last_variable {
                y_grad: f32 = 1.0; //最初はdy/dy = 1より1に設定
                last_variable = false;
            } else {
                y_grad = y.0.borrow().grad.as_ref().unwrap().clone();
            }
            let x_grad = f_borrowed.backward(&y_grad);
            
            x.0.borrow_mut().grad = Some(x_grad); //gradの更新
            
            if x.borrow().creator.is_none() { 
                break;
            }; //xのcreatorがない時はその変数がはじめの変数なので終了

            funcs.push(Rc::clone(x.borrow().creator.as_ref().unwrap()));
        }  // ↑ xのcreatorのFunction構造体を新たにfuncsリストに追加
    }
}

このbackwardメソッドではVariableのcreatorから関数を受け取り、その関数の入力を取り出した後にbackwardメソッドを呼び出すという処理をループを使って実装しています。

let mut funcs: Vec<Rc<RefCell<Functions>>> = vec![Rc::clone(self.creator.as_ref().unwrap())];

最初の要素として、現在のVariable（self）を生み出したFunctionのcreatorを取得します。unwrap()は、このVariableが既に何らかの計算結果であるという前提（終端ではない）なので使えます。
let x = f.borrow().get_input(); では取り出したFunctionの入力変数（Variable）を取得します。

今回はlast_variableによって最初に処理するVariableか途中のVariableかで場合分けします。それは最初のoutputにあたる変数はこのbackward()を実行しているVariableであり、すでにselfが借用されているためです。つまり最初のVariableでborrow()を使ってgradを取り出すことができないことを指しています。最初に処理するVariableは必ずdy/dyなので１.0です。一方途中のVariableは以前のFunctionから勾配を取得する必要があります。y_grad: Functionの出力側から伝わってきた微分の値です。

f.borrow().backward(y_grad) ではFunction のbackwardメソッドに出力の微分の値を渡し、微分の値を計算します。

x.borrow_mut().grad = Some(...) では 計算された微分の値を入力側のVariable のgradフィールドに書き込みます。borrow_mut()が使われているのは、内部のgradフィールドを変更するためです。

if x.borrow().creator.is_none() がtrueの場合、Variable の親のFunctionが存在しないユーザーが作成したVariableであり、これ以上遡る必要のないノードです。つまり微分がすべて終わったことを意味します。ここでループを終了する必要があるのです。

funcs.push(...) では、次のFunctionをfuncsに追加します。これにより、次のループでそのFunctionを処理することができます。

forward,backwardのRcVariable対応

可変長への拡張

今までの関数(exp,square)は一変数関数で入力と出力の個数が一個です。しかし、今後実装するaddやmulといった演算子の関数はa+b,a*bというように入力は2個です。なのでこれからは二変数関数、ひいてはそれ以上の多変数に対応できるように拡張していきます。具体的にはFunctionトレイト・構造体のフィールド、そしてVariableのbackward()メソッドの二つを変更していきます。

trait Function: Debug {
    fn call(&mut self,input: &RcVariable) -> RcVariable;
    fn forward(&self, x: &[RcVariable]) -> RcVariable;
    fn backward(&self, gy: &RcVariable) -> Vec<RcVariable>;
    fn get_inputs(&self) -> [Option<RcVariable>; 2];
    fn get_outputs(&self) -> [Option<RcVariable>; 2];
}

はじめにFunctionトレイトを変更します。注目すべきところはトレイトで定義された関数の引数と戻り値の型です。今までは戻り値の型をf32、RcVariableだけでしたが、ここではスライス型として渡し、出力します。これにより可変個の入力を渡せるようになりました。ここでは配列をvec型としましたが、ここにはRustの動的、静的といった考えがあります。

Rustの複数のデータを保持する型の種類として主に２種類あります。それはVec型と配列型です。これらの違いは保持するデータの個数、すなわち長さがいつ決まるかということです。Vec型ではコードの実行中に、配列ではコンパイル時に決まります。いうなればVecは可変の長さであり、実行中に自由に変更することができます。一方配列はすでに決まった長さ、つまり不変の長さです。一見するとVecの方が自由度が高くて便利だと思いますが、長さが実行中に決まるので、決定するのに多少時間がかかってしまいます。要するにVecは自由度が高い代わりに、パフォーマンスを犠牲にしているということです。逆に配列は自由度が低い分、パフォーマンス的には良いです。まとめると、自由度とパフォーマンス、どっちをとるかということです。ただし実際はこれら二つのパフォーマンスに差が出るのは要素が数万個ぐらいのものを処理する場合であり、数個の場合はまったくと言っていいほど差はありません。よってここでは柔軟性が高いVec型を採用しています。このような静的と動的の考えはRustでは非常に重要になってきます。

add関数の実装

可変長に対応することができたので、実際に2変数関数であるadd関数を実装してみましょう。

#[derive(Debug, Clone)]
struct AddF {
    inputs: Vec<Rc<RefCell<Variable>>>,
    output: Option<Weak<RefCell<Variable>>>,
    id: usize,
}

impl Function for AddF {
    fn call(&mut self) -> RcVariable {
        let inputs = &self.inputs;
        if inputs.len() != 2 {
            panic!("Addは二変数関数です。inputsの個数が二つではありません。")
        }

        let output = self.forward(inputs);

        if get_grad_status() == true {
            
            //  outputを弱参照(downgrade)で覚える
            self.output = Some(output.downgrade());

            let self_f: Rc<RefCell<dyn Function>> = Rc::new(RefCell::new(self.clone()));

            //outputsに自分をcreatorとして覚えさせる
            output.0.borrow_mut().set_creator(self_f.clone());
        }

        output
    }

    fn forward(&self, xs: &[RcVariable]) -> f32 {
        let x0 = &xs[0];
        let x1 = &xs[1];
        let y_data = x0.data() + x1.data();

        y_data
    }

    fn backward(&self, gy: &RcVariable) -> Vec<f32> {
        let gx0 = gy.clone();
        let gx1 = gy.clone();
        
        let gxs = vec![gx0, gx1];
        gxs
    }

    fn get_inputs(&self) -> &[RcVariable] {
        &self.inputs
    }

    fn get_output(&self) -> RcVariable {
        let output;
        output = self
            .output
            .as_ref()
            .unwrap()
            .upgrade()
            .as_ref()
            .unwrap()
            .clone();

        RcVariable(output)
    }

    fn get_id(&self) -> usize {
        self.id
    }
}
impl AddF {
    fn new(inputs: &[RcVariable]) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Self {
            inputs: inputs.to_vec(),
            output: None,
            id: id_generator(),
        }))
    }
}

pub fn add(xs: &[RcVariable]) -> RcVariable {
    AddF::new(xs).borrow_mut().call()
}

Add関数の微分ですが、ここで初めて一つの変数の一変数関数ではなく、複数の変数の多変数関数が登場しました。今までは一変数でつながりも一直線でしたが、多変数関数の場合は分岐します。ではどのように微分を求めるのでしょうか。ここで必要な知識は偏微分です。

graph RL
 z(("$$z$$")) -->|"$$\frac{\partial L}{\partial z}$$"| mul["$$Mul(x,y)$$"]
 mul --> z
 x --> mul
 y(("$$y$$")) --> mul
 mul -->|"$$\frac{\partial L}{\partial z}\cdot 1$$"| y
 mul -->|"$$\frac{\partial L}{\partial z}\cdot 1$$"| x(("$$x$$"))

ここではz(x,y)の関数を考えてみましょう。\(z=x+y\)とすると、\(z\)は独立する\(x\)と\(y\)によって定まるので二変数関数です。これは足し算の関数ですが、\(x、y\)のそれぞれの\(z\)に対する偏微分を考えます。すると\(\frac{\partial z}{\partial x} = 1\)、\(\frac{\partial z}{\partial y} = 1\)となります。よって上から流れて来た微分の値、図で言うなら∂L/∂z、コードで言うならgyをそのまま流すことになります。このように偏微分を駆使して実装すれば多変数関数でも正しく逆伝播を行うことができます。

Add構造体の名前がAddFなのは、今後これを演算子のオーバーロードで用いるのですが、すでにAdd構造体がもともと存在するので、それと名前が重複するのを防ぐためです。

同じ変数を繰り返し使う

TODO: 今後書く予定

微分の理論（複雑な関数）

これまでわたしたちは下のグラフのような一直線に並ぶ計算グラフを扱ってきました。

graph LR
 A((X)) --> B[A]
 B --> C(( a ))
 C --> D[B]
 D --> E((b))
 E --> F[C]
 F --> G((y))

しかし変数と関数の「つながり」はそのような一直線とは限りません。これまでの実装により私たちはadd関数のように２変数への拡張を行ってきました。それによってより複雑な「つながり」を作ることができました。しかし今のフレームワークでは複雑な「つながり」の逆伝播を正しくすることができません。 今のフレームワークにどのような問題があるのか調べるために１つのシンプルな計算グラフについて考えましょう。

graph LR
 X((X)) --> A[A]
 A --> a((a))
 a --> B[B]
 B --> b((b))
 a --> C[C]
 C --> c((c))
 b --> D[D]
 c --> D
 D --> y((y))

この計算グラフで注目したい点は変数aです. 前の章で紹介したようにaの微分にはaの出力側から伝播する２つの微分が必要になります。その点を注力してまずはじめに正しい計算グラフの流れを一つ一つ表示してみます。（矢印は逆伝播を表しています。
※ここではb,cへの逆伝播の順番は問題ではないので、そこの流れは省略しています。

1

graph RL
 y((y)) --> D[D]
 D --> b((b))
 D --> c((c))

2

graph RL
 y((y)) --> D[D]
 D --> b((b))
 D --> c((c))
 c --> C[C]
 C --> a((a))

3

graph RL
 y((y)) --> D[D]
 D --> b((b))
 D --> c((c))
 b --> B[B]
 c --> C[C]
 B --> a((a))
 C --> a

4

graph RL
 y((y)) --> D[D]
 D --> b((b))
 D --> c((c))
 b --> B[B]
 c --> C[C]
 B --> a((a))
 C --> a
 a --> A[A]
 A --> X((X))

今示した順番で伝播していますが関数の目線からいうとD、B、C、Aの順番で逆伝播を行っています。ここでB,Cの順番はどちらでもかまわないのですが、大切なのは関数B、Cの逆伝播を終わらせてから、関数Aの逆伝播をおこなうということです。

では私たちのプログラムは実際にこのように正しく関数を取り出してくれるのでしょうか。下のコードで実行してみましょう。
TODO: コード追加予定

すると私たちの計算グラフは次の順番で処理します。

1

graph RL
 y((y)) --> D[D]
 D --> b((b))
 D --> c((c))

2

graph RL
 y((y)) --> D[D]
 D --> b((b))
 D --> c((c))
 c --> C[C]
 C --> a((a))

3

graph RL
 y((y)) --> D[D]
 D --> b((b))
 D --> c((c))
 c --> C[C]
 C --> a((a))
 a --> A[A]
 A --> X((X))

4

graph RL
 y((y)) --> D[D]
 D --> b((b))
 D --> c((c))
 b --> B[B]
 c --> C[C]
 B --> a((a))
 C --> a
 a --> A[A]
 A --> X((X))

5

graph RL
 y((y)) --> D[D]
 D --> b((b))
 D --> c((c))
 b --> B[B]
 c --> C[C]
 B --> a((a))
 C --> a
 a --> A[A]
 a --> A
 A --> X((X))
 A --> X

ではここで、この処理がどのように行われるのか考えてみましょう。重要なのは、Variableのメソッドのbackwardの中のfuncsです。逆伝播で遡るにあたって、一つ前の関数をこのリスト(正しくはベクタ)に入れ、そこから関数を取り出してbackwardを行い、その関数のinputのcreatorである関数を再び追加するという流れですが、このような多変数関数になると、リストの中には複数の関数が含まれる状態になります。なので、いかに正しい順番で関数を取り出せるかが重要となってきます。逆に言えば、今までの一変数関数の場合は、リストに常に関数は一つしかないので、取りだし方を考えなくてもよかったのです。

現状の設計ではこのように、本来B、Cの関数を先に呼び出した後、Aを呼び出すべきであるにもかかわらず、実際はCの後、先にAをとりだしてしまっています。さらにBを取り出した後、Bのinputであるaが認識されることで、再びAが取り出され、backwardするという処理を行ってしまいます。ではこの場合のfunsベクタの要素の挙動を先ほどの逆伝播の流れとともに考えてみましょう。

vecの追加、取り出しは、どちらとも後ろから追加、取り出しです。

  1. 最初にfuncsは初期化されるので、funcs = []　つまり空です。

  2. はじめにyのcreatorである関数Dが追加されます。funcs = [D]

  3. 次にDが取り出され、DのinputであるbのcreatorのBが追加されます。funcs = [B]

  4. Dのもう一つのinputであるcのcreatorのCが追加されます。funcs = [B,C]

  5. Cが取り出され、CのinputであるaのcreatorのAが追加されます。funcs = [B,A]

  6. Aが取り出され、Aのinputであるxのcreatorはないので何も追加されません。funcs = [B]

  7. Bが取り出され、BのinputであるaのcreatorのAが追加されます。funcs = [A]

  8. Aが取り出され、Aのinputであるxのcreatorはないので何も追加されません。funcsは空になったので終了。funcs = []

今までのbackwardメソッドでは計算グラフが直線的だったため、何も考えずに前の関数を呼んでいましたが、今後は適切な順番に関数を取り出すことが求められます。ではどのように正しい順番で取り出すのか。まずは正しい順番とは何かを考えて、そこから一般化してみましょう。

先ほどの例に出した関数で考えてみましょう。正しい順番としては、D→(C→B)→Aでした。ここでCとBの間で括弧をしたのは、どちらの順番でも良いという意味を示すためです。ではなぜBとCは順番がどちらでもよいのでしょうか。それは、同じ変数をinputとしているからです。aに向かう微分、(ここでは∂b/∂a,∂c/∂aを指す)が求まってはじめて、変数aのVariableとしての微分gradが定まるからです。これが定まらないと、xの微分が求まりません。なぜならxの微分を求めるAのbackwardはaの微分を引数としているおり、aの微分に依存しているからです。つまり、同じinputを持っている関数同士では、順位は同じであり、同じ順位の関数が全て呼び出されてはじめて次の関数につなげられるということです。このことから、関数の取り出す優先順位は、inputの順番で決まると考えられます。

次にinput、すなわち変数となるVariableの順番を考えましょう。これは源流であるxに近い方が優先順位が低いと言えます。なぜなら、源流の真反対である終端のyから源流にたどっていくからです。今回はx←a←(b←c)←yなので、優先順位は左に行くほど低く、右に行くほど高いです。このVariableの優先順位は、変数が生まれた順番で定義することができます。 例えばxが一番はじめに生まれるので0とし、次に生まれるaを1、その要領で数字を割り振っていけば、数字が高いほど優先順位が高いと定義できます。

以上のことをまとめると、関数の優先順位はinputの順番で決まり、そして、そのinputの順番はそのVariableが生まれた順番で決まるということが分かりました。あとはVariableの生まれた順番をどう求めるかです。それは順伝播でつながりを構築している最中で求めることができます。ある関数fが変数xをinputとして受け取ったとき、xの優先順位を読み取り、outputとなる変数を生み出す際にその優先順位より1優先順位が高い値を持たせればいいのです。ではこれらの処理を実際にプログラムに取り込むために、改めて一から処理を考えてみましょう。

まず優先順位の基準となる値を世代としましょう。この値は自然数です。はじめの変数(ここではx)には0を持たせます。この変数を関数に渡す際、関数の世代をinputの変数の世代と同じに設定します。そうすれば、関数の優先順位はinputの順番で決まるというルールを従うことができます。そして、関数がoutputを出力する際、その変数の世代の値を自身の世代に1足して持たせれば、新たな世代を生み出すことができます。言葉だけの説明ではイメージが湧かないので、世代というものがどのように与えられるか可視化します。

forward

graph LR
 subgraph 世代:0
 X((X)) --> A[A]
 end
 subgraph 世代:1
 A --> a((a))
 a --> B[B]
 a --> C[C]
 end
 subgraph 世代:2
 B --> b((b))
 C --> c((c))
 b --> D[D]
 c --> D
 end
 D --> y((y))

backward

graph RL
 y((y)) --> D[D]
 subgraph 世代:2
 D --> b((b))
 D --> c((c))
 end
 b --> B[B]
 c --> C[C]
 subgraph 世代:1
 B --> a((a))
 C --> a
 end
 subgraph 世代:0
 a --> A[A]
 A --> X((X))
 end

forwardの場合、例えばxが世代1と設定されているので、Aも世代0と設定されます。そして、Aのoutputであるaは世代が0に1足されて1となります。これを繰り返せば、グラフのような世代の関係が生まれます。 その後逆伝播では、順伝播で作られた世代の関係の大きい順に関数を取り出せば、正しい順番で処理できます。先ほどの問題では、Cの後にBではなくAを取り出してしまいましたが、世代の順番に従えば、世代が0のAより、世代が1であるBが先に取り出されるようになります。

ではこれらの理論をもとに、次の章で実際に実装していきましょう。

微分の実装（複雑な関数）

idの設定

世代を持たせる前に、idというものを設定します。さらに関数が複雑になり、変数(Variable)と関数(Function)が増えてくると、管理が大変なため、今後のためにそれぞれの構造体にidをつけます。idを個々の構造体に付与することで、後の複雑な関数の自動微分でバグが起きないように安全に処理することができます。

idの生成

idを生成する NEXT_ID というグローバル変数を設定します。この変数はどこのプログラムからでもアクセスできる変数です。扱いやすいですが、同時にアクセスし、変更することができるので、安全ではありません。そこで使われるのがAtomic型です。複数のスレッドで同時に使用される変数の場合に使われる型です。AtomicUsizeはその中のusize型を扱うものです。これを用いてグローバル変数を設定します。この変数の現在の値に1を加算し、生成されたVariableや関数のidとして渡します。そのようにして構造体が作成されるたびに新しいidが付与されます。イメージとしては整理券です。NEXT_IDは整理券を発行していて、券を取るたびに番号が1ずつ増えていくのです。このいわば券を取り、番号を1足す作業を id_generator() として関数にします。この関数を呼び出せば、その時点のidが返されます。次のVariableやFunction構造体の変更のところで用います。

use std::sync::atomic::{AtomicBool, AtomicUsize, Ordering};
static NEXT_ID: AtomicUsize = AtomicUsize::new(0);
/// idを生成する関数。構造体のコンストラクタを作成する際に、呼び出して、idを付ける
fn id_generator() -> usize {
    NEXT_ID.fetch_add(1, Ordering::SeqCst)
}

//例
fn main() {
    let first = id_generator();
    println!("first = {}",first);
    let second = id_generator();
    println!("second = {}",second);
}

Variableの変更

Variableに id というフィールドを持たせましょう。このidの値をusize型として保持します。

struct Variable {
    data: f32,
    grad: Option<f32>,
    creator: Option<Rc<RefCell<dyn Function>>>,
    name: Option<String>,
    id: usize,
}

impl Variable {
    fn new_rc(data: ArrayD<f32>) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Variable {
            data: data,
            grad: None,
            creator: None,
            name: None,
            id: id_generator(),
        }))
    }
}

この時、idは先ほど作成した id_generator() 関数を用いることで、正しいidが付与されます。このidは呼び出されるたびに新たなidを返すので、重複することはありません。

Functionトレイトおよび構造体の変更

はじめにFunctionトレイトを変更します。具体的にはFunction構造体のidを返す関数を追加します。

trait Function: Debug {
    fn call(&mut self) -> RcVariable;

    
    fn forward(&self, x: &[RcVariable]) -> RcVariable;
    fn backward(&self, gy: &RcVariable) -> Vec<RcVariable>;

    //　関数クラス.inputs, .outputではvariableのbackwardの中でアクセスできないので、関数にして取得
    fn get_inputs(&self) -> &[RcVariable];
    fn get_output(&self) -> RcVariable;
    fn get_id(&self) -> usize;  //  <--  今回追加するもの 
}

トレイトにidを返す関数を追加する理由は、5.3で説明したget_input 、get_output と同じです。idフィールドもinputsフィールド、outputフィールドと同様に保持しているか不明なので、トレイト内で関数として定義しているということです。これは後のgeneration でも同じことなので、同じようにget_generation を追加することになります。

次にFunction構造体の変更です。Exp構造体を例にして変更しますので、その他のものはそれに従って変更してください。

struct Exp {
    inputs: Vec<RcVariable>,
    output: Option<Weak<RefCell<Variable>>>,
    id:usize,
}

impl Function for Exp {
    fn get_id(&self) -> usize {
        self.id
    }
}

先ほどのget_id()をFunctionトレイトして実装します。

以上により、VariableとFunctionトレイト・構造体の id への対応ができました。

世代（ジェネレーション）の保持

では前節で説明した世代を持たせ、それに従って処理するよう変更していきます。

Variableの変更

まずVariableに generation というフィールドを持たせましょう。このgenerationには世代の値を保持します。例えば前のグラフのXは第0世代なので0という値を持ちます。

struct Variable {
    data: f32,
    grad: Option<f32>,
    creator: Option<Rc<RefCell<dyn Function>>>,
    name: Option<String>,
    id: usize,
    generation: i32,

}

impl Variable {
    pub fn new_rc(data: ArrayD<f32>) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Variable {
            data: data,
            grad: None,
            creator: None,
            name: None,
            generation: 0,
            id: id_generator(),
        }))
    }
}

generationの値は0以上の整数のみを扱うので、i32型に設定します。はじめは0として設定し、次に初期化したgenerationを正しいgenerationに変更する処理を加えます。

impl Variable {
    pub fn new_rc(data: ArrayD<f32>) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Variable {
            data: data,
            grad: None,
            creator: None,
            name: None,
            generation: 0,
            id: id_generator(),
        }))
    }

    fn set_creator(&mut self, func: Rc<RefCell<dyn Function>>) {
        self.creator = Some(Rc::clone(&func));
        self.generation = func.borrow().get_generation() + 1;
    }
}

Variableのメソッドとして set_creator がありますが、これを変更します。Variableが自分のcreatorを覚える関数ですが、その時に自分の世代を、creatorの世代に1足して設定します。この作業は関数がoutputを出力するとき、自身の世代に1足してoutputに持たせる作業を指します。

Function構造体の変更

次にFunction構造体の変更です。Exp構造体を例にして変更しますので、その他のものはそれに従って変更してください。

struct Exp {
    inputs: Vec<RcVariable>,
    output: Option<Weak<RefCell<Variable>>>,
    generation: i32,
    id:usize,
}

impl Function for Exp {
    fn call(&mut self) -> RcVariable {
        let inputs = &self.inputs;
        if inputs.len() != 1 {
            panic!("Expは一変数関数です。inputsの個数が一つではありません。")
        }

        let output = self.forward(inputs);

        //inputのgenerationで一番大きい値をFuncitonのgenerationとする
        self.generation = inputs.iter().map(|input| input.generation()).max().unwrap();

        //  outputを弱参照(downgrade)で覚える
        self.output = Some(output.downgrade());

        let self_f: Rc<RefCell<dyn Function>> = Rc::new(RefCell::new(self.clone()));

        //outputsに自分をcreatorとして覚えさせる
        output.0.borrow_mut().set_creator(self_f.clone()); //先ほどset_creator()を変更したので、Variableの世代はここ関数によって変更される。
        output
    }

    fn get_generation(&self) -> i32 {
        self.generation
    }
}

はじめにVariableと同じく、構造体にgenerationフィールドを持たせましょう。初期値もVariableと同じく0です。

次にcallの中で、世代をinputのVariableの世代に設定します。多変数関数の場合、inputが複数存在するので、inputの世代で最も値が大きいものを採用します。

また、自身の世代を取り出す関数(get_generation)を作成します。この関数は　FunctionTrait で定義しましたが、その理由は後の 関数を取り出す処理の変更 のところで説明します。

以上により、VariableとFunction構造体のgenerationへの対応ができました。

関数を取り出す処理の変更

前節の説明の通りに、Variableのbackwardの funcsベクタ から関数を取り出す処理を変更します。

impl Variable {
    fn backward(&self, double_grad: bool) {
        let mut funcs: Vec<Rc<RefCell<dyn Function>>> =
            vec![Rc::clone(self.creator.as_ref().unwrap())];

        let mut seen_set = HashSet::new();

        fn add_func(
            funcs_list: &mut Vec<Rc<RefCell<dyn Function>>>,
            seen_set: &mut HashSet<usize>,
            f: Rc<RefCell<dyn Function>>,
        ) {
            if seen_set.insert(f.borrow().get_id()) {
                funcs_list.push(Rc::clone(&f));
                funcs_list.sort_by(|a, b| {
                    a.borrow()
                        .get_generation()
                        .cmp(&b.borrow().get_generation())
                });
            }
        }

        //&selfで最初の変数はborrowされるので場合分け
        let mut last_variable = true;
        while let Some(f_rc) = funcs.pop() {
            //println!("f = {:?}\n", get_struct_name(&f_rc.borrow()));
            let f_borrowed = f_rc.borrow();
            let xs = f_borrowed.get_inputs();
            let y = f_borrowed.get_output();

            let y_grad: RcVariable;

            if last_variable {
                y_grad = ArrayD::<f32>::ones(self.data.shape()).rv();

                last_variable = false;
            } else {
                //関数の出力は一つだけなので、[1]は必要なし
                y_grad = y.0.borrow().grad.as_ref().unwrap().clone();
            }

            let xs_grad = f_borrowed.backward(&y_grad);

            for (x, x_grad) in zip(xs, &xs_grad) {
                // gradをすでに保持しているなら、元のgradに新たなgradを足す。
                // gradをまだ持っていないならそれを持たせる。
                if let Some(current_grad_data) = x.grad() {
                    x.0.borrow_mut().grad = Some(current_grad_data + x_grad.clone());
                } else {
                    x.0.borrow_mut().grad = Some(x_grad.clone());
                }

                // creatorがあるならその関数をfuncsに追加
                if let Some(func_creator) = &x.0.borrow().creator {
                    add_func(&mut funcs, &mut seen_set, func_creator.clone());
                }
            }
        }
    }
}

backwardの中に add_func() という関数を追加します。この関数がまさに、funcsベクタ から関数を正しく取り出すための処理です。この add_funcs 関数が行っていることは主に二つです。一つ目は 関数の重複確認 です。

はじめに backward 内で funcs とは別に新たな配列 seen_set というものを用意します。これは HashSet という型です。HashSet とは普通の配列とは異なり、値の重複を検知することができます。これは funcs に関数が追加される際、今までに追加された関数と同じもの、すなわち重複したものが間違えて入っていないかを確かめるためのものです。前節の 微分の理論 を思い出してください。間違えた処理の場合、関数Aを2回取り出してしまいました。そこで funcs に追加された関数のidを記憶しておくことで、後に追加されるものが今までのものと重複しないかを確認し、もし重複していたら funcs に追加しないという処理を行えばよいのです。このような処理を加えることで、バグの温床を減らす設計にすることができます。 if seen_set.insert(f.borrow().get_id()) のところでidを用いて重複しているかどうか確認します。もし重複していなかったら、seen_set にidを追加し、funcs に関数を追加します。もし重複していたら、今説明した処理は行われません。

二つ目は世代順への並び替えです。 funcs から関数が取り出されるときは、一番後ろから取り出されるので、世代の小さい順に並び変えてあげれば、世代の小さい方が先に取り出されるということを防げます。これを、Vec型 で提供される sort_by() を用いて並べ替えます。

HashSet や、sort_by についてはgithubリポジトリのREFERENCES.mdの文献をご覧ください。

これらの設定により、 add_funcs を用いて funcs に関数を追加していけば、正しく関数を取り出せるようになりました。なので、あとはコードの最後の //creatorがあるならその関数をfuncsに追加 のところをadd_funcsに変更すればよいだけです。

TODO: 並び替えや取り出す処理を簡単に説明するコード追加予定
TODO: 複雑な関数がうまくバックプロパゲーションしてくれるか確かめるコード追加予定

利便性を高める拡張

前のステップでは実数による自動微分を実装しました。しかしまだ演算子に対応しておらず計算をするための関数をいちいち書かなくてはいけません。このステップの目標としては＋や * などの演算子に対応することです。例をあげるなら、aとbをVariable構造体として、a,bを掛けるとき、**mul(a,b)**と書かなければなりません。a*b と書けると便利になります。これから＋や*など演算子が扱えるようにVariable構造体を拡張していきます。

四則演算を行う関数の作成

まず演算子を使えるようにVariable構造体を拡張する前に私たちのフレームワークに掛け算や割り算をする関数を実装する必要があります。足し算をするadd関数はすでに実装したので、これを参考に引き算のsub関数、掛け算のmul関数、割り算のdiv関数を実装してみましょう。 基本的に+を-、*、/に変更するだけですが、backwardのところは偏微分が関わっているので、微分の計算式を載せます。5.3のadd関数のバックプロパゲーションのグラフを参考に考えてみましょう。答えはgithubのコードで確認しましょう。

  • 足し算（add）
    forward: \(z = x_0 + x_1\)
    backward: \(\partial z/\partial x_0 = 1、 \partial z/\partial x_1 = 1\)
    \(x_0,x_1\)どちらとも偏微分の値は1なので、上流からきた微分の値に1をかける,つまりそのまま流せばよい

  • 掛け算（mul）
    forward: \(z = x_0 \cdot x_1\)
    backward: \(\partial z/\partial x_0 = x_1、 \partial z/\partial x_1 = x_0\)
    \(x_0\)の偏微分は\(x_1\)なので、上流からきた微分の値に\(x_1\)をかけて流す。
    \(x_1\)の偏微分は\(x_0\)なので、上流からきた微分の値に\(x_0\)をかけて流す。

  • 引き算（sub）
    forward: \(z = x_0 - x_1\)
    backward: \(\partial z/\partial x_0 = 1、 \partial z/\partial x_1 = -1\)
    \(x_0\)の偏微分は1なので、上流からきた微分の値をそのまま流す。
    \(x_1\)の偏微分は-1なので、上流からきた微分の値をマイナスにして流す。

  • 割り算（div）
    forward: \(z = x_0 / x_1\)
    backward: \(\partial z/\partial x_0 = 1/x_1、 \partial z/\partial x_1 = -x_0/x_1^2\)
    \(x_0\)の偏微分\(1/x_1\)なので、上流からきた微分の値に\(1/x_1\)をかけて流す。
    \(x_1\)の偏微分は\(-x_0/x_1^2\)なので、上流からきた微分の値に\(-x_0/x_1^2\)をかけて流す。

  • 負数（neg） ・・・この関数は一変数関数なので、普通の微分です。
    forward: \(z = -x_0\)
    backward: \(\partial z/\partial x_0 = -1\)
    \(x\)の微分は-1なので、上流からきた微分の値をマイナスにして流す。

また、四則演算と呼べるのかは定かではありませんが、ここで定義しておくと便利なので累乗を計算する関数 pow も実装しましょう。pow関数は \(z = x^c\) で表され、cは定数とします。すると、微分は、\(\partial z/\partial x = cx^{c-1}\)となります。

表にまとめると、このようになります。

関数ForwardBackward (\(\partial z/\partial x_0\))Backward (\(\partial z/\partial x_1\))
add\(z = x_0 + x_1\)\(\partial z/\partial x_0 = 1\)\(\partial z/\partial x_1 = 1\)
mul\(z = x_0 \cdot x_1\)\(\partial z/\partial x_0 = x_1\)\(\partial z/\partial x_1 = x_0\)
sub\(z = x_0 - x_1\)\(\partial z/\partial x_0 = 1\)\(\partial z/\partial x_1 = -1\)
div\(z = x_0 / x_1\)\(\partial z/\partial x_0 = 1/x_1\)\(\partial z/\partial x_1 = -x_0/x_1^2\)
neg\(z = -x_0\)\(\partial z/\partial x_0 = -1\)なし
pow\(z = x^c\)\(\partial z/\partial x = cx^{c-1}\)なし

TODO: 四則演算関数の微分が正しく行われるか試すコードを追加予定

演算子のオーバーロード

この章ではいままで実装してきた演算子の関数をオーバーロードしていきます。Add関数とNeg関数を例にして説明します。

use std::ops::{Add, Div, Mul, Neg, Sub};

impl Add for RcVariable {
    type Output = RcVariable;
    fn add(self, rhs: RcVariable) -> Self::Output {
        // add_op関数はRc<RefCell<Variable>>を扱う
        let add_y = add(&[self.clone(), rhs.clone()]);
        add_y
    }
}

impl Neg for RcVariable {
    type Output = RcVariable;
    fn neg(self) -> Self::Output {
        let neg_y = neg(&[self.clone()]);
        neg_y
    }
}

そもそもオーバーロードとは何でしょうか？実際のところ、オーバーロードは関数オーバーロードというものを指すことが多く、pythonやC系言語などに用いられるものであり、Rustにはそのような関数を再定義する機能は一般的ではありません。その代わり、トレイトを用いて似たような実装ができます。そして、演算子という一部の機能ではオーバーロードといった関数の再定義が行えます。ここではイメージがしやすいように、より一般的なオーバーロードという言葉を用いて説明していきます。

オーバーロードにおいて大事なところは、はじめのstdでインポートした{Add、Div …}です。これらはある自分で実装した構造体に演算子を実装したいときに用いるものです。コードを見ると、これはAddトレイトというstdに含まれる提供されたトレイトを、自作したRcVariableに実装しています。また、type Output、といったあまり見慣れないコードもありますが、これはオーバーロードするトレイトの独自の機能です。ではここでAddの部分のコードがどのような振る舞いをするのか見てみましょう。

はじめに型を指定します。足し算といった2項演算子の場合、z = x・yという形になりますが、ここでself は左の x 、rhs は右の y 、Output は z の型を示します。これらは全て RcVariable なので、コードのように指定します。self については RcVariable にオーバーロードするので型を明示する必要はありません。そのあとは左と右のRcVariableを add関数 に渡すという処理です。これは、+ 演算子を用いると、add関数が自動で呼び出されるということです。 これらのオーバーロードを用いると、このようにコードをシンプルに書けるようになります。

fn main() {
    let a = RcVariable::new(1.0);
    let b = RcVariable::new(2.0);

    // 今まで足し算
    let c = add(a.clone(),b.clone());

    // オーバーロードを用いた足し算
    let c = a.clone() + b.clone();
}

二つのコードを比べてみてください。コードの量自体はあまり変わっていませんが、明らかに計算の構造が理解しやすくなりました。これはまだただの足し算なのであまり良さを感じませんが、今後複雑な関数を実装していくことで、このオーバーロードの恩恵を実感することになります。 この要領で他の演算子もオーバーロードしていきましょう。

neg は一変数関数なので、注意してください。また、オーバーロードに対応しているのはadd,mul,sub,div,negのみなので、pow はできません。よくわからない場合はstucrsのプログラムのlibをご覧ください。

TODO: RcVariableの追加実装の説明追加予定

RcVariableの追加実装

ここではよりRcVariableを便利に、そして可読性が増すように、RcVariableに独自のメソッドを追加していきます。

データからRcVariableを直接生成する関数

RcVariable を生成する際、今まではRcVariableのメソッドである new(data) という関数をもちいていました。実はこれをもっと簡単に、シンプルに生成する方法があります。前に用いた オーバーロード を用いるのです。※ここで言うオーバーロードも厳密にはオーバーロードではありません。

具体的にはf32型に自身のデータを保持するRcVariableを生成する関数をf32にオーバーロードするのです。試しにコードを作成しましょう。

pub trait F32ToRcVariable {
    fn rv(&self) -> RcVariable;
}

impl F32ToRcVariable for f32 {
    fn rv(&self) -> RcVariable {
        RcVariable::new(self.clone()) // <- 中の処理はRcVariableのnew()と同じ
    }
}

はじめにトレイトで関数を定義します。rv() という関数名です。続いてこの関数の処理を書き、f32に実装します。関数の処理はただのRcVariableの new() と同じです。これにより、f32型独自のメソッドを実装できました。ではこの関数を用いてRcVariableを生成してみましょう。

fn main() {
    let a = RcVariable::new(1.0); // <- 今までの生成方法
    let b = 2.0.rv();             // <- 新たな生成方法
}

この二つの生成を比べると、変数となる a,b がどのような値を持つRcVariableか一目で判断がつきます。今後このメソッドを多用しますので、忘れずに実装しておきましょう。

行列計算への前準備

いままで私たちは変数として主に「 スカラ 」を扱ってきました。しかし機械学習ではベクトルや行列などの「 テンソル 」を扱ってきます。本章ではテンソルを使うための前準備としてフレームワークを行列へと拡張する方法について説明していきます。 今回は多次元配列を扱うライブラリとしてndarrayというライブラリを使用します。 これはpythonでいうところの Numpy に相当します。使い方に関しては上のリンクを参照してください。

f32からArrayDへの変更

いままではf32型を用いてStuCrsフレームワークを構築していました。しかしこれから私たちは「テンソル」を扱いたいです。そこで「ndarray」ライブラリのArrayD型を使用しようと思います。

ArrayD型 とは実行中に次元数が決定する「テンソル」を扱う 動的な型です。 なぜArray1型やArray2型などの型を使わずArrayD型を使用するかというと私たちのこれから複雑な計算をさせる関数を実装するときに、ある関数では入力データが４次元であったり、また他の関数では入力データが２次元であったりするためです。これらのようにコンパイル時に次元数がわからず、様々な次元数を用いる関数を共通して管理したいためArrayD型を用います。 では実際にf32型をArrayD型に変更してみましょう。

Variableの変更

はじめに構造体の変更です。dataフィールド の型を f32 から ArrayD<f32> に変更します。

use ndarray::{array,ArrayD};
struct Variable {
    pub data: ArrayD<f32>,     // <- f32から変更
    grad: Option<RcVariable>,
    creator: Option<Rc<RefCell<dyn Function>>>,
    pub name: Option<String>,
    pub generation: i32,
    pub id: usize,
}

今回用いるArrayDの要素の型はf32に指定しました。ジェネリック型を用いれば、f64など、他の型にも対応できるようになりますが、コードが少し複雑になるため、f32に限定します。

impl Variable {
    pub fn new_rc(data: ArrayD<f32>) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Variable {
            data: data,
            grad: None,
            creator: None,
            name: None,
            generation: 0,
            id: id_generator(),
        }))
    }

    pub fn set_creator(&mut self, func: Rc<RefCell<dyn Function>>) {
        self.creator = Some(Rc::clone(&func));
        self.generation = func.borrow().get_generation() + 1;
    }

    fn backward(&self, double_grad: bool) {
        let mut funcs: Vec<Rc<RefCell<dyn Function>>> =
            vec![Rc::clone(self.creator.as_ref().unwrap())];

        let mut seen_set = HashSet::new();

        fn add_func(
            funcs_list: &mut Vec<Rc<RefCell<dyn Function>>>,
            seen_set: &mut HashSet<usize>,
            f: Rc<RefCell<dyn Function>>,
        ) {
            if seen_set.insert(f.borrow().get_id()) {
                funcs_list.push(Rc::clone(&f));
                funcs_list.sort_by(|a, b| {
                    a.borrow()
                        .get_generation()
                        .cmp(&b.borrow().get_generation())
                });
            }
        }
        //let first_grad = ArrayD::<f32>::ones(self.data.shape()).rv();

        //&selfで最初の変数はborrowされるので場合分け
        let mut last_variable = true;

        while let Some(f_rc) = funcs.pop() {

            let f_borrowed = f_rc.borrow();
            if double_grad == true {
                set_grad_true();
            } else {
                set_grad_false();
            }

            let xs = f_borrowed.get_inputs();

            let y = f_borrowed.get_output();

            let y_grad: RcVariable;

            if last_variable {
                y_grad = ArrayD::<f32>::ones(self.data.shape()).rv(); // <- 変更点

                last_variable = false;
            } else {
                //関数の出力は一つだけなので、[1]は必要なし
                y_grad = y.0.borrow().grad.as_ref().unwrap().clone();
            }

            let xs_grad = f_borrowed.backward(&y_grad);

            for (x, x_grad) in zip(xs, &xs_grad) {
                // gradをすでに保持しているなら、元のgradに新たなgradを足す。
                // gradをまだ持っていないならそれを持たせる。
                if let Some(current_grad_data) = x.grad() {
                    x.0.borrow_mut().grad = Some(current_grad_data + x_grad.clone());
                } else {
                    x.0.borrow_mut().grad = Some(x_grad.clone());
                }

                // creatorがあるならその関数をfuncsに追加
                if let Some(func_creator) = &x.0.borrow().creator {
                    add_func(&mut funcs, &mut seen_set, func_creator.clone());
                }
            }
        }
        
    }
}

続いてVariableに実装したメソッドの変更点です。変更する点は二つです。一つ目は new_rc() の引数です。 data の型をf32からArrayDにしてください。二つ目は y_grad です。矢印で示したところをArrayDにして対応させてください。

y_grad を変更する前に、下の説明の rv() の変更を先に行ってください。

f32からArrayDへの変更であまり変更するところがない理由は、データの値をVariableで扱っているため、Variableのデータさえ変更してしまえばあとは比較的同じように動作してくれるからです。 ArrayD型にも演算子などといった様々なメソッドが定義されているため、f32と同じように同じふるまいをしてくれるというわけです。

ArrayDからRcVariableを生成する関数

前回、f32からRcVariableを生成する関数、rv() を実装しましたが、これも ArrayD から RcVariable を生成する関数も実装しましょう。

//array型からRcVariable型を生成
pub trait ArrayDToRcVariable {
    fn rv(&self) -> RcVariable;
}
//arrayは任意の次元に対応
impl<D: Dimension> ArrayDToRcVariable for ArrayBase<OwnedRepr<f32>, D> {
    fn rv(&self) -> RcVariable {
        RcVariable::new(self.view().into_dyn())
    }
}

pub trait F32ToRcVariable {
    fn rv(&self) -> RcVariable;
}

//f32からarray型に変換し、rv()でRcVariableを生成
impl F32ToRcVariable for f32 {
    fn rv(&self) -> RcVariable {
        let array = array![*self as f32]; // <- f32型をArray型に変換
        array.rv()                        // <- Array型を通してRcVariableを生成
    }
}

ArrayD型にもf32型と同様に、 rv() メソッドを定義し、オーバーロードします。この際、ArrayDはf32とは異なりプリミティブ型ではないので、ただ単に渡せばよいものではありません。 view() や into_dyn() に関しては ndarray のドキュメントを参照してください。Array型からRcVariableへの変換を基準にして、f32からの変換の場合は、f32型をArray型に変換してからRcVariableを生成するという流れにします。

fn main() {
    let a = array![2.0,3.0].rv(); // <- 新たな生成方法
}

するとこのようにArray型の行列からRcVariableを生成できるようになりました。

各関数の行列への対応

私たちはこれまで様々な関数をStuCrsに実装してきました。例えばadd,mul,sin関数などです。しかしながら私たちはそれらの関数を実装するにあたって入力と出力がすべて「スカラ」であることを想定してきました。なので「テンソル」として各関数を行列へと拡張していきましょう。今回はsin関数を例にとって説明したいと思います。

まずは理論的なところから説明します。今まではxという「スカラ」にsin関数を適用する場合sin(x)とすればよかったわけです。もしｘが「テンソル」の場合、たとえば行列の場合はどうなるでしょうか。その場合は要素ごとにsin関数が適用されます。

例えばｘ＝[1,2,3,4,5,6]という2次元行列にsin関数を適用した場合、下に表したように要素ごとにsin関数が適用されるわけです。このことを念頭において実装してみましょう。

$$ Sin(X): \begin{bmatrix} 1 & 2 & 3\\ 4 & 5 & 6 \end{bmatrix} \mapsto \begin{bmatrix} sin(1) & sin(2) & sin(3)\\ sin(4) & sin(5) & sin(6) \end{bmatrix} $$

※変更する関数のみ表示します。

impl Function for Sin {
    fn call(&mut self) -> RcVariable {
        let inputs = &self.inputs;
        if inputs.len() != 1 {
            panic!("Sinは一変数関数です。inputsの個数が一つではありません。")
        }

        let output = self.forward(inputs);

        if get_grad_status() == true {
            //inputのgenerationで一番大きい値をFuncitonのgenerationとする
            self.generation = inputs.iter().map(|input| input.generation()).max().unwrap();

            //  outputを弱参照(downgrade)で覚える
            self.output = Some(output.downgrade());

            let self_f: Rc<RefCell<dyn Function>> = Rc::new(RefCell::new(self.clone()));

            //outputsに自分をcreatorとして覚えさせる
            output.0.borrow_mut().set_creator(self_f.clone());
        }

        output
    }

    fn forward(&self, xs: &[RcVariable]) -> RcVariable {
        let x = &xs[0];
        let y_data = x.data().mapv(|x| x.sin()); // <- mapv()を用いて要素ごとに計算

        y_data.rv()
    }

    fn backward(&self, gy: &RcVariable) -> Vec<RcVariable> {
        let x = &self.inputs[0];
        let gx = cos(x) * gy.clone(); // <- Cos構造体を呼び出してRcVariableとして処理
        let gxs = vec![gx];
        gxs
    }
}

Function構造体でのArrayD型への変更は実はほとんどありません。なぜなら、Function構造体はデータをRcVariableですべて管理しているからです。前回Variable、RcVariableへの対応は済みましたので、特に変更する点はありません。唯一変更する点は forward の計算処理の変更です。forwardの中はArray型で計算処理されるので、変更しなければならないところが生まれます。Array型のメソッドの mapv() を用いて計算します。backwardの処理は全て、RcVariableで処理されるので、変更点はありません。

一方で四則演算の関数はどうでしょうか。実はforwardの処理のところでさえ変更する必要はありません。なぜなら、Array型は演算子に対応しているからです。 つまり、Array型でそのまま+や-が使えるということです。これは私たちが先ほど RcVariable に実装した演算子のオーバーロードと同じです。あの演算子のメソッドが自動で呼び出されるので、f32の記法のままで同じ処理をしてくれます。

ここで重要なのは微分する際のbackwardのときのcosの計算も要素ごとに行わなければならないということです。そこで、cosも同じように構造体を変更して用いれば、正しく微分を行えます。図で表すとこのようになります。

$$ Forward: \begin{bmatrix} x_0 & x_1 & x_2\\ x_3 & x_4 & x_5 \end{bmatrix} \xrightarrow{Sin} \begin{bmatrix} sin(x_0) & sin(x_1) & sin(x_2)\\ sin(x_3) & sin(x_4) & sin(x_5) \end{bmatrix} $$

$$ Backward: \begin{bmatrix} cos(x_0) \cdot gy_0 & cos(x_1) \cdot gy_1 & cos(x_2) \cdot gy_2\\ cos(x_3) \cdot gy_3 & cos(x_4) \cdot gy_4 & cos(x_5) \cdot gy_5 \end{bmatrix} \xleftarrow{Sin’} \begin{bmatrix} gy_0 & gy_1 & gy_2\\ gy_3 & gy_4 & gy_5 \end{bmatrix} $$

今までスカラーでバックプロパゲーションを行っていたのが、行列を導入することで、要素ごとのバックプロパゲーションを実現することができました。

ではSin関数と同じように他の関数も同様に変更してみましょう。四則演算の関数は変更するところはありません。

行列特有の関数の実装

Reshape関数の実装

graph LR
 arrx["$$x:\begin{pmatrix}x_0 & x_1 & x_2 \\\ x_3 & x_4 & x_5\end{pmatrix}$$"] --> reshape[Reshape]
 reshape --> arry["$$y:\begin{pmatrix}x_0 & x_1 & x_2 & x_3 & x_4 & x_5\end{pmatrix}$$"]
 
 style arrx stroke-width:0px
 style arry stroke-width:0px

graph RL
 arrgy["$$gy:\begin{pmatrix}gy_0 & gy_1 & gy_2 & gy_3 & gy_4 & gy_5\end{pmatrix}$$"] --> reshape["Reshape'"]
 reshape --> arrgx["$$gx:\begin{pmatrix}gy_0 & gy_1 & gy_2 \\\ gy_3 & gy_4 & gy_5\end{pmatrix}$$"] 

 style arrgx stroke-width:0px
 style arrgy stroke-width:0px

reshape関数は行列の形状だけを変更する関数です。つまり数値計算を行わず、値をそのまま流すだけなので、微分の値は1であり、上流からきた微分の値をそのまま流します。このような微分が1の関数は今後も行列計算では出てきます。ここで重要なことは、返す微分の行列はinputの行列と同じ形状でなければならないことです。 つまり、xとgx、yとgyの形状はそれぞれ一致していなければなりません。 なので、backwardの場合、形状を戻して渡す必要があるのです。

struct Reshape {
    inputs: Vec<RcVariable>,
    output: Option<Weak<RefCell<Variable>>>,
    shape: IxDyn,
    generation: i32,
    id: usize,
}

impl Function for Reshape {
    fn call(&mut self) -> RcVariable {
        let inputs = &self.inputs;
        if inputs.len() != 1 {
            panic!("Reshapeは一変数関数です。inputsの個数が一つではありません。")
        }

        let output = self.forward(inputs);

        if get_grad_status() == true {
            //inputのgenerationで一番大きい値をFuncitonのgenerationとする
            self.generation = inputs.iter().map(|input| input.generation()).max().unwrap();

            //  outputを弱参照(downgrade)で覚える
            self.output = Some(output.downgrade());

            let self_f: Rc<RefCell<dyn Function>> = Rc::new(RefCell::new(self.clone()));

            //outputsに自分をcreatorとして覚えさせる
            output.0.borrow_mut().set_creator(self_f.clone());
        }

        output
    }

    fn forward(&self, xs: &[RcVariable]) -> RcVariable {
        let x = &xs[0];
        let y_shape = self.shape.clone();
        let y_data = x.data().to_shape(y_shape).unwrap().to_owned();

        y_data.rv()
    }

    fn backward(&self, gy: &RcVariable) -> Vec<RcVariable> {
        let x = &self.inputs[0];
        let x_shape = IxDyn(x.data().shape());
        let gx = reshape(gy, x_shape);
        let gxs = vec![gx];

        gxs
    }

    fn get_inputs(&self) -> &[RcVariable] {
        &self.inputs
    }

    fn get_output(&self) -> RcVariable {
        let output;
        output = self
            .output
            .as_ref()
            .unwrap()
            .upgrade()
            .as_ref()
            .unwrap()
            .clone();

        RcVariable(output)
    }

    fn get_generation(&self) -> i32 {
        self.generation
    }
    fn get_id(&self) -> usize {
        self.id
    }
}
impl Reshape {
    fn new(inputs: &[RcVariable], shape: IxDyn) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Self {
            inputs: inputs.to_vec(),
            output: None,
            shape: shape,
            generation: 0,
            id: id_generator(),
        }))
    }
}

fn reshape_f(xs: &[RcVariable], shape: IxDyn) -> RcVariable {
    Reshape::new(xs, shape).borrow_mut().call()
}

pub fn reshape(x: &RcVariable, shape: IxDyn) -> RcVariable {
    let y = reshape_f(&[x.clone()], shape);
    y
}

backward で形状を戻すために、inputの形状を覚えておく必要がありますが、inputのデータは保持しており、そこから呼び出せばよいので特に考えません。

変形させたい形状を引数として受け取るために、呼び出す関数の引数に shape を設定します。このshapeの型は IxDyn で動的な形状を指します。また、それを保持するために、shapeフィールド を持ちます。フィールドとして保持し、forwardでinputのArrayのデータの形状を変形します。backwardでは上流からきた微分の値である gy をinputの形状に変形させるために、xを呼び出して形状を求め、それによって変形させます。backward内の reshape 関数は今まさに実装している Reshape構造体 です。

では実装した Reshape 関数をテストしてみましょう。微分の値や形状などに着目してください。

#[test]
    fn reshape_test() {
        use crate::core_new::ArrayDToRcVariable;

        let a = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]].rv();
        let b = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]].rv();

        let mut y0 = reshape(&a, Dim(IxDyn(&[1, 6])));

        let mut y1 = b.reshape(Dim(IxDyn(&[1, 6])));

        println!("y0 = {}", y0.data()); //[[1,2,3,4,5,6]] shape(1,6)
        println!("y1 = {}", y1.data()); //[[1,2,3,4,5,6]] shape(1,6)

        y0.backward(false);
        y1.backward(false);
        println!("a_grad = {:?}", a.grad().unwrap().data()); // [[1.0,1.0,1.0],[1.0,1.0,1.0]] shape(2,3)
        println!("b_grad = {:?}", b.grad().unwrap().data()); // [[1.0,1.0,1.0],[1.0,1.0,1.0]] shape(2,3)
    }

すると、上の図のようにaとa_gradの形状が一致しているのがわかります。

Transpose関数の実装

graph LR
 arrx["$$x:\begin{pmatrix}x_0 & x_1 & x_2 \\\ x_3 & x_4 & x_5\end{pmatrix}$$"] --> f[Transpose]
 f --> arry["$$y:\begin{pmatrix}x_0 & x_3 \\\  x_1 & x_4\\\ x_2 & x_5 \end{pmatrix}$$"]
 
 style arrx stroke-width:0px
 style arry stroke-width:0px

graph RL
 arrgy["$$gy:\begin{pmatrix}gy_0 & gy_3 \\\  gy_1 & gy_4\\\ gy_2 & gy_5\end{pmatrix}$$"] --> reshape["Transpose'"]
 reshape --> arrgx["$$gx:\begin{pmatrix}gy_0 & gy_1 & gy_2 \\\ gy_3 & gy_4 & gy_5\end{pmatrix}$$"] 
 

 style arrgx stroke-width:0px
 style arrgy stroke-width:0px

Transpose関数はinputの行列の転置行列を返す関数です。この関数も形状を変更するだけの関数であり、微分といった原理はReshapeと同じです。繰り返し話しますが、行列計算におけるバックプロパゲーションの重要なところは、形状が一致するように戻すことです。

struct Transpose {
    inputs: Vec<RcVariable>,
    output: Option<Weak<RefCell<Variable>>>,
    generation: i32,
    id: usize,
}

impl Function for Transpose {
    fn call(&mut self) -> RcVariable {
        let inputs = &self.inputs;
        if inputs.len() != 1 {
            panic!("Transposeは一変数関数です。inputsの個数が一つではありません。")
        }

        let output = self.forward(inputs);

        if get_grad_status() == true {
            //inputのgenerationで一番大きい値をFuncitonのgenerationとする
            self.generation = inputs.iter().map(|input| input.generation()).max().unwrap();

            //  outputを弱参照(downgrade)で覚える
            self.output = Some(output.downgrade());

            let self_f: Rc<RefCell<dyn Function>> = Rc::new(RefCell::new(self.clone()));

            //outputsに自分をcreatorとして覚えさせる
            output.0.borrow_mut().set_creator(self_f.clone());
        }

        output
    }

    fn forward(&self, xs: &[RcVariable]) -> RcVariable {
        let x = &xs[0];
        let y_data = x.data().t().to_owned();

        y_data.rv()
    }

    fn backward(&self, gy: &RcVariable) -> Vec<RcVariable> {
        let gxs = vec![gy.t().to_owned()];

        gxs
    }

    fn get_inputs(&self) -> &[RcVariable] {
        &self.inputs
    }

    fn get_output(&self) -> RcVariable {
        let output;
        output = self
            .output
            .as_ref()
            .unwrap()
            .upgrade()
            .as_ref()
            .unwrap()
            .clone();

        RcVariable(output)
    }

    fn get_generation(&self) -> i32 {
        self.generation
    }
    fn get_id(&self) -> usize {
        self.id
    }
}
impl Transpose {
    fn new(inputs: &[RcVariable]) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Self {
            inputs: inputs.to_vec(),
            output: None,
            generation: 0,
            id: id_generator(),
        }))
    }
}

fn transpose_f(xs: &[RcVariable]) -> RcVariable {
    Transpose::new(xs).borrow_mut().call()
}

pub fn transpose(x: &RcVariable) -> RcVariable {
    let y = transpose_f(&[x.clone()]);
    y
}

転置行列の形状変更の手段は一つに決まっており、backward でも上流からきた微分の値を転置させればinputと同じ形状に戻るので特に形状を覚えるといった操作は必要ありません。

forwardでinputのArrayのデータの形状を転置行列として変形します。backwardでは上流からきた微分の値である gy を再び転置します。backward内の gy.t() は Transpose構造体 をRcVariableのメソッドとして呼び出しています。

では実装した Transpose 関数をテストしてみましょう。微分の値や形状などに着目してください。

#[test]
    fn transpose_test() {
        use crate::core_new::ArrayDToRcVariable;

        let a = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]].rv();
        let b = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]].rv();

        let mut y0 = transpose(&a);

        let mut y1 = b.t();

        println!("y0 = {}", y0.data()); //[[1,4],[2,5],[3,6]] shape(3,2)
        println!("y1 = {}", y1.data()); //[[1,4],[2,5],[3,6]] shape(3,2)

        y0.backward(false);
        y1.backward(false);
        println!("a_grad = {:?}", a.grad().unwrap().data()); // [[1.0,1.0,1.0],[1.0,1.0,1.0]] shape(2,3)
        println!("b_grad = {:?}", b.grad().unwrap().data()); // [[1.0,1.0,1.0],[1.0,1.0,1.0]] shape(2,3)
    }

すると、上の図のようにaとa_gradの形状が一致しているのがわかります。

Transpose関数 として軸を指定し、軸を入れ替える関数として定義されているものもあります。今回のこの転置という処理は軸0と軸1を入れ替えたものと捉えることもできます。これが3次元の行列などになると、3つの軸を入れ替えるという作業も必要となります。具体的には CNN の関数で使われます。そのような処理も正しくバックプロパゲーションできる関数を今後 PermuteAxes として別の関数で定義します。本来ならば、軸の入れ替えを統合的に行える処理をこのTranspose 関数で実装すべきでしたが、このような処理は複雑なため、分離させた方が良いと考えました。

Sum関数の実装

graph LR
 arrx["$$x:\begin{pmatrix}0 & 1 & 2 \\\ 3 & 4 & 5\end{pmatrix}$$"] --> f[Sum]
 f --> arry["$$y:\begin{pmatrix}21 \end{pmatrix}$$"]
 
 style arrx stroke-width:0px
 style arry stroke-width:0px

graph RL
 arrgy["$$gy:\begin{pmatrix}gy\end{pmatrix}$$"] --> reshape["Sum'"]
 reshape --> arrgx["$$gx:\begin{pmatrix}gy & gy & gy \\\ gy & gy & gy\end{pmatrix}$$"] 
 

 style arrgx stroke-width:0px
 style arrgy stroke-width:0px

Sum関数はinputの行列のを合計値を計算し、スカラーとして出力する関数です。この関数は足し算として数値計算を行う関数ですが、\(y=x_0+x_1+x_2+x_3+x_4+x_5\)なので、\(\frac{\partial y}{\partial x_k}=1\)となり、gy はすべて1となります。微分の原理はReshapeと同じです。繰り返し話しますが、行列計算におけるバックプロパゲーションの重要なところは、形状が一致するように戻すことです。

Sum関数 は行列の要素の和を足し合わせてそれを返す関数です。このとき要素数が減るので、自然と形状は変化しますが、ここで重要なのがどの要素をどのように足すかということです。具体的には軸を指定することでいくつかの足し方を表すことができます。
ここで軸を指定したときのいくつかの例を見てみましょう。

① 軸指定なし

graph LR
 arrx["$$x:\begin{pmatrix}1 & 2 & 3 \\\ 4 & 5 & 6\end{pmatrix}$$"] --> f[Sum]
 f --> arry["$$y:\begin{pmatrix}21 \end{pmatrix}$$"]
 
 style arrx stroke-width:0px
 style arry stroke-width:0px

graph RL
 arrgy["$$gy:\begin{pmatrix}1\end{pmatrix}$$"] --> reshape["Sum'"]
 reshape --> arrgx["$$gx:\begin{pmatrix}1 & 1 & 1 \\\ 1 & 1 & 1\end{pmatrix}$$"] 
 

 style arrgx stroke-width:0px
 style arrgy stroke-width:0px

② 軸(axis) = 0

graph LR
 arrx["$$x:\\downarrow\\begin{pmatrix}1 & 2 & 3 \\\ 4 & 5 & 6\end{pmatrix}$$"] --> f[Sum]
 f --> arry["$$y:\begin{pmatrix}5 & 7 & 9 \end{pmatrix}$$"]
 
 style arrx stroke-width:0px
 style arry stroke-width:0px

graph RL
 arrgy["$$gy:\begin{pmatrix}1 & 1& 1\end{pmatrix}$$"] --> reshape["Sum'"]
 reshape --> arrgx["$$gx:\begin{pmatrix}1 & 1 & 1 \\\ 1 & 1 & 1\end{pmatrix}$$"] 
 

 style arrgx stroke-width:0px
 style arrgy stroke-width:0px

③ 軸(axis) = 1

graph LR
 arrx["$$x:\underrightarrow{\begin{pmatrix}1 & 2 & 3 \\\ 4 & 5 & 6\end{pmatrix}}$$"] --> f[Sum]
 f --> arry["$$y:\begin{pmatrix}6 \\\15 \end{pmatrix}$$"]
 
 style arrx stroke-width:0px
 style arry stroke-width:0px

graph RL
 arrgy["$$gy:\begin{pmatrix}1 \\\1 \end{pmatrix}$$"] --> reshape["Sum'"]
 reshape --> arrgx["$$gx:\begin{pmatrix}1 & 1 & 1 \\\ 1 & 1 & 1\end{pmatrix}$$"] 
 

 style arrgx stroke-width:0px
 style arrgy stroke-width:0px

ここでは2次元の行列におけるSum関数を例に挙げました。軸の指定とは、その軸に沿って要素を足すことです。上の図の x の行列には矢印が書かれていますが、その矢印の向きに足すということです。例えば、軸=0のとき、縦を表すのでその軸に沿って足すと形状は、(2,3) → (1,3)となります。軸1も同じ考え方です。

軸の指定の方法について、今後も多く扱うため、このドキュメントに関して統一しておきたいと思います。ある行列の形状が(m,n)の場合、この配列のインデックスで軸を決めます。軸＝0ならmを、1ならnを指します。このように定義すれば、これから3次元以上となったとしても一般化して規則的に扱うことができます。

struct Sum {
    inputs: Vec<RcVariable>,
    output: Option<Weak<RefCell<Variable>>>,
    axis: Option<u16>,
    generation: i32,
    id: usize,
}

impl Function for Sum {
    fn call(&mut self) -> RcVariable {
        let inputs = &self.inputs;
        if inputs.len() != 1 {
            panic!("Sumは一変数関数です。inputsの個数が一つではありません。")
        }

        let output = self.forward(inputs);

        if get_grad_status() == true {
            //inputのgenerationで一番大きい値をFuncitonのgenerationとする
            self.generation = inputs.iter().map(|input| input.generation()).max().unwrap();

            //  outputを弱参照(downgrade)で覚える
            self.output = Some(output.downgrade());

            let self_f: Rc<RefCell<dyn Function>> = Rc::new(RefCell::new(self.clone()));

            //outputsに自分をcreatorとして覚えさせる
            output.0.borrow_mut().set_creator(self_f.clone());
        }

        output
    }

    fn forward(&self, xs: &[RcVariable]) -> RcVariable {
        let x = &xs[0];
        let axis = self.axis;

        let y_data;

        if let Some(axis_data) = axis {
            if axis_data != 0 && axis_data != 1 {
                panic!("axisは0か1の値のみ指定できます")
            }

            y_data = x
                .data()
                .sum_axis(Axis(axis_data as usize))
                .insert_axis(Axis(axis_data as usize));
        } else {
            let scalar_y = x.data().sum();
            y_data = array![scalar_y].into_dyn();
        }

        y_data.rv()
    }

    fn backward(&self, gy: &RcVariable) -> Vec<RcVariable> {
        let x = &self.inputs[0];
        let x_shape = IxDyn(x.data().shape());
        let gx = broadcast_to(gy, x_shape);
        let gxs = vec![gx];

        gxs
    }

    fn get_inputs(&self) -> &[RcVariable] {
        &self.inputs
    }

    fn get_output(&self) -> RcVariable {
        let output;
        output = self
            .output
            .as_ref()
            .unwrap()
            .upgrade()
            .as_ref()
            .unwrap()
            .clone();

        RcVariable(output)
    }

    fn get_generation(&self) -> i32 {
        self.generation
    }
    fn get_id(&self) -> usize {
        self.id
    }
}
impl Sum {
    fn new(inputs: &[RcVariable], axis: Option<u16>) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Self {
            inputs: inputs.to_vec(),
            output: None,
            axis: axis,

            generation: 0,
            id: id_generator(),
        }))
    }
}

fn sum_f(xs: &[RcVariable], axis: Option<u16>) -> RcVariable {
    Sum::new(xs, axis).borrow_mut().call()
}

pub fn sum(x: &RcVariable, axis: Option<u16>) -> RcVariable {
    let y = sum_f(&[x.clone()], axis);
    y
}

このsum関数の軸指定は Array型 にsumメソッドとして標準装備されているので、そのまま使用します。軸指定の値は Option<u16> としてフィールドで保持します。 None が軸指定なし、0,1という整数の数字が軸の値を指します。ここで軸の数値について、軸の個数、すなわち次元数以上の軸は指定できないので、それ以上の値は渡せないように panicやResult型 で設定します。Sum関数は2次元の行列を対象としています。今回の場合、2次元なので0から数えて1までの二つまで指定できます。よって2以上の値は指定できません。

では実装した Sum 関数を上の3つの場合に分けてテストしてみましょう。微分の値や形状などに着目してください。

#[test]
    fn sum_test() {
        use crate::core_new::ArrayDToRcVariable;

        let a = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]].rv();
        let b = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]].rv();
        let c = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]].rv();

        let mut y0 = sum(&a, None);
        let mut y1 = sum(&b, Some(0));
        let mut y2 = sum(&c, Some(1));

        println!("y0 = {}", y0.data()); // 21.0
        println!("y1 = {}", y1.data()); //
        println!("y2 = {}", y2.data()); //

        y0.backward(false);
        y1.backward(false);
        y2.backward(false);

        println!("a_grad = {:?}", a.grad().unwrap().data()); //
        println!("b_grad = {:?}", b.grad().unwrap().data()); //
        println!("c_grad = {:?}", c.grad().unwrap().data()); //
    }

すると、上の図のように求めるgxと一致しているのがわかります。

MatMul関数の実装

行列計算において非常に重要な計算として行列積があります。行列積はベクトルによるベクトル積の拡張であり、深層学習の線形変換、ひいては全結合型のニューラルネットワークを構築するという重要な役割を果たしています。 このドキュメントを読んでいる方はすでに行列の扱い方に詳しいと思いますが、ここで改めて行列積について考え、さらに行列積のバックプロパゲーションの原理を導き出しましょう。

はじめに、行列積を深層学習でなぜ用いるのか考えてみましょう。

graph LR
 x0(("$$x_0$$")) -->|"$$w_0$$"| y0(("$$y_0$$"))
 x0(("$$x_0$$")) -->|"$$w_1$$"| y1(("$$y_1$$"))
 x1(("$$x_1$$")) -->|"$$w_2$$"| y0
 x1(("$$x_1$$")) -->|"$$w_3$$"| y1

これはニューラルネットワークの基本となるニューロンの図式です。これらの関係は、 入力の\(X = (x_0,x_1)\)、重みの\(W\)、出力の\(Y = (y_0,y_1)\)が存在しますが、出力と入力の関係は、\(y_0 = x_0\cdot w_0 + x_1\cdot w_2\)、\(y_1 = x_0\cdot w_1 + x_1\cdot w_3\) というような数式で表せます。これはまさに、

$$(y_0,y_1) = (x_0,x_1)* \begin{pmatrix} w_0 & w_1 \\ w_2 & w_3 \end{pmatrix} $$
という行列積の計算です。このとき$X$はベクトルですが、入力の値が複数あった場合、$X$を行列として拡張すれば、行列積という一つの関係式で表せることができます。これらのニューロンは一つでは表現できる式は浅いものとなってしまいます。このニューロンがたくさんつながりあい、深い層となって深い表現が可能となり、精度向上につながるのです。そして、たくさんあるニューロンの計算を効率よく行う手段として行列積を用いるということです。

ではこの行列積を実現するFunction構造体を実装していきます。はじめに行列積を扱う関数MatMulの式の関係はこのようになります。

graph RL
 L(("$$L$$")) --> ...
 ... --> Y
 Y --> ...
 ... --> L
 Y(("$$Y$$")) -->|"$$\frac{\partial L}{\partial Y}$$"| matmul[Matmul]
 matmul --> Y
 X(("$$X$$")) --> matmul
 matmul -->|"$$\frac{\partial L}{\partial X}$$"| X
 matmul -->|"$$\frac{\partial L}{\partial W}$$"| W(("$$W$$")) 
 W --> matmul

まず行列積\(Y = X\cdot W\)という関数を考えます。行列積という演算は非可換(\(a\cdot b \ne b\cdot a\))なので、かける順番を\(X,W\)という順番で正確に定義しておきます。なので引数として行列を渡す際、順番を考えて渡さなくてはなりません。そうでなければ、順伝播、逆伝播においても計算結果が異なってしまいます。

この行列積を行う計算処理はArray型で装備されているので、試しに使用してみましょう。

TODO: Arrayで行列積を行うコード追加

Array型で提供されるdotメソッドは2次元の行列までしか対応していません。よって、Matmul関数も2次元の行列のみに対応します。3次元の行列積、いわゆるテンソル積 はCNNといった機能で必要となるため、今後TensorDot として別のFunction構造体を定義します。

では次に行列積のバックプロパゲーションを考えてみます。はじめに各変数の偏微分の式はこのようになります。

$$ \begin{align} \frac{\partial L}{\partial X} = \frac{\partial L}{\partial Y} \cdot W^\mathsf{T}\\ \frac{\partial L}{\partial W} = X^\mathsf{T} \cdot \frac{\partial L}{\partial Y} \end{align} $$

転置行列などが用いられるなど、直感に反する面があります。この式に登場する\(\frac{\partial L}{\partial Y}\)は上流からきた微分の値を指します。

本来ならば、この偏微分の式を一からこの場で解説したいところですが、筆者である私自身、LaTeXの記法に慣れていないので、この説明に関してはより詳しく、そしてわかりやすく解説しているサイトにゆだねたいと思います。これに関しては、一度自分の手で実際に数式などを考え、偏微分で転置行列が出てくるところなどを試していただきたいです。

偏微分の数学的な説明を省略してしまいましたが、ここで確認してほしい重要な点は、形状の確認です。行列積はあらゆる行列の組み合わせにおいて行える計算処理ではありません。二つの行列には、形状においてある条件を満たしていなければいけません。

$$(y_0,y_1) = (x_0,x_1)* \begin{pmatrix} w_1 & w_2 \\ w_3 & w_4 \end{pmatrix} $$ この時、\(y_0\)はいわばベクトル\((x_0,x_1)\)とベクトル\((w_1,w_3)\)のベクトル積です。ベクトル積で求めるということではこの二つのベクトルの長さは等しくなければなりません。つまり前の行列の列数と後の行列の行数が一致していなければならないという条件が必要です。これがまさに、行列積を行える形状の条件です。これを一般化してみます。

形状が\((k,l)\)の行列\(A\)と、\((m,n)\)の行列\(B\)における行列積では、\(A\)の列数と\(B\)の行数、すなわ\(l\)と\(m\)が一致していなければなりません。またこの条件を満たしたうえで、行列積を行った場合の出力の行列\(C\)の形状は\((k,n)\)となります。出力の形状がこのようになるのは明らかなので証明は省略します。

この考えを用いれば、行列のバックプロパゲーションで重要な変数のdataとgradの形状の一致も先ほどの偏微分の式から分かります。つまり、\(\frac{\partial L}{\partial X}\)は\(X\)と、\(\frac{\partial L}{\partial W}\)は\(W\)と形状が一致するということです。

前置きが長くなりましたが、それではMatMul構造体を実装してみます。

struct MatMul {
    inputs: Vec<RcVariable>,
    output: Option<Weak<RefCell<Variable>>>,
    generation: i32,
    id: usize,
}

impl Function for MatMul {
    fn call(&mut self) -> RcVariable {
        let inputs = &self.inputs;
        if inputs.len() != 2 {
            panic!("Matmulは二変数関数です。inputsの個数が二つではありません。")
        }

        let output = self.forward(inputs);

        if get_grad_status() == true {
            //inputのgenerationで一番大きい値をFuncitonのgenerationとする
            self.generation = inputs.iter().map(|input| input.generation()).max().unwrap();

            //  outputを弱参照(downgrade)で覚える
            self.output = Some(output.downgrade());

            let self_f: Rc<RefCell<dyn Function>> = Rc::new(RefCell::new(self.clone()));

            //outputsに自分をcreatorとして覚えさせる
            output.0.borrow_mut().set_creator(self_f.clone());
        }

        output
    }

    fn forward(&self, xs: &[RcVariable]) -> RcVariable {
        //xs[0]の方をX, xs[1]の方をWとする
        let x = &xs[0];
        let w = &xs[1];

        let x_data = x.data();
        let w_data = w.data();

        //match以降の場合分けを関数にしたい
        let y_data = array_matmul(&x_data.view(), &w_data.view());

        y_data.rv()
    }

    fn backward(&self, gy: &RcVariable) -> Vec<RcVariable> {
        let x = &self.inputs[0];
        let w = &self.inputs[1];

        let gx = matmul(gy, &w.t());
        let gw = matmul(&x.t(), gy);
        let gxs = vec![gx, gw];

        gxs
    }

    fn get_inputs(&self) -> &[RcVariable] {
        &self.inputs
    }

    fn get_output(&self) -> RcVariable {
        let output;
        output = self
            .output
            .as_ref()
            .unwrap()
            .upgrade()
            .as_ref()
            .unwrap()
            .clone();

        RcVariable(output)
    }

    fn get_generation(&self) -> i32 {
        self.generation
    }
    fn get_id(&self) -> usize {
        self.id
    }
}
impl MatMul {
    fn new(inputs: &[RcVariable]) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Self {
            inputs: inputs.to_vec(),
            output: None,
            generation: 0,
            id: id_generator(),
        }))
    }
}

pub fn array_matmul(x_array: &ArrayViewD<f32>, w_array: &ArrayViewD<f32>) -> ArrayD<f32> {
    let y = match (x_array.ndim(), w_array.ndim()) {
        // 1D × 1D → スカラー出力
        (1, 1) => {
            let x = x_array.clone().into_dimensionality::<Ix1>().unwrap();
            let w = w_array.clone().into_dimensionality::<Ix1>().unwrap();

            let y = x.dot(&w);
            ArrayD::from_elem(ndarray::IxDyn(&[]), y) // スカラーとして返す
        }

        // 2D × 1D
        (2, 1) => {
            let x = x_array.clone().into_dimensionality::<Ix2>().unwrap();
            let w = w_array.clone().into_dimensionality::<Ix1>().unwrap();
            let y = x.dot(&w);
            y.into_dyn()
        }

        // 1D × 2D
        (1, 2) => {
            let x = x_array.clone().into_dimensionality::<Ix1>().unwrap();
            let w = w_array.clone().into_dimensionality::<Ix2>().unwrap();
            let y = x.dot(&w);
            y.into_dyn()
        }

        // 2D × 2D
        (2, 2) => {
            let x = x_array.clone().into_dimensionality::<Ix2>().unwrap();
            let w = w_array.clone().into_dimensionality::<Ix2>().unwrap();
            let y = x.dot(&w);
            y.into_dyn()
        }

        _ => {
            panic!("3次元以上の行列積は未実装");
        }
    };

    y
}

fn matmul_f(xs: &[RcVariable]) -> RcVariable {
    MatMul::new(xs).borrow_mut().call()
}

pub fn matmul(x: &RcVariable, w: &RcVariable) -> RcVariable {
    let y = matmul_f(&[x.clone(), w.clone()]);
    y
}

基本的には Add や Mul といった二変数関数の変数を X,W に変更し、Array型の dotメソッド を用いて計算を行うというのが主な変更点です。注意したい点は backward と dot の処理方法です。

backwardの方は先ほどの偏微分の式に当てはめて計算します。\(\frac{\partial L}{\partial X} = \frac{\partial L}{\partial Y} \cdot W^\mathsf{T}\)の場合、上流からきた微分の値である\(\frac{\partial L}{\partial Y}\)と\({W}^\mathsf{T}\)の行列積を行います。\(\frac{\partial L}{\partial W}\)も同様です。

Array型のメソッドの dot は実は静的な形状のArray型にしか対応していません。私はこのニューラルネットワークにおいて行列は動的な形状である ArrayD として扱ってきたため、そのまま代入することはできません。よって fn array_matmul() 関数定義して用いることで動的な行列を静的に変更します。静的にするため、場合分けし、 into_dimensionality() で静的に変換します。私たちは動的な行列を基本としｔますので、dot計算したあとは into_dyn() で再び動的に戻します。これらの変換はパフォーマンスにはほとんど影響を及ぼしません。

では実装した Matmul 関数をテストしてみましょう。二変数ということに注目してください。

#[test]
    fn matmul_test() {
        use crate::core_new::ArrayDToRcVariable;

        let a = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]].rv();

        let b = array![[7.0, 8.0], [9.0, 10.0], [11.0, 12.0]].rv();

        let mut y = matmul(&a, &b);

        println!("y = {}", y.data()); // [[58,64],[139,154]]
        y.backward(false);

        println!("a_grad = {:?}", a.grad().unwrap().data());
        println!("b_grad = {:?}", b.grad().unwrap().data());
    }

ブロードキャスト対応

行列計算においてブロードキャストという機能があります。これは簡単に言えば、行列の演算が行われるよう、自動で計上を変形させる機能です。
ここでブロードキャストの例を見てみます。

①

graph LR
 arr_a("$$A:\begin{pmatrix}1 \end{pmatrix}$$")  --> 
 arr_ab["$$A':\begin{pmatrix}1 & 1 & 1 \\\ 1 & 1 & 1\end{pmatrix}$$"]
 arr_ab --> Add[Add]
 arr_b["$$B:\begin{pmatrix}1 & 2 & 3 \\\ 4 & 5 & 6\end{pmatrix}$$"] --> Add
 Add --> arr_c["$$C:\begin{pmatrix}2 & 3 & 4 \\\ 5 & 6 & 7\end{pmatrix}$$"]
 
 style arr_a stroke-width:0px
 style arr_ab stroke-width:0px
 style arr_b stroke-width:0px
 style arr_c stroke-width:0px

この行列計算を考えます。以前の説明により、要素ごとの計算は同じ形状である必要がありました。なので、この行列の四則演算も同じように、A と B は同じ形状でなければなりません。しかし、四則演算は多用するためユーザービリティの観点から 自動で形状を揃えることで、異なる形状の演算を可能にします。そして、この形状を揃える機能はArray型 に装備されており、暗黙的に変換されます。上の計算の場合、AをBに揃えるような形で自動でA’に拡張することで、足し算を実現します。

次の例も見てみましょう。

②

graph LR
 arr_a("$$A:\begin{pmatrix}1 &2 &3\end{pmatrix}$$")  --> 
 arr_ab["$$A'':\begin{pmatrix}1 & 2 & 3 \\\ 1 & 2 & 3\end{pmatrix}$$"]
 arr_ab --> Add[Add]
 arr_b["$$B:\begin{pmatrix}1 & 2 & 3 \\\ 4 & 5 & 6\end{pmatrix}$$"] --> Add
 Add --> arr_c["$$C:\begin{pmatrix}2 & 4 & 6 \\\ 5 & 7 & 9\end{pmatrix}$$"]
 
 style arr_a stroke-width:0px
 style arr_ab stroke-width:0px
 style arr_b stroke-width:0px
 style arr_c stroke-width:0px

この場合、Aは下にデータを拡張することでBと形状を揃えます。この拡張は上下左右にデータをコピーすることで拡張します。

ではこの場合はどうでしょうか。

③

graph LR
 arr_a["$$A:\begin{pmatrix}1 & 2 & 3 \\\ 1 & 2 & 3\end{pmatrix}$$"]
 arr_a --> Add[Add]
 arr_b["$$B:\begin{pmatrix}1 & 2 \\\ 3 & 4 \\\ 5 & 6\end{pmatrix}$$"] --> Add

 style arr_a stroke-width:0px
 style arr_b stroke-width:0px

この場合、AまたはBをブロードキャストすることはできません。もちろん列数の2,3の最小公倍数を考えれば、AとBが同じような形状に拡張できます。しかし、ブロードキャストは基本的に一方の行列の行または列を整数倍するのみなので、 この場合は対応できず、エラーを吐きます。

TODO: Arrayのブロードキャストのテキスト

このArray型による自動でのブロードキャストは大変便利ですが、一方で暗黙的に形状が変化してしまうため、気づかないうちに計算グラフの中で計算データのズレが生じる恐れがあります。 ①の例を改めてみてみます。

①

graph LR
 arr_a("$$A:\begin{pmatrix}1 \end{pmatrix}$$")  -->|"暗黙的な変換"| arr_ab["$$A':\begin{pmatrix}1 & 1 & 1 \\\ 1 & 1 & 1\end{pmatrix}$$"]
 arr_ab --> Add[Add]
 arr_b["$$B:\begin{pmatrix}1 & 2 & 3 \\\ 4 & 5 & 6\end{pmatrix}$$"] --> Add
 Add --> arr_c["$$C:\begin{pmatrix}2 & 3 & 4 \\\ 5 & 6 & 7\end{pmatrix}$$"]
 
 style arr_a stroke-width:0px
 style arr_ab stroke-width:0px
 style arr_b stroke-width:0px
 style arr_c stroke-width:0px

AからA’への 暗黙的な変換 はArray型が自動で行うものであり、私たちが実装したAdd関数といったRcVariableを扱うFunction構造体はこの変換を認識できません。よってAdd関数の場合、inputのAの行列の形状を(2,3)と認識するため、バックプロパゲーションにおいても(2,3)の微分の値として返してしまいます。しかし、本来のAの形状は(1,1)、もしくはスカラーであり形状がずれています。このズレがバックプロパゲーションで伝搬してしまうのは防がなければなりません。 ではこの暗黙的なブロードキャスト機能をどのように明示的な機能にすればよいでしょうか。それは、AをA’に変換した処理を見逃さず、この処理によって生じた変化をうまく補正するFunction構造体を実装することです。ではこの関数を考えてみます。

BroadcastTo関数の実装

はじめに、ブロードキャストで形状を正しく変更する関数をBroadcastTo とします。 このブロードキャストの機能はArray型のメソッドである broadcast() によって成り立っています。このメソッドには引数として演算する際のもう片方の行列の形状を渡します。そして、その渡した形状が自身の行列の形状とブロードキャストで対応するか確かめ、うまく合う場合は自身の形状を変形させます。うまく合わない場合はエラーを返します。このメソッドを用いて、行列を変形させます。

BroadcastTo関数 を実装する前にひとつ確認しておかなくてはならないことがあります。それはBackward です。ここに二つの例を挙げました。

BroadcastTo関数のForwardとBackward

Forward

graph LR
 arrx["$$x:\begin{pmatrix}2 \end{pmatrix}$$"] --> f[BroadcastTo]
 f --> arry["$$y:\begin{pmatrix}2 & 2 & 2 \\\ 2 & 2 & 2\end{pmatrix}$$"]
 
 style arrx stroke-width:0px
 style arry stroke-width:0px

Backward

graph RL
 arrgy["$$gy:\begin{pmatrix}1 & 1 & 1 \\\ 1 & 1 & 1\end{pmatrix}$$"] --> reshape["SumTo ( = BroadcastTo')"]
 reshape --> arrgx["$$gx:\begin{pmatrix}6\end{pmatrix}$$"] 
 
 style arrgx stroke-width:0px
 style arrgy stroke-width:0px

もう一つの場合

Forward

graph LR
 arrx["$$x:\begin{pmatrix}3 & 3 & 3\end{pmatrix}$$"] --> f[BroadcastTo]
 f --> arry["$$y:\begin{pmatrix}3 & 3 & 3 \\\ 3 & 3 & 3\end{pmatrix}$$"]
 
 style arrx stroke-width:0px
 style arry stroke-width:0px

Backward

graph RL
 arrgy["$$gy:\begin{pmatrix}1 & 1 & 1 \\\ 1 & 1 & 1\end{pmatrix}$$"] --> reshape["SumTo ( = BroadcastTo')"]
 reshape --> arrgx["$$gx:\begin{pmatrix}2 & 2 &2 \end{pmatrix}$$"] 
 
 style arrgx stroke-width:0px
 style arrgy stroke-width:0px

この二つは違う形状の input を同じ形状に変形させています。注目すべき点はgx の行列の値です。BroadcastTo でinputの行列の値をいわばコピーして拡張しています。なので、 Backward の場合はそのコピーを戻す形で足し合わせることで gx に戻します。

ここでBroadcastToで拡張した際、バックプロパゲーションでは要素を足すことを数学的に確かめます。 $$ X:\begin{pmatrix} x_0 & x_1 & x_2 \end{pmatrix} \xrightarrow{\text{BroadCastTo}} Y:\begin{pmatrix} y_0 & y_1 & y_2\\ y_3 & y_4 & y_5 \end{pmatrix} = \begin{pmatrix} x_0 & x_1 & x_2\\ x_0 & x_1 & x_2 \end{pmatrix} $$

$$ \frac{\partial L}{\partial X}:\begin{pmatrix} \frac{\partial L}{\partial y} \cdot \frac{\partial y}{\partial x_0} & \frac{\partial L}{\partial y} \cdot \frac{\partial y}{\partial x_1} & \frac{\partial L}{\partial y} \cdot \frac{\partial y}{\partial x_2} \end{pmatrix} \xleftarrow{\text{SumTo}} \frac{\partial L}{\partial Y}:\begin{pmatrix} \frac{\partial L}{\partial y_0} & \frac{\partial L}{\partial y_1} & \frac{\partial L}{\partial y_2} \\ \frac{\partial L}{\partial y_3} & \frac{\partial L}{\partial y_4} & \frac{\partial L}{\partial y_5} \end{pmatrix} $$ この際 \(x_0\) の偏微分、$$\frac{\partial L}{\partial x_0} = \frac{\partial L}{\partial y} \cdot \frac{\partial y}{\partial x_0}$$で考えてみます。\(y\)での偏微分ですが、この\(y\)は\(x_0\)に関係する\(y\)の要素なので、\(y_0,y_3\)を指します。よってこの部分は正確には $$\frac{\partial L}{\partial x_0} = \frac{\partial L}{\partial y_0} \cdot \frac{\partial y_0}{\partial x_0} + \frac{\partial L}{\partial y_3} \cdot \frac{\partial y_3}{\partial x_0}$$ となります。ここで、\(\frac{\partial y_0}{\partial x_0}\)と\(\frac{\partial y_3}{\partial x_0}\)はただ値をコピーしただけ関係なのでそれぞれ1です。なので最終的に$$\frac{\partial L}{\partial x_0} = \frac{\partial L}{\partial y_0} + \frac{\partial L}{\partial y_3}$$ となります。これは\(gy\)の要素の和を取ったものです。

そして、このバックプロパゲーションを実装するにあたって、gy の要素を正しく足してgx を返す関数が必要となります。これがもう一つの関数、SumTo です。つまり、SumTo関数 を実装しなければ BroadcastTo関数 のバックプロパゲーションが実装できないのです。なので、下にあるBroadcastToのコードは次に説明する SumTo関数 を見てからにしてください。

struct BroadcastTo {
    inputs: Vec<RcVariable>,
    output: Option<Weak<RefCell<Variable>>>,
    shape: IxDyn,
    generation: i32,
    id: usize,
}

impl Function for BroadcastTo {
    fn call(&mut self) -> RcVariable {
        let inputs = &self.inputs;
        if inputs.len() != 1 {
            panic!("BroadcastToは一変数関数です。inputsの個数が一つではありません。")
        }

        let output = self.forward(inputs);

        if get_grad_status() == true {
            //inputのgenerationで一番大きい値をFuncitonのgenerationとする
            self.generation = inputs.iter().map(|input| input.generation()).max().unwrap();

            //  outputを弱参照(downgrade)で覚える
            self.output = Some(output.downgrade());

            let self_f: Rc<RefCell<dyn Function>> = Rc::new(RefCell::new(self.clone()));

            //outputsに自分をcreatorとして覚えさせる
            output.0.borrow_mut().set_creator(self_f.clone());
        }

        output
    }

    fn forward(&self, xs: &[RcVariable]) -> RcVariable {
        let x = &xs[0];

        let y_shape = self.shape.clone();

        // 実際の形状を `IxDynImpl` からスライスとして抽出

        let y_data = x.data().broadcast(y_shape).unwrap().mapv(|x| x.clone());

        y_data.rv()
    }

    fn backward(&self, gy: &RcVariable) -> Vec<RcVariable> {
        let x = &self.inputs[0];
        let x_shape = IxDyn(x.data().shape());

        let gx = sum_to(gy, x_shape);
        let gxs = vec![gx];

        gxs
    }

    fn get_inputs(&self) -> &[RcVariable] {
        &self.inputs
    }

    fn get_output(&self) -> RcVariable {
        let output;
        output = self
            .output
            .as_ref()
            .unwrap()
            .upgrade()
            .as_ref()
            .unwrap()
            .clone();

        RcVariable(output)
    }

    fn get_generation(&self) -> i32 {
        self.generation
    }
    fn get_id(&self) -> usize {
        self.id
    }
}
impl BroadcastTo {
    fn new(inputs: &[RcVariable], shape: IxDyn) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Self {
            inputs: inputs.to_vec(),
            output: None,
            shape: shape,
            generation: 0,
            id: id_generator(),
        }))
    }
}

fn broadcast_to_f(xs: &[RcVariable], shape: IxDyn) -> RcVariable {
    BroadcastTo::new(xs, shape).borrow_mut().call()
}


pub fn broadcast_to(x: &RcVariable, shape: IxDyn) -> RcVariable {
    let y = broadcast_to_f(&[x.clone()], shape);
    y
    let y;
    let x_shape = IxDyn(x.data().shape());
    if x_shape == shape {
        y = x.clone();
    } else {
        y = broadcast_to_f(&[x.clone()], shape);
    }

    y
}

pub fn sum_to(x: &RcVariable, shape: IxDyn) -> RcVariable {
    let y;
    let x_shape = IxDyn(x.data().shape());
    if x_shape == shape {
        y = x.clone();
    } else {
        y = sum_to_f(&[x.clone()], shape);
    }

    y
}

この下の説明は次の SumTo関数 の説明を読んでから進めてください。

引数としてこの形状に変形させたいという形状を渡し、 shape フィールドとして保持します。forwardで broadcast メソッドにその形状を渡すことで、ブロードキャストを実現します。ブロードキャストで返された行列は参照を渡すので、要素ごとにクローンして実体としてのデータを求めます。

SumTo関数も実装したら、 BroadcastTo 関数を二つの場合でテストしてみましょう。

TODO: コード検証必要

#[test]
    fn broadcast_to_test() {
        use crate::core_new::ArrayDToRcVariable;

        let x = array![3.0, 3.0, 3.0].rv();

        let mut y = broadcast_to(&x,Dim(IxDyn(&[2, 3])));

        println!("y = {}", y.data()); 

        y.backward(false);

        println!("x_grad = {:?}", x.grad().unwrap().data()); // 
    }

SumTo関数の実装

先ほど、ブロードキャストで形状を正しく変更する関数をBroadcastTo を説明しました。その過程で、BroadcastTo関数のBackwardにおいて新たな関数、SumTo が必要になりました。ここではこの関数について解説します。

ではこのSumTo関数 のFowwardとBackwardを見てみましょう。すると、もうお気づきかもしれませんが、実は、BroadcastToとSumTo は ForwardとBackwardにおいて表裏一体の関係 なのです。

SumTo関数のForwardとBackward

Forward

graph LR
 arrx["$$x:\begin{pmatrix}2 & 2 & 2 \\\ 2 & 2 & 2\end{pmatrix}$$"] --> f[SumTo]
 f --> arry["$$y:\begin{pmatrix}12 \end{pmatrix}$$"]
 
 style arrx stroke-width:0px
 style arry stroke-width:0px

Backward

graph RL
 arrgy["$$gy:\begin{pmatrix}1\end{pmatrix}$$"] --> reshape["BroadcastTo ( = SumTo')"]
 reshape --> arrgx["$$gx:\begin{pmatrix}1 & 1 & 1 \\\ 1 & 1 & 1\end{pmatrix}$$"] 
 
 style arrgx stroke-width:0px
 style arrgy stroke-width:0px

もう一つの場合

Forward

graph LR
 arrx["$$x:\begin{pmatrix}3 & 3 & 3 \\\ 3 & 3 & 3\end{pmatrix}$$"] --> f[SumTo]
 f --> arry["$$y:\begin{pmatrix}6 & 6 & 6\end{pmatrix}$$"]
 
 style arrx stroke-width:0px
 style arry stroke-width:0px

Backward

graph RL
 arrgy["$$gy:\begin{pmatrix}1 & 1 &1 \end{pmatrix}$$"] --> reshape["BroadcastTo ( = SumTo')"]
 reshape --> arrgx["$$gx:\begin{pmatrix}1 & 1 & 1 \\\ 1 & 1 & 1\end{pmatrix}$$"] 
 
 style arrgx stroke-width:0px
 style arrgy stroke-width:0px

これらのグラフと前のBroadcastToの時のグラフを比べて見ると、 input をブロードキャストで拡張したのを、SumTo で戻し、逆にSumTo で和を取ったものを、BroadcastTo で拡張して戻すというお互いに依存しあう関係になっています。なので、これら二つの関数を実装する際は同時に実装する必要があります。

struct SumTo {
    inputs: Vec<RcVariable>,
    output: Option<Weak<RefCell<Variable>>>,
    shape: IxDyn,
    generation: i32,
    id: usize,
}

impl Function for SumTo {
    fn call(&mut self) -> RcVariable {
        let inputs = &self.inputs;
        if inputs.len() != 1 {
            panic!("SumToは一変数関数です。inputsの個数が一つではありません。")
        }

        let output = self.forward(inputs);

        if get_grad_status() == true {
            //inputのgenerationで一番大きい値をFuncitonのgenerationとする
            self.generation = inputs.iter().map(|input| input.generation()).max().unwrap();

            //  outputを弱参照(downgrade)で覚える
            self.output = Some(output.downgrade());

            let self_f: Rc<RefCell<dyn Function>> = Rc::new(RefCell::new(self.clone()));

            //outputsに自分をcreatorとして覚えさせる
            output.0.borrow_mut().set_creator(self_f.clone());
        }

        output
    }

    fn forward(&self, xs: &[RcVariable]) -> RcVariable {
        let x = &xs[0];
        let y_shape = self.shape.clone();
        let y_data = array_sum_to(&x.data().view(), y_shape);

        y_data.rv()
    }

    fn backward(&self, gy: &RcVariable) -> Vec<RcVariable> {
        let x = &self.inputs[0];

        let x_shape = IxDyn(x.data().shape());

        let gx = broadcast_to(gy, x_shape);
        let gxs = vec![gx];

        gxs
    }

    fn get_inputs(&self) -> &[RcVariable] {
        &self.inputs
    }

    fn get_output(&self) -> RcVariable {
        let output;
        output = self
            .output
            .as_ref()
            .unwrap()
            .upgrade()
            .as_ref()
            .unwrap()
            .clone();

        RcVariable(output)
    }

    fn get_generation(&self) -> i32 {
        self.generation
    }
    fn get_id(&self) -> usize {
        self.id
    }
}
impl SumTo {
    fn new(inputs: &[RcVariable], shape: IxDyn) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Self {
            inputs: inputs.to_vec(),
            output: None,
            shape: shape,
            generation: 0,
            id: id_generator(),
        }))
    }
}

fn array_sum_to(x_array: &ArrayViewD<f32>, shape: IxDyn) -> ArrayD<f32> {
    let x_shape = x_array.shape();

    let mut axes_to_sum = HashSet::new();

    // 合計する軸を特定する
    for i in 0..x_shape.len() {
        if i >= shape.ndim() || x_shape[i] != shape[i] {
            axes_to_sum.insert(i);
        }
    }

    let mut y = x_array.to_owned();

    // HashSetの要素をVecに収集し、ソートして逆順にイテレーションする
    let mut sorted_axes: Vec<_> = axes_to_sum.into_iter().collect();
    sorted_axes.sort_unstable();

    // 特定した軸を合計する
    for &axis in sorted_axes.iter().rev() {
        y = y.sum_axis(Axis(axis)).insert_axis(Axis(axis));
    }

    y
}

fn sum_to_f(xs: &[RcVariable], shape: IxDyn) -> RcVariable {
    SumTo::new(xs, shape).borrow_mut().call()
}

pub fn sum_to(x: &RcVariable, shape: IxDyn) -> RcVariable {
    let y;
    let x_shape = IxDyn(x.data().shape());
    if x_shape == shape {                 // 形状が変化しないならそのまま流す
        y = x.clone();
    } else {
        y = sum_to_f(&[x.clone()], shape);
    }

    y
}

基本的にはBroadcastTo関数と同様に、引数として求める形状を渡します。しかし、違う点はBroadcastTo を実際に処理するメソッド、broadcast() がArray型 には装備されていたのに対して、sumto というメソッドは現状実装されていません。つまり、自分でその処理を書かなければなりません。その処理が array_sum_to() です。この処理については詳しくは解説しませんが、簡単な説明として、inputの行列のある軸に対して和をとる際、求める形状の軸に対してどこが違うのか、すなわちどの軸を対象とすべきかという軸を求め、それを配列として保持し、その軸においてsum_axis() メソッドで和をとります。また、sum_to() 関数では、形状が変化しない場合は、そのまま流すようにします。

BroadcastToも実装したら、 SumTo関数関数を二つの場合でテストしてみましょう。

TODO: コード検証必要

#[test]
    fn sum_to_test() {
        use crate::core_new::ArrayDToRcVariable;

        let x = array![[3.0, 3.0, 3.0],[ 3.0, 3.0,3.0]].rv();

        let mut y = sum_to(&x,Dim(IxDyn(&[1, 3])));

        println!("y = {}", y.data()); 

        y.backward(false);

        println!("x_grad = {:?}", x.grad().unwrap().data()); 
    }

ブロードキャストの四則演算対応

ブロードキャストによる形状の変換をFunctionトレイトの関数として実装することで、この変形を正しく逆伝播することができるようになりました。ではこのブロードキャストの関数を実際に組み込んでいきます。

このブロードキャストははじめに説明した通り、行列の四則演算を簡単に行うためのものでした。つまり、このブロードキャストという機能は四則演算、すなわちAdd,Mul,Sub,Div らの関数を呼び出す際に生じる機能です。なので、この関数の中でブロードキャストが起きるかどうかを調べ、もし起きる場合は私たちが先ほど実装したBroadcastTo 関数を用いてブロードキャストを行えば、形状変換が正しく認識されるのです。

ここではAdd,Mul 関数を例に挙げてどう組み込むのか説明します。

– Forward

graph LR
 arr_a("$$A:\begin{pmatrix}1 \end{pmatrix}$$")  -->|Array型のブロードキャストを認識|arr_ab["$$A':\begin{pmatrix}1 & 1 & 1 \\\ 1 & 1 & 1\end{pmatrix}$$"]
 arr_ab --> Add[Add]
 arr_b["$$B:\begin{pmatrix}1 & 2 & 3 \\\ 4 & 5 & 6\end{pmatrix}$$"] --> Add
 Add --> arr_c["$$C:\begin{pmatrix}2 & 3 & 4 \\\ 5 & 6 & 7\end{pmatrix}$$"]
 
 style arr_a stroke-width:0px
 style arr_ab stroke-width:0px
 style arr_b stroke-width:0px
 style arr_c stroke-width:0px

Backward

TODO:数値ではなく、記号に変更する予定

graph RL

 arr_c["$$C:\begin{pmatrix}2 & 3 & 4 \\\ 5 & 6 & 7\end{pmatrix}$$"] --> Add'
 Add' --> arr_ab["$$A':\begin{pmatrix}1 & 1 & 1 \\\ 1 & 1 & 1\end{pmatrix}$$"]
 arr_ab -->|SumTo関数で元に戻す| arr_a("$$A:\begin{pmatrix}6 \end{pmatrix}$$")
 Add' --> arr_b["$$B:\begin{pmatrix}1 & 2 & 3 \\\ 4 & 5 & 6\end{pmatrix}$$"] 
 arr_b -->|"SumTo関数($$B$$は変形しないためそのまま流す)"| arr_b2("$$B:\begin{pmatrix}1 & 2 & 3 \\\ 4 & 5 & 6\end{pmatrix}$$")


 style arr_a stroke-width:0px
 style arr_ab stroke-width:0px
 style arr_b stroke-width:0px
 style arr_b2 stroke-width:0px
 style arr_c stroke-width:0px

※ここで重要なのはブロードキャストで変形させた形状をどう戻すかであり、数値は関係ありません。

順伝播の方はArray型 のブロードキャスト機能で自動で変形させます。重要なのは、逆伝播の際、inputの形状と微分して出力した行列の形状が異なっているかどうか調べることです。ここでの場合、\(A\)と\(A’\)、また\(B\)ともう一つの\(B\)の形状が違うかどうかです。もし異なるなら、順伝播の際にブロードキャストが起きたと認識することができます。そしてここで SumTo関数に元のinputの形状を渡し、通すことでブロードキャストでの拡張を元に戻すことができます。SumTo には形状が同じ場合はそのまま流すようになっているので、\(B\)の場合はそのまま流します。

Add関数 の変更点

impl Function for AddF {

    fn backward(&self, gy: &RcVariable) -> Vec<RcVariable> {
        let mut gx0 = gy.clone();
        let mut gx1 = gy.clone();

        let x0 = &self.inputs[0];
        let x1 = &self.inputs[1];

        let x0_shape = IxDyn(x0.data().shape()); // <- inputの形状を取得
        let x1_shape = IxDyn(x1.data().shape()); // <- 

        if x0_shape != x1_shape {
            gx0 = sum_to(&gx0, x0_shape); // <- SumTo関数に通す。
            gx1 = sum_to(&gx1, x1_shape); // 形状が同じならそのまま流す。
        }

        let gxs = vec![gx0, gx1];

        gxs
    }
}

Mul関数 の変更点

impl Function for MulF {

    fn backward(&self, gy: &RcVariable) -> Vec<RcVariable> {
        let x0 = &self.inputs[0];
        let x1 = &self.inputs[1];

        let mut gx0 = x1.clone() * gy.clone();
        let mut gx1 = x0.clone() * gy.clone();

        let x0_shape = IxDyn(x0.data().shape());
        let x1_shape = IxDyn(x1.data().shape());

        if x0_shape != x1_shape {
            gx0 = sum_to(&gx0, x0_shape);
            gx1 = sum_to(&gx1, x1_shape);
        }

        let gxs = vec![gx0, gx1];
        gxs
    }
}

関数の変更は backward() メソッドの中で sum_to() の処理を追加することです。ただし、Mul の場合はAdd とは異なり、微分の値を計算するので、\(gx_0,gx_1\) を計算してから sum_to に通します。

ではこのように他の四則演算の関数も変更します。

ではブロードキャストに対応した Add 関数をテストしてみましょう。

#[test]
    fn add_with_broadcast_test() {
        use crate::core_new::ArrayDToRcVariable;

        let a = array![1.0, 1.0, 1.0, 1.0, 1.0].rv();

        let b = array![2.0].rv();

        let mut c = a.clone() + b.clone();

        println!("c = {}", c.data()); // [3,3,3,3,3]

        c.backward(false);

        println!("a_grad = {:?}", a.grad().unwrap().data()); // [1.0,1.0,1.0,1.0,1.0]
        println!("b_grad = {:?}", b.grad().unwrap().data()); // [5.0]
    }

パッケージ化

私たちは今まで様々な関数を実装してきました。これらのたくさんの機能をまとめ、整ったフレームワークにするために、複数のファイルに分割しましょう。

ライブラリクレートの作成

私たちが開発するフレームワークを外部から読み込むライブラリクレートとして実装します。cargo new - lib で作成しましょう。今回はフレームワークの名前からstucrs とします。

次にこのようなディレクトリ構成にします。

alt text
※ 一部ファイルを省略しています。

今まで実装してきたmain.rs ファイルはtestフォルダ の中にありますが、それとは別にstucrs ライブラリのためをクレートを先ほどの説明で用意します。その後、stucrsの中のsrcにfunctionsフォルダ、core.rsファイルを入れます。また、functionsフォルダに math.rs, matrix.rs, mod.rs ファイルを入れます。

ではこれらの追加したフォルダ、ファイルについて説明します。functions フォルダは今まで実装してきたFunction構造体である関数を入れるところであり、たくさんの種類があるため、種類ごとに関数を分けて保存します。今のところ、math.rs, matrix.rs を用意していますが、今後も関数を実装していくにつれ、新たなファイルが必要になるため、追加していきます。mod.rs は追加していったファイル(math.rs, matrix.rs)などを管理するためのファイルです。
core.rs ファイルはフレームワークの核となる機能を保存します。具体的にはVariable、RcVariable,Functionトレイト などです。なお、四則演算の関数はFunction構造体ですが、基本的な関数なので、core.rs に入れます。

core.rsのファイルは現在のリポジトリではcore_new.rs という名前になっています。

また、lib.rs にはRcVariable にオーバーロードするコードを移します。

lib.rs の中身

use core_new::RcVariable;
use core_new::{add, div, mul, neg, sub};
use std::ops::{Add, Div, Mul, Neg, Sub};

//演算子のオーバーロード

impl Add for RcVariable {
    type Output = RcVariable;
    fn add(self, rhs: RcVariable) -> Self::Output {
        // add_op関数はRc<RefCell<Variable>>を扱う
        let add_y = add(&[self.clone(), rhs.clone()]);
        add_y
    }
}

impl Mul for RcVariable {
    type Output = RcVariable;
    fn mul(self, rhs: RcVariable) -> Self::Output {
        let mul_y = mul(&[self.clone(), rhs.clone()]);
        mul_y
    }
}

impl Sub for RcVariable {
    type Output = RcVariable;
    fn sub(self, rhs: RcVariable) -> Self::Output {
        let sub_y = sub(&[self.clone(), rhs.clone()]);
        sub_y
    }
}

impl Div for RcVariable {
    type Output = RcVariable;
    fn div(self, rhs: RcVariable) -> Self::Output {
        let div_y = div(&[self.clone(), rhs.clone()]);
        div_y
    }
}

impl Neg for RcVariable {
    type Output = RcVariable;
    fn neg(self) -> Self::Output {
        let neg_y = neg(&[self.clone()]);
        neg_y
    }
}


pub mod core;
pub mod functions;  
　　　　　　　　　　　// <- 今後新しいファイルができたらここに追加

mod.rs の中身

pub mod math;
pub mod matrix; 
　　　　　　　　　　// <- 今後新しいファイルができたらここに追加

基本的には今まで実装してきた関数のコードをそのままファイルに貼り付けるだけですが、一つだけ変更しなければならない点があります。

pubをつける

今まで実装して関数を外部クレートとして読み込むようになるため、一部の関数を外からアクセスできるようpub を付けます。これはユーザーがライブラリのあらゆる関数にアクセスできるのを防ぐために、外からアクセスできるものを設定するためのものです。

色々説明してきましたが、移行の仕方は実際のコードを見るのが一番手っ取り早いです。ここからリポジトリのstucrsを見てどう移すか、またどの関数にpubを付けるのかを確認してください。

手動によるニューラルネットワークの構築

ではいよいよ今まで造ってきた関数を組み合わせてニューラルネットワークを構築していきます。

私たちはニューラルネットワークを自動で構築するフレームワークを実装していくわけですが、最初は手動で組み立て、どのように自動化していくか考えていきます。

はじめにニューラルネットワークの構造を簡単に説明します。深層学習の基本的な原理は全結合層 による線型変換、活性化関数 による非線型変換 を互いにパイのように重ね、データを線型変換→非線型変換→線型変換→非線型変換の順番で次々と変換する処理です。

graph LR
 x((Input)) --> A[線型変換]
 A --> B[非線型変換]
 B --> C[線型変換]
 C --> D[非線型変換]
 D --> E[...]

この二つの種類の変換を深い層にし、多く用いることで多様な表現を行う、すなわち複雑な問題にも対応できるようになります。この深い層を深層学習(DeepLearning)と言います。

なぜこの2種類の変換を用いるかというと、それは線型、非線型の2種類のデータに対応できるようになるからです。例えば直線 \(y = ax\) は線型のデータです。一方、 \(y = sinx\) は非線型のデータです。このように多様なデータに対応するために、線型、非線型を組み合わせるのです。

では次に線型変換、非線型変換について説明します。

線型変換は以前に実装した行列の積 を基本とする処理です。MatMul の説明の際にも触れましたが、inputのデータの\(X\)と重みの\(W\) を行列積でとることで、ニューラルネットワークの全結合層を表すことができます。この時の数式は \(Y = X \cdot W + b\) となります。この\(b\) はバイアスと言います。これは先ほどの直線の式 \(y = ax+b\)と似ています。スカラー版の直線の式を行列に拡張したというイメージです。一方非線型変換はシグモイド関数や\(tanhx\) などといった活性化関数に使われます。

機械学習の分野における線型性と非線型性の違いはよく、直線となめらかな曲線で説明されます。線型性のみの場合は、分類において直線による線引きしかできませんが、非線型性を持つと、その線を曲線として表現でき、より複雑な分類を行えるようになります。

次に誤差逆伝播法を実装しますが、説明は後の誤差逆伝播法のページで行います。

ではこれらの線型変換、非線型変換、そして誤差逆伝播法のための関数を簡易的に実装し、手動でニューラルネットワークを構築していきます。

Linear関数

では線型変換を行うLinear関数を実装します。線型変換は前に確認した通り、行列積を用いた\(Y = X \cdot W + b\)という式になります。正確には線型変換は\(Y = X \cdot W \)の行列積のみを指し、バイアスの\(b\)を加えると非線型変換になってしまいますが、深層学習の分野ではこの式を線型変換と捉えるのが一般的です。またこのような式は深層学習のフレームワークにおける全結合層の処理を指します。

この式をグラフにするとこのようになります。

graph LR
 x((X)) --> matmul[MatMul]
 w((W)) --> matmul
 matmul --> t((t))
 t --> add[Add]
 b((b)) --> add
 add --> y((Y))

バイアスのbは省略されることもあります。 私たちはこれらの関数(Matmul,Add)を実装してきたので、この式のバックプロパゲーションを自動で行うことができます。ではこの式を関数として定義します。

pub fn linear_simple(x: &RcVariable, w: &RcVariable, b: &Option<RcVariable>) -> RcVariable {
    let t = matmul(&x, &w);

    let y;
    if let Some(b_rc) = b {
        y = t + b_rc.clone();
    } else {
        y = t;
    }

    y
}

先ほどの式をそのままにして実装しました。このようにシンプルで可読性の高いコードを書くことができるのは、私たちが今までの間、つながりを作り、バックプロパゲーションを行う構造を自動化し、それをオーバーロードしたことによるものです。このコードは、xとwの行列積をとってtを出力し、バイアスはオプション的な存在なので、必要か必要でないかで場合分けします。使用しない場合、tをそのまま流します。

今後私たちは主に2種類の関数を実装していきますが、役割が全く異なり、混同してしまうと危険なのでここで整理しておきます。

今私たちが作成した関数は処理をまとめたものです。つまり、Linearの処理を行うために、MatmulやAdd といったFunction構造体の関数を呼び出すようにまとめたものです。もっと分かりやすく言うならば、このlinear_simple関数はただFunction構造体順序正しく呼び出しているだけであり、実際の計算をしているのは今まで実装してきたFunction構造体ということです。このただ呼び出す関数とFunction構造体の関数の違いを理解しておくことは、今後のニューラルネットワークの構築において重要になってきます。これは今後のLayerの実装のところで理解できるので、心配せず進んでください。

活性化関数

続いて非線型変換の関数です。こちらは先ほどの線形変換とは逆の非線型の変換を行う関数です。深層学習のフレームワークではこのような関数は活性化関数(activation functions)と呼ばれています。ではその中でも有名はシグモイド関数を紹介し、実装します。

Sigmoid関数

TODO:グラフ追加

$$y = \frac{1}{1 + \exp(-x)}$$

式はこのようにネイピア数の指数関数を含んでいて、少し複雑な関数ですが、私たちは分数、すなわち割り算、Exp で指数にも対応していますので、このように直感的にコードを書くことができます。

pub fn sigmoid_simple(x: &RcVariable) -> RcVariable {
    let mainasu_x = -x.clone();
    let y = 1.0f32.rv() / (1.0f32.rv() + exp(&mainasu_x));
    y
}

活性化関数は他にもTanhやReLUなど様々な種類があり、それぞれ学習のスタイルに対して向き不向きがあります。自分の実行する学習に合わせた活性化関数を選ぶ必要があります。他の活性化関数も同様に実装していき、選択の多い実用的なフレームワークを目指します。

誤差逆伝播法

最後に誤差逆伝播法についてです。ニューラルネットワークは誤差逆伝播法を用いて自動で重みを調整するのですが、今までバックプロパゲーションを実装してい来たのはまさにこの誤差逆伝播法を実現するためのものです。

ニューラルネットワークの学習とはいわば関数の最適化です。その時の関数はモデルの出力、すなわちそのモデルの予測値と答えとなる教師データの 誤差を求める関数 です。

誤差が小さいモデルを目指すわけですからその関数の値が小さくなるよう重みである パラメーター を調節します。ではどのようにパラメーターを調整するのでしょうか。ここで、変数の微分値を用います。

パラメーターの更新

私たちが今まで実装してきた変数である RcVariable は grad としてその変数の微分の値を保持しています。最終的な変数は誤差を示す値なので L とすると、あるパラメーター W の grad は\(\frac{\partial L}{\partial W}\) となります。これは L をW の関数と見立てたときの勾配と言えます。この勾配を用いれば、L を最小化するW を求めることできます。そして、その値に近づくように更新するのをすべてのパラメーターで行うことが、深層学習の学習です。この最適化するアルゴリズムは最適化関数によるものであり、いくつか種類がありますが、単に勾配の値の負の方、つまり逆の方に値を近づけていくことで誤差を最小にしていくアルゴリズムを 勾配降下法 と言います。今回はこの一番基本的でシンプルな勾配降下法を使用してみます。勾配降下法のパラメーター更新の方法を数式にすると、$$W \xleftarrow{} W - \alpha\cdot \frac{\partial L}{\partial W}$$となります。\(\alpha\)は学習を表しています。 数式で表すと少し難しく感じますが、プログラムで書くと、

let current_grad_w = W.grad().unwrap().data();
W.0.borrow_mut().data = w_data - lr * current_grad_w;

と理解しやすくなると思います。

この最適化関数についてはその後のOptimizer のところで再び触れます。

ニューラルネットワークのイメージ loss変数からバックプロパゲーションを行うことで、すべての変数の勾配が求まる。それを用いてパラメーターを更新していく。

graph LR
 X((Input)) --> linear1
 subgraph レイヤー1
 w1((W)) --> linear1
 b1((b)) --> linear1
 linear1 --> y1((y))
 end
 
 subgraph レイヤー2
 Sigmoid --> y2((y))
 end
 subgraph レイヤー3
 w2((W)) --> linear2
 b2((b)) --> linear2
 end

y1 --> Sigmoid
y2 --> linear2
linear2 --> pre((予測値))
pre --> Loss
t((答え)) --> Loss
Loss --> loss((L))

誤差を求める関数

続いて先ほど触れた誤差を求める関数です。前に説明した学習時のモデルの予測値と、答えとなる教師データの正解ラベルの誤差を求める関数を 損失関数 と言います。誤差を求めると言いますが、予測値と正解ラベルがどれくらい違うかを評価する誤差を導き出す方法は様々あります。ここでは代表的な関数である 二乗平均誤差 について説明します。
この関数の処理は名前の通り、モデルの予測値と答えの差を2乗し、足し合わせて平均を取るというものです。\(N\)はバッチ数、つまりデータ数を表します。予測値と答えの差が大きいと、\(L\) は大きく、逆に差が小さくなるほど\(L\)は小さくなります。

$$L = \frac{1}{N}\sum_{i=1}^N (y_i-t_i)^2$$

ではこの関数を実装してみます。

fn mean_squared_error(x0:&RcVariable,x1:&RcVariable) ->RcVariable{
    let diff =x0.clone()-x1.clone();
    let len = diff.len() as f32;
    println!("len = {:?}",len);
    
    let error = sum(&diff.pow(2.0), None) /len.rv();
    
    error
}

len はデータ数の\(N\) を表しています。ここで予測値と正解ラベルの行列は同じ形状であるため(答えはスカラーなのに、予測値が行列だというように、答えと形式が異なるなら、もともとニューラルネットワークの構造が正しくありません。)

$$ Y:\begin{pmatrix} y_0 \\ y_1 \\ y_2\\ \vdots\\ y_n \end{pmatrix} - T:\begin{pmatrix} t_0 \\ t_1 \\ t_2\\ \vdots\\ t_n \end{pmatrix} = diff:\begin{pmatrix} y_0 - t_0 \\ y_1 - t_1 \\ y_2 - t_2\\ \vdots\\ y_n - t_n \end{pmatrix} $$

$$ diff:\begin{pmatrix} y_0 - t_0 \\ y_1 - t_1 \\ y_2 - t_2\\ \vdots\\ y_n - t_n \end{pmatrix} \xrightarrow{\text{Square}} \begin{pmatrix} (y_0 - t_0)^2 \\ (y_1 - t_1)^2 \\ (y_2 - t_2)^2\\ \vdots\\ (y_n - t_n)^2 \end{pmatrix} \xrightarrow{\text{Sum}} error:\begin{pmatrix} L \end{pmatrix} $$ このように今までの行列の関数を用いて二乗平均誤差 を求めることができます。

以上でニューラルネットワークを構築する最低限の材料がそろいました。では次のページでニューラルネットワークを構築します。

手動による学習

では手動でニューラルネットワークを構築していきます。使用するのは先ほどまで実装してきたLinear関数、活性化関数、二乗平均誤差関数です。これらを組み合わせ、一つの関数としてバックプロパゲーションを行い、勾配を用いてパラメーターを更新していきます。

学習の準備

では学習するに当たってデータを用意します。今回は非線型性のデータの代表的な\(sin\)関数を用意します。ここでは\(y = sin(2\pi\cdot x) + b\)とします。\(b\)はランダムな値で、\(y\)のデータにノイズを与えます。 

use ndarray::{array, Array, IxDyn};
use std::array;
use std::time::Instant;
use stucrs::core_new::{F32ToRcVariable, RcVariable};
use stucrs::functions_new::{linear_simple,  sigmoid_simple};

use std::f32::consts::PI;

use stucrs::core_new::ArrayDToRcVariable;
use stucrs::functions_new::mean_squared_error;

use ndarray_rand::rand_distr::{StandardNormal, Uniform};
use ndarray_rand::RandomExt;

fn main() {
    let x_data = Array::random((100, 1), Uniform::new(0.0f32, 1.0));
    //let x2 = array![[11.0f32, 12.0, 13.0], [14.0, 15.0, 16.0]].rv();

    let y_data =
        x_data.mapv(|x| (2.0 * PI * x).sin()) + Array::random((100, 1), Uniform::new(0.0f32, 1.0));

    let x = x_data.rv();

    let y = y_data.rv();

    let I = 1;
    let H = 10;
    let O = 1;

    let mut W1 = (0.01f32 * Array::random((I, H), StandardNormal)).rv();

    let mut B1 = Array::zeros(H).rv();

    let mut W2 = (0.01f32 * Array::random((H, O), StandardNormal)).rv();

    let mut B2 = Array::zeros(O).rv();

    fn predict(
        x: &RcVariable,
        w1: &RcVariable,
        b1: &Option<RcVariable>,
        w2: &RcVariable,
        b2: &Option<RcVariable>,
    ) -> RcVariable {
        let t0 = linear_simple(&x, &w1, &b1);
        let t1 = sigmoid_simple(&t0);
        let y = linear_simple(&t1, &w2, &b2);
        y
    }

    //let y_data = 2.0*x_data.clone() +5.0 + Array::random((100, 1), Uniform::new(0.0f32, 1.0));
    //println!("y_data = {:?}",y_data.clone());

    //let x=x_data.rv();

    //let y=y_data.rv();

    //let mut w = Array::zeros((1,1)).rv();
    //let mut b = Array::zeros(1).rv();
    let lr = 0.2;
    let iters = 1000;
    for i in 0..iters {
        let y_pred = predict(&x, &W1, &Some(B1.clone()), &W2, &Some(B2.clone()));

        

        let mut loss = mean_squared_error(&y, &y_pred);

       

        W1.cleargrad();
        W2.cleargrad();
        B1.cleargrad();
        B2.cleargrad();

        loss.backward();
        

        let w1_data = W1.data();
        let b1_data = B1.data();
        let w2_data = W2.data();
        let b2_data = B2.data();

        let current_grad_w1 = W1.grad().unwrap().data();
        let current_grad_b1 = B1.grad().unwrap().data();
        let current_grad_w2 = W2.grad().unwrap().data();
        let current_grad_b2 = B2.grad().unwrap().data();

        

        W1.0.borrow_mut().data = w1_data - lr * current_grad_w1;
        B1.0.borrow_mut().data = b1_data - lr * current_grad_b1;
        W2.0.borrow_mut().data = w2_data - lr * current_grad_w2;
        B2.0.borrow_mut().data = b2_data - lr * current_grad_b2;
        
        if i % 100 == 0 {
            println!("i = {:?}", i);
            println!("loss = {:?}\n", loss.data());
        }
    }
}

はじめにこのニューラルネットワークで用いる変数(パラメーター)を用意します。\(x\)は学習データ、\(y\)は正解ラベル(xとyのペアのデータが教師データです。)、\(w\)と\(b\)は重み(パラメーター)です。この時パラメーターは乱数を用いて最初の重みを設定します。行列でランダムな値を持つ行列を作成するためにArray型で提供されるndarray-randを使用します。これは依存関係を別途追加しておく必要がありますが、これによりArrayにランダムな値を生成する機能が自動で追加されます。

パラメーターの初期値は0～1の範囲で乱数を用いて設定します。もし初期値をある値(0.0など)に設定すると、すべてのパラメーターが画一的なふるまいをしてしまいます。これでは深層学習の複雑な分類は行えません。

TODO:データ、エポック数の説明追加

次に関数を構築します。Linearと活性化関数のsigmoidを通すだけで、自動でモデルとなる関数が生まれます。それを正解ラベルとともにmean_squared_error関数に流すことで、誤差を求める関数となります。この関数を predict() 関数としてまとめ、推論を行います。

predictで出力したモデルの予測値\(y\)と正解ラベル\(y_pred\)をmean_squared_error関数に渡し、誤差をlossとして求めます。その後、この誤差の変数である loss を起点にバックプロパゲーションを行います。それが loss.backward() です。これにより各パラメーターの勾配を導出します。

最後に求めたパラメーターの勾配をもとに各パラメーターを一つずつ更新します。今回は勾配降下法を用います。更新方法は前の誤差逆伝播法で確認してください。

以上より手動による学習を行うことができます。最後のコードにある println!(“loss = {:?}\n”, loss.data()); より、学習する度に loss すなわち予測値と正解との誤差が小さくなるのがわかります。

線形、非線型変換を繰り返し行うことで、sin関数といった非線型な曲線に対応することができました。

今の状態では、パラメーターの初期化、更新を自分たちで行わなければなりません。これがものすごく深い層となった場合、私たちの手に負えなくなります。これから先は、学習にあたり今までの処理を自動化し、フレームワークとしての機能を加えていきます。

Layerの実装

では先ほどの手動でニューラルネットワークの学習を一つずつ自動化していきます。最初は Layer についてです。Layerは主に二つの役割があります。一つ目はデータを計算、処理すること。二つ目はパラメーターを管理することです。

先のコードの関数である Linear を見てみます。するとこの関数は\(w\)と\(b\)の二つのパラメーターを使用します。そして、このパラメーターを外で一括で更新しています。(current_gradなどを使っているところです。)このLayerというものは簡単言えば、関数と関数に関するパラメーターを保持、管理させるためのものです。そうすれば、自分たちが行ってきたパラメーターの初期化や、更新をこのLayerに丸投げすることができます。ではそれを可能にする構造体を実装していきます。

pub fn linear_simple(x: &RcVariable, w: &RcVariable, b: &Option<RcVariable>) -> RcVariable {
    let t = matmul(&x, &w);

    let y;
    if let Some(b_rc) = b {
        y = t + b_rc.clone();
    } else {
        y = t;
    }
    y
}

Layerトレイト

はじめに Layerトレイト を実装します。Layerと言っても様々種類があるので、それらにLayerとしての共通のふるまいを持たせます。考え方は Functionトレイト と同じです。トレイトをもとに、いくつかの構造体を実装していきます。

では core.rs ファイルと同じ階層に layer.rs ファイルを追加します。このファイルは Functions と同じように、Layerをまとめるためのファイルです。新たにファイルを追加したので、モジュールとして認識してもらうよう、 lib.rs、mod.rs に layer.rs の名前を追加しておきます。これに関してはパッケージ化で説明してあります。

では Layerトレイト を layer.rs に追加します。

pub trait Layer: Debug {
    fn set_params(&mut self, param: &RcVariable);
    fn get_input(&self) -> RcVariable;
    fn get_output(&self) -> RcVariable;
    fn call(&mut self, input: &RcVariable) -> RcVariable;
    fn get_generation(&self) -> i32;
    fn get_id(&self) -> usize;
    fn params(&mut self) -> &mut FxHashMap<usize, RcVariable>;
    fn cleargrad(&mut self);
    fn has_params(&self) -> bool;
}

今後、Function構造体 と同じように、 Layerトレイト を継承した構造体を Layer構造体 と呼ぶことにします。

get_input()、get_output()、call()、get_generation()、get_id() に関しては Function構造体 と同様です。ただし、callにおいてレイヤー構造体はinputもoutputも弱参照で保持します。これは単に関数も互いに参照しているため、複雑になってしまうからです。

基本的に、Layer構造体 は Function構造体 にパラメーターを管理させるよう拡張させたものなので、 Functionトレイト と近しいものとなります。このトレイトについては、次に実装する Linear構造体 を参考にして理解していただきたいです。

Layerトレイト の大きな違いはパラメーターの扱いにあります。set_params()、params()、cleargrad()、has_params()　はパラメーターを管理するLayer特有の関数です。まず前提として、Layer構造体はフィールドとして パラメーターの配列を保持しています。 この配列については Linearレイヤー で説明します。

  • set_params()
    この関数はレイヤー構造体にパラメーターを持たせる関数です。引数として RcVariable を渡し、フィールドであるパラメーターの配列に追加します。この時、RcVariable(正確にはRc型の参照)とそのid二つを持たせます。idを持たせる理由は今後の拡張を踏まえ、より管理しやすくするためです。

  • params()
    これは単にフィールドとして保持している配列を返す関数です。これは後のOptimizerでOptimizer構造体にパラメーターを渡して更新してもらう際に使用します。

  • cleargrad()
    これは各パラメーターの勾配、 grad の値を None にする関数です。保持するパラメーターをすべてイテレーターとして取り出し、Noneにします。

  • has_params()
    これはレイヤー構造体がパラメーターを保持しているかを bool型 で返す関数です。実はこの後、パラメーターを保持しないレイヤー構造体を実装します。(もしくは今後の拡張で保持しないレイヤーを実装する可能性があります。)なので、そのための設定です。Linearの場合はもちろん保持するのでtrueを返します。これは個別に変更します。

ではLinearレイヤーを踏まえて説明します。

Linearレイヤー

下はLinearレイヤー構造体のコードです。まずはじめにパラメーターの配列について説明します。これはレイヤーがフィールドとして持つもので、フィールド名は基本的に params です。型は FxHashMap です。この型はいわゆる辞書型というもので、パラメーターの RcVariable とidをペアとして保存します。これは普通の HashMap とは異なりますが、それは、イテレーターとして返す際、HashMapはネットワークの安全性を保つためにあるアルゴリズムを処理して返すのに対し、FxHashMapはただ単に順番通りに返すため、処理が少し速いとされています。今回はパラメーターをただ返す単純な処理で十分なので、こちらの型を用意します。外部クレートなので、各自依存関係を追加してください。

FxHashMapの仕様はHashMapとほぼ同じです。詳しくはFxHashMapをご覧ください。こちらはリポジトリの参考文献に載せているものです。

では Linearレイヤー をもとにレイヤーの計算処理について説明していきます。これに関しては説明が長くなるので、分割して説明していきます。

レイヤー構造体はFunction構造体に対して、パラメーターの初期化、管理という重要な機能があります。まずはパラメーターの初期化についてです。コンストラクタとして構造体を作成し、値を流すという処理の中に、自動でパラメーターの作成し、保持するという処理を加えます。ではパラメーターの初期化の方法を改めて考えます。

手動による学習のところで最初にパラメーターを作成しました。これをもとに考えていくわけですが、行列を作成するにあたって重要なことは形状です。コードを見るとI,H,Oと私たちが形状のデータを与えていますが、今回は自動で行うため、引数として形状の値を渡します。そして行列をnew()またはforward()中で作成します。パラメーターは正規分布に従ってランダムな値にするため、standardNormal といったオプションを用います。((1.0f32 / i_f32).sqrt()) はスケールを合わせる計算です。

また行列計算をするにあたり、列数と行数が一致しているという条件があるので、それを用いれば、ユーザーが行数を渡さなくても、自動で行列を作成することができます。これがforward()のなかで処理させることです。

そして最後にforward()で関数の計算を行います。この処理は普通のFunction構造体の処理です。

TODO: パラメーターのスケール設定
TODO: 説明増やす。特に行列の初期化、forward

pub struct Linear {
    input: Option<Weak<RefCell<Variable>>>,
    output: Option<Weak<RefCell<Variable>>>,
    out_size: u32,
    w_id: Option<usize>,
    b_id: Option<usize>,
    params: FxHashMap<usize, RcVariable>,
    generation: i32,
    id: usize,
}

impl Layer for Linear {
    fn set_params(&mut self, param: &RcVariable) {
        self.params.insert(param.id(), param.clone());
    }
    fn get_input(&self) -> RcVariable {
        let input = self
            .input
            .as_ref()
            .unwrap()
            .upgrade()
            .as_ref()
            .unwrap()
            .clone();
        RcVariable(input)
    }

    fn get_output(&self) -> RcVariable {
        let output;
        output = self
            .output
            .as_ref()
            .unwrap()
            .upgrade()
            .as_ref()
            .unwrap()
            .clone();

        RcVariable(output)
    }

    fn call(&mut self, input: &RcVariable) -> RcVariable {
        let output = self.forward(input);
        self.input = Some(input.downgrade());
        self.output = Some(output.downgrade());

        output
    }

    fn get_generation(&self) -> i32 {
        self.generation
    }
    fn get_id(&self) -> usize {
        self.id
    }
    fn params(&mut self) -> &mut FxHashMap<usize, RcVariable> {
        &mut self.params
    }

    fn cleargrad(&mut self) {
        for (_id, param) in self.params.iter_mut() {
            param.cleargrad();
        }
    }

    fn has_params(&self) -> bool {
        true
    }
}

impl Linear {
    fn forward(&mut self, x: &RcVariable) -> RcVariable {
        if let None = &self.w_id {
            let i = x.data().shape()[1];
            let o = self.out_size as usize;
            let i_f32 = i as f32;

            let w_data: ArrayBase<OwnedRepr<f32>, Dim<[usize; 2]>> =
                &Array::random((i, o), StandardNormal) * ((1.0f32 / i_f32).sqrt());

            let w = w_data.rv();

            self.w_id = Some(w.id());
            self.set_params(&w.clone());
        }

        // フィールドでパラメータのidを保持しているので、idでパラメータを呼び出す
        let w_id = self.w_id.unwrap();
        let w = self.params.get(&w_id).unwrap();

        //bはoption型なので、場合分け
        let b;
        if let Some(b_id_data) = self.b_id {
            b = self.params.get(&b_id_data).cloned();
        } else {
            b = None;
        }

        let y = linear_simple(&x, &w, &b);

        y
    }

    pub fn new(out_size: u32, biased: bool, opt_in_size: Option<u32>) -> Self {
        let mut linear = Self {
            input: None,
            output: None,
            out_size: out_size,
            w_id: None,
            b_id: None,
            params: FxHashMap::default(),
            generation: 0,
            id: id_generator(),
        };

        //in_sizeが設定されていたら、ここでWを作成
        //されていない場合は後で作成
        if let Some(in_size) = opt_in_size {
            let i = in_size as usize;
            let o = out_size as usize;

            let i_f32 = in_size as f32;

            let w_data: ArrayBase<OwnedRepr<f32>, Dim<[usize; 2]>> =
                &Array::random((i, o), StandardNormal) * ((1.0f32 / i_f32).sqrt());

            let w = w_data.rv();

            linear.w_id = Some(w.id());
            linear.set_params(&w.clone());
        }

        if biased == true {
            let b = Array::zeros(out_size as usize).rv();
            linear.b_id = Some(b.id());
            linear.set_params(&b.clone());
        }

        linear
    }
}

レイヤー構造体が行うのは実質パラメーターの初期化し、保持しながらFunction構造体に渡すことだけです。グラフのように レイヤーは関数をあくまで包むイメージ です。なので順伝播、バックプロパゲーションといった計算はFunction構造体の関数に任せる構造となっています。

graph LR
 X((Input)) --> linear1
 subgraph レイヤー1
 w1((W)) --> linear1
 b1((b)) --> linear1
 linear1 --> y1((y))
 end
 
 subgraph レイヤー2
 Sigmoid --> y2((y))
 end
 subgraph レイヤー3
 w2((W)) --> linear2
 b2((b)) --> linear2
 end

y1 --> Sigmoid
y2 --> linear2
linear2 --> pre((予測値))
pre --> Loss
t((答え)) --> Loss
Loss --> loss((L))

Layer構造体も Function構造体 のように今後様々なレイヤーを実装していきます。その際は、このLinear構造体を参考にして実装してください。例えばLinearに活性化関数をセットした Denseレイヤー、活性化関数のみを扱う activationレイヤー などです。これらのレイヤーの実装に関しては補足のところで説明したいと思います。

TODO: 補足説明追加

Modelの実装

続いて Model の実装です。前のLayerがパラメーターを管理するなら、このModelはレイヤーを管理するものです。このModelも同様に、Modelトレイト を実装し、様々な構造体に継承させるので、Modelトレイトを実装した構造体をModel構造体 と呼ぶことにします。

先ほどと同じように core.rs ファイルと同じ階層に model.rs ファイルを追加します。モジュールとして認識してもらうよう、 lib.rs、mod.rs に model.rs の名前を追加しておきます。

pub trait Model {
    fn stack(&mut self, layer: impl Layer + 'static);
    fn forward(&mut self, x: &RcVariable) -> RcVariable;
    fn layers(&self) -> Rc<RefCell<Vec<Box<dyn Layer + 'static>>>>;
    fn layers_mut(&mut self) -> Rc<RefCell<Vec<Box<dyn Layer + 'static>>>>;
    fn cleargrad(&mut self);
}

Modelトレイトを実装しましたが、考え方は先ほどのLayerと似ています。先ほどと同じように基本となるModel構造体の BaseModel を交えて説明します。

graph LR

 X((Input)) --> layer1
 subgraph Model
 layer1 --> activation1[活性化関数]
 activation1 --> layer2
 layer2 --> activation2[活性化関数]
 activation2 --> ...
 end

 ... --> pre((予測値))
pre --> Loss
t((答え)) --> Loss
Loss --> loss((L))

ModelがフィールドとしてLayerの配列を保持します。これは先ほどのLayerがパラメーターを配列で保持するのと同じことです。この際、フィールドの配列にLayerを追加する処理を stack() としています。そのように設定すれば、下のコードのようにmodelにレイヤーを持たせることができます。次に forward() で保持しているLayerたちをイテレーターで取り出し、データを順に流していきます。最後に cleargrad() でLayerを取り出し、Layerが保持するパラメーターの微分を初期化することで、すべてのパラメーターに関して微分を初期化できます。ModelとしてLayerを扱うことで、手動によるパラメーターの管理をすべてModelに任せることができました。

let mut model = BaseModel::new();
model.stack(L::Linear::new(10, true, None));

また、私たちはニューラルネットワークの層をブロックと見立て、それをモデルに積んでいくイメージを持てるような設計にしました。 model.stack() によってLayerというブロックを一つずつ積んでいく感覚を、今後より大規模なニューラルネットワークを構築する際に感じることができると思います。

pub struct BaseModel {
    input: Option<Weak<RefCell<Variable>>>,
    output: Option<Weak<RefCell<Variable>>>,
    layers: Rc<RefCell<Vec<Box<dyn Layer + 'static>>>>,
}

impl Model for BaseModel {
    fn stack(&mut self, layer: impl Layer + 'static) {
        let add_layer = Box::new(layer);
        self.layers.borrow_mut().push(add_layer);
    }

    fn cleargrad(&mut self) {
        for layer in self.layers.borrow_mut().iter_mut() {
            layer.cleargrad();
        }
    }
    fn layers(&self) -> Rc<RefCell<Vec<Box<dyn Layer + 'static>>>> {
        self.layers.clone()
    }

    fn layers_mut(&mut self) -> Rc<RefCell<Vec<Box<dyn Layer + 'static>>>> {
        self.layers.clone()
    }

    fn forward(&mut self, x: &RcVariable) -> RcVariable {
        let mut y = x.clone();

        for layer in self.layers.borrow_mut().iter_mut() {
            let t = y;
            y = layer.call(&t);
        }

        y
    }
}

//loss: Loss, optimizer: Optimizer, learning_rate: f32

impl BaseModel {
    pub fn new() -> Self {
        BaseModel {
            input: None,
            output: None,
            layers: Rc::new(RefCell::new(Vec::new())),
        }
    }

    pub fn call(&mut self, input: &RcVariable) -> RcVariable {
        // inputのvariableからdataを取り出す

        let output = self.forward(input);

        //　inputsを覚える
        self.input = Some(input.downgrade());

        //  outputを弱参照(downgrade)で覚える
        self.output = Some(output.downgrade());

        output
    }
}

Model構造体も Function構造体 のように今後様々なレイヤーを実装していきます。また複雑で独自なレイヤーを自動で初期化し、保持するといったオリジナルのモデル構造体も実装することができます。ぜひこの BaseModel を参考にして実装してみてください。

Optimizerの実装

では最後に Optimizer の実装です。Optimizerはパラメーターの更新を扱う構造体です。誤差逆伝播法のところでも触れましたが、パラメーターを更新する最適化関数の種類は様々あり、それを選択できるように様々な最適化を行う構造体を用意する必要があります。なので、こちらも今までと同様に、 Optimizerトレイト を実装し、このトレイトを実装した構造体を Optimizer構造体 と呼ぶことにします。


先ほどと同じように core.rs ファイルと同じ階層に optimizer.rs ファイルを追加します。モジュールとして認識してもらうよう、 lib.rs、mod.rs に optimizer.rs の名前を追加しておきます。

pub trait Optimizer {
    fn setup(&mut self, target_model: &impl Model);
    fn update(&self);
    fn update_param(&self, param: &RcVariable) -> ArrayD<f32>;
}
graph LR

 X((Input)) --> layer1
 subgraph Model
 layer1 --> activation1[活性化関数]
 activation1 --> layer2
 layer2 --> activation2[活性化関数]
 activation2 --> ...
 end

 ... --> pre((予測値))
pre --> Loss
t((答え)) --> Loss
Loss --> loss((L))
Optimizer --- Model
Model ---|パラメーターを共有| Optimizer

Opimizerの基本的な構造は Model構造体 からパラメーター(正確に言えば可変な参照)を受け取り、パラメーターをその構造体の最適化手法で更新します。

ではトレイトの関数について説明します。

  • setup()
    Model構造体の参照を渡すことで、Modelの保持するパラメーター(正確にはLayerを保持し、間接的にパラメーターにアクセスする)を共有することができます。

  • update()
    Modelと共有したパラメーターを更新する関数です。Function構造体 でいうと call() と同じ立ち位置です。

  • update_param()
    パラメーターのデータを実際に更新する関数です。Function構造体 でいうと forward() と同じ立ち位置です。

ではいままで説明してきた 勾配降下法 に近い最適化手法である SGD を例にして実装します。

pub struct SGD {
    lr: f32,
    layers: Option<Rc<RefCell<Vec<Box<dyn Layer + 'static>>>>>,
}

impl Optimizer for SGD {
    fn setup(&mut self, target_model: &impl Model) {
        self.layers = Some(target_model.layers().clone());
    }

    fn update(&self) {
        for layer in self
            .layers
            .as_ref()
            .expect("SGDにModelが設定されていません")
            .borrow_mut()
            .iter_mut()
            .filter(|layer| layer.has_params())
        {
            for (_id, param) in layer.params() {
                let new_param = self.update_param(&param);
                param.0.borrow_mut().data = new_param;
            }
        }
    }
    fn update_param(&self, param: &RcVariable) -> ArrayD<f32> {
        let current_param_data = param.data();
        let param_grad = param
            .grad()
            .as_ref()
            .expect("SGDで更新中のパラメータにgradがありません")
            .data();

        let new_param_data = current_param_data - self.lr * param_grad;
        new_param_data
    }
}

impl SGD {
    pub fn new(lr: f32) -> Self {
        Self {
            lr: lr,
            layers: None,
        }
    }
}

setup() で対象のModelを渡しながら、ModelのLayerを配列で保持します。update() と update_param() は Function構造体 の call() と forward() に対応していると考えると理解しやすいと思います。

今回はOptimizer構造体とModel構造体を分離させて、パラメーターを共有する設計となっています。この設計は役割がはっきり分離しているため、理解しやすい構造となっていますが、Layerをコピーしたりと、最適な設計とは言えません。例えば、Model構造体がOptimizer構造体をフィールドとして保持すれば、Model構造体の中で完結し、よりシンプルになります。

ここではSGD法を実装しましたが、他にも最適化関数は存在します。代表的なものは

  • Momentum
  • NAG
  • Adagrad
  • Adadelta
  • Rmsprop

など数多くあります。最適化の方法を調べ、SGDと同じように実装してみてください。そうすれば、フレームワークで用いることができます。

実際に使用するイメージはこうです。

let mut model = BaseModel::new();
model.stack(Dense::new(1000, true, None, Activation::Relu));
model.stack(Linear::new(10, true, None));

let mut optimizer = SGD::new(lr); //optimizer構造体の初期化
optimizer.setup(&model);          // modelの参照を保持
.
.

    model.cleargrad();      // パラメーターの勾配を初期化
    loss.backward(false);   // バックプロパゲーションを行い、パラメーターは微分を保持
    optimizer.update();     //パラメーターの微分をもとに更新

ではいよいよ今までの構造体を組み合わせ、自動化されたニューラルネットワークを構築しましょう。

Configの実装

ここではフレームワーク全体の設定を管理するConfigファイルを実装します。

今までの同様にに core.rs ファイルと同じ階層に config.rs ファイルを追加し、lib.rs、mod.rs に名前を追加しておきます。

use std::sync::atomic::{AtomicBool, AtomicUsize, Ordering};
/// Variableや関数たちにidを付けるための値
pub static NEXT_ID: AtomicUsize = AtomicUsize::new(0);

/// 微分するかしないかというフラグ
/// 推論するときなど、微分する必要がないときに切り替える
pub static GRAD_CONFIG: AtomicBool = AtomicBool::new(true);

/// idを生成する関数。構造体のコンストラクタを作成する際に、呼び出して、idを付ける
pub fn id_generator() -> usize {
    NEXT_ID.fetch_add(1, Ordering::SeqCst)
}

pub fn set_grad_true() {
    GRAD_CONFIG.store(true, Ordering::SeqCst);
}

pub fn set_grad_false() {
    GRAD_CONFIG.store(false, Ordering::SeqCst);
}

pub fn get_grad_status() -> bool {
    GRAD_CONFIG.load(Ordering::SeqCst)
}

はじめに id に関する設定を行います。微分の実装（複雑な関数）で実装した id を扱う NEXT_ID、 id_generator() をここに移します。また高階微分を行うかというフラグ GRAD_CONFIG の設定もここで行います。このフラグの状態、または変更する関数を実装します。

高階微分に関しては補足のところで解説します。

TODO:高階微分補足の場所で説明

ニューラルネットワークの学習

では早速今まで実装してきた構造体を組み合わせて、はじめの手動でのニューラルネットワークを自動化して学習させます。

use ndarray::{array, Array};
use std::array;
use std::cell::RefCell;
use std::rc::Rc;
use std::time::Instant;
use stucrs::core_new::{F32ToRcVariable, RcVariable};
use stucrs::functions_new as F;
use stucrs::layers::{self as L, Dense, Layer, Linear};
use stucrs::models::{BaseModel, Model};
use stucrs::optimizers::{Optimizer, SGD};

use std::f32::consts::PI;

use stucrs::core_new::ArrayDToRcVariable;

use ndarray_rand::rand_distr::{StandardNormal, Uniform};
use ndarray_rand::RandomExt;

fn main() {
    let x_data = Array::random((100, 1), Uniform::new(0.0f32, 1.0));
    //let x2 = array![[11.0f32, 12.0, 13.0], [14.0, 15.0, 16.0]].rv();

    let y_data =
        x_data.mapv(|x| (2.0 * PI * x).sin()) + Array::random((100, 1), Uniform::new(0.0f32, 1.0));

    let x = x_data.rv();

    let y = y_data.rv();

    let l1 = L::Dense::new(10, false, None, L::Activation::Sigmoid);
    let l2 = L::Linear::new(10, false, None);

    let mut model = BaseModel::new();
    model.stack(l1);
    model.stack(l2);

    let lr = 0.1;

    let mut optimizer = SGD::new(lr);

    optimizer.setup(&model);


    let iters = 10000;
    for i in 0..iters {
        let y_pred = model.call(&x);

        let mut loss = F::mean_squared_error(&y, &y_pred);

        model.cleargrad();
        loss.backward(false);

        optimizer.update();

        if i % 1000 == 0 {
            println!("i = {:?}", i);
            println!("loss = {:?}\n", loss.data());
        }
    }
}

前回のコードを見てもらうとわかるように、パラメーターを処理するコードが劇的に少なくなりました。また predict() という関数をModelとして設定し、Optimizerを設定することで、パラメーター更新も update() のみですべてのパラメーターを更新することができました。では改めて、同じデータで学習させてみましょう。


ここまでで、 ニューラルネットワークの骨格となるシステムが構築しました。 あとはよりニューラルネットワークのフレームワークとして機能を仕上げていきます。次に実装する機能は 損失関数 です。

損失関数

私たちは実装した損失関数は 二乗誤差関数 は主に回帰に適した関数であり、このように損失関数は様々な種類があり、それぞれある学習に特化した関数です。続いて私たちは分類に特化した損失関数である クロスエントロピー誤差 、並びにそれとよくセットで用いられる softmax関数 も実装していきます。

クロスエントロピー誤差とsoftmax関数は多値分類を行うモデルで一般的にセットで使用されます。なので、既存のフレームワークの中ではクロスエントロピー誤差とsoftmaxを一つにまとめているものもあります。

クロスエントロピー誤差とはどのような損失関数なのでしょうか。

graph LR

 X((Input)) --> layer1
 subgraph Model
 layer1 --> ...
 end

 ... --> pre((予測値))
pre --> Softmax
Softmax --> Cross_entropy
t((答え)) --> Cross_entropy
Cross_entropy --> loss((L))

先ほどの説明のようにはじめにModelが求めた予測値をSoftmax関数に通した値を Cross_entropy に流すことで分類を行うことができるという仕組みです。ここで softmax　は予測値を確率に変換する役割があります。次に Cross_entropy を処理を表します。

$${L = -\sum_{j} t_j\cdot\log p_j}$$


\(t_j\)は正解ラベルを表しますが、今回実装する Cross_entropy は正解ラベルを one_hotベクトル で渡す設計にしますです。これにより、処理が簡潔になります。 one_hotベクトル とはこのような行列データです。例えばあるデータを0～2の3つに分類するとします。この時、クラス数は3なので one_hotベクトル の行列の列数は3です。左の正解ラベルの値のインデックスを1にし、他の値が全て0にすることで、かけることにより、正解ラベルの時の情報だけを維持し、他のデータは捨てることができます。またモデルの出力はクラス数の長さの行列をバッチ数の数分の行列を返すので、この one_hotベクトル の行列と形状が一致します。


$$ T:\begin{pmatrix} 0 \\ 1 \\ 2\\ 1\\ 2\\ 0\\ \vdots\\ \end{pmatrix} \xrightarrow{\text{one-hot-vector}} T’:\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \\ \vdots & \vdots & \vdots \\ \end{pmatrix} $$


$$ Y:\begin{pmatrix} y_0 & y_1 & y_2\\ y_3 & y_4 & y_5 \end{pmatrix} \xrightarrow{\text{Softmax}} P:\begin{pmatrix} p_0 & p_1 & p_2\\ p_3 & p_5 & p_6 \end{pmatrix} \xrightarrow{\text{Log}} P_{log}:\begin{pmatrix} \log p_0 & \log p_1 & \log p_2\\ \log p_3 & \log p_4 & \log p_5 \end{pmatrix} $$


$$ -P_{log}:\begin{pmatrix} -\log p_0 & -\log p_1 & -\log p_2\\ -\log p_3 & -\log p_4 & -\log p_5 \end{pmatrix} \cdot T’:\begin{pmatrix} t_0 & t_1 & t_2\\ t_3 & t_4 & t_5 \end{pmatrix} \xrightarrow{\text{}} -P_{log}:\begin{pmatrix} -t_0\cdot\log p_0 & -t_1\cdot\log p_1 & -t_2\cdot\log p_2\\ -t_3\cdot\log p_3 & -t_4\cdot\log p_4 & -t_5\cdot\log p_5 \end{pmatrix} $$


先ほどの数式、

$${L = -\sum_{j} t_j\cdot\log p_j}$$

を改めて考えてみます。\(t_j\)はone_hotベクトル の行列の要素なので、値は0か1です。つまり、掛け算をすることで、この\(t_j\)が1の時、すなわち正しく予測できた場合のみ加算されます。誤差を表す\(L\) は正解したときのデータのみを足し合わせた値をマイナスにするので、正解が多いほどこの値自体は小さくなります。これによって誤差を表しているのです。

でははじめにクロスエントロピー誤差を実装するにあたって先に softmax関数、Log関数を実装していきます。

Softmax関数

ではsoftmax関数の数式、そしてどのような処理なのか見ていきます。

$${p_i = \frac{exp(y_i)}{\sum_{j=1}^N exp(y_j)}}$$

ではデータ数が2の場合を見てみましょう。


$$ X:\begin{pmatrix} x_0 & x_1 & x_2 \\ x_3 & x_4 & x_6 \end{pmatrix} \xrightarrow{\text{予測}} Y:\begin{pmatrix} y_0 & y_1 & y_2\\ y_3 & y_4 & y_5 \end{pmatrix} $$


$$ Y:\begin{pmatrix} y_0 & y_1 & y_2\\ y_3 & y_4 & y_5 \end{pmatrix} \xrightarrow{\text{Softmax}} P:\begin{pmatrix} \frac{e^{y_0}}{e^{y_0}+e^{y_1}+e^{y_2}} & \frac{e^{y_1}}{e^{y_0}+e^{y_1}+e^{y_2}} & \frac{e^{y_2}}{e^{y_0}+e^{y_1}+e^{y_2}}\\ \frac{e^{y_3}}{e^{y_3}+e^{y_4}+e^{y_5}} & \frac{e^{y_4}}{e^{y_3}+e^{y_4}+e^{y_5}} & \frac{e^{y_5}}{e^{y_3}+e^{y_4}+e^{y_5}} \end{pmatrix} $$


文字が小さくなってしまいましたが、要は同じ行の中でネイピア数の指数にかけて平均を取っているようなものです。ではなぜ指数を取るのかというとそれは負の値であっても処理できるからです。指数は正負かかわらず必ず正になるので、どんな値でも処理することができます。

ではsoftmax関数を実装していきます。この処理は二乗平均誤差と同じように、分母となる値を sum で求め、それを行列として割ればよいのです。

$$ Y:\begin{pmatrix} y_0 & y_1 & y_2\\ y_3 & y_4 & y_5 \end{pmatrix} \xrightarrow{\text{Exp}} Y_{exp}:\begin{pmatrix} e^{y_0} & e^{y_1} & e^{y_2}\\ e^{y_3} & e^{y_4} & e^{y_5} \end{pmatrix} \xrightarrow{\text{Sum(axis=1)}} Y_{sum-exp}:\begin{pmatrix} e^{y_0} + e^{y_1} + e^{y_2} \\ e^{y_3} + e^{y_4} + e^{y_5} \end{pmatrix} $$


$$ Y_{exp}:\begin{pmatrix} e^{y_0} & e^{y_1} & e^{y_2}\\ e^{y_3} & e^{y_4} & e^{y_5} \end{pmatrix} \div Y_{sum-exp}:\begin{pmatrix} e^{y_0} + e^{y_1} + e^{y_2} \\ e^{y_3} + e^{y_4} + e^{y_5} \end{pmatrix} \xrightarrow{\text{Softmax}} P:\begin{pmatrix} \frac{e^{y_0}}{e^{y_0}+e^{y_1}+e^{y_2}} & \frac{e^{y_1}}{e^{y_0}+e^{y_1}+e^{y_2}} & \frac{e^{y_2}}{e^{y_0}+e^{y_1}+e^{y_2}}\\ \frac{e^{y_3}}{e^{y_3}+e^{y_4}+e^{y_5}} & \frac{e^{y_4}}{e^{y_3}+e^{y_4}+e^{y_5}} & \frac{e^{y_5}}{e^{y_3}+e^{y_4}+e^{y_5}} \end{pmatrix} $$

Exp で行列の要素全体に指数を施し、Sumで軸を1に指定して合計することで、分母となる値を持つ行列を計算します。あとはその行列で割れば、求める処理を実現できます。ではこれをコードで実装します。

pub fn softmax_simple(x: &RcVariable) -> RcVariable {
    let exp_y = exp(&x);

    let sum_y = sum(&exp_y, Some(1)); //sumで軸1指定

    let y = exp_y.clone() / sum_y.clone();
    y
}

では実装した softmax関数 をテストします。
TODO: softmax関数テストコード

このようにsoftmax関数は出力されたデータを確率に変換します。この変換はニューラルネットワークの分類において重要な処理です。

では続いてクロスエントロピー誤差を実装するにあたり必要な Log関数 を Function構造体 として実装します。

Log関数の実装

ではクロスエントロピー誤差に必要な関数である Log関数 を実装していきます。この関数は以前の Function構造体 なので同じように実装します。ただし、オプションが多いので、微分の式など処理が多少複雑なので、説明します。


まず数学の分野では明示しなければLog関数の底は\(e\)というのが暗黙の了解なので、sum関数 の時ど同様、引数がNone の時を底\(e\)とします。そして底が他の値の場合はSome(base) として渡すようにします。もちろんですが、base の型はもちろん普通の数値と同じ f32 です。

また、この底の値で微分のふるまいが異なります。正確に言えば、統一した方法で計算するは可能ですが、分けた方が、数学的な理解としても、パフォーマンス的にしても良いので(\(e\)以外の底の値を使用することはまれなので)、分けます。

では実装する前にlogの微分を考えてます。今回は底で場合分けして考えます。まず\(y = \log x\)、もしくは\(y = \ln x\)の時の微分は

$$\frac{dy}{dx} = \frac{1}{x}$$

となります。次に底が指定された場合、つまり\(y = \log_a x\)の時の微分は

$$\frac{dy}{dx} = \frac{1}{x\cdot \ln a}$$

となります。ではこれらをもとに Log関数 を実装していきます。

struct Log {
    inputs: Vec<RcVariable>,
    output: Option<Weak<RefCell<Variable>>>,
    base: Option<f32>,
    generation: i32,
    id: usize,
}

impl Function for Log {
    fn call(&mut self) -> RcVariable {
        let inputs = &self.inputs;
        if inputs.len() != 1 {
            panic!("Logは一変数関数です。inputsの個数が一つではありません。")
        }

        let output = self.forward(inputs);

        if get_grad_status() == true {
            //inputのgenerationで一番大きい値をFuncitonのgenerationとする
            self.generation = inputs.iter().map(|input| input.generation()).max().unwrap();

            //  outputを弱参照(downgrade)で覚える
            self.output = Some(output.downgrade());

            let self_f: Rc<RefCell<dyn Function>> = Rc::new(RefCell::new(self.clone()));

            //outputsに自分をcreatorとして覚えさせる
            output.0.borrow_mut().set_creator(self_f.clone());
        }

        output
    }

    fn forward(&self, xs: &[RcVariable]) -> RcVariable {
        let base = self.base;
        let x = &xs[0];
        let y_data;

        //baseがeか他の値かで場合分け(eの場合、baseはNone)
        if let Some(base_data) = base {
            y_data = x.data().mapv(|x| x.log(base_data));
        } else {
            y_data = x.data().mapv(|x| x.ln());
        }
        y_data.rv()
    }

    fn backward(&self, gy: &RcVariable) -> Vec<RcVariable> {
        let x = &self.inputs[0];
        let gx;

        let base = self.base;

        //baseがeか他の値かで場合分け(eの場合、baseはNone)
        if let Some(base_data) = base {
            gx = 1.0.rv() / (x.clone() * base_data.ln().rv()) * gy.clone();
        } else {
            gx = (1.0.rv() / x.clone()) * gy.clone();
        }
        let gxs = vec![gx];
        gxs
    }

    fn get_inputs(&self) -> &[RcVariable] {
        &self.inputs
    }

    fn get_output(&self) -> RcVariable {
        let output;
        output = self
            .output
            .as_ref()
            .unwrap()
            .upgrade()
            .as_ref()
            .unwrap()
            .clone();

        RcVariable(output)
    }

    fn get_generation(&self) -> i32 {
        self.generation
    }
    fn get_id(&self) -> usize {
        self.id
    }
}
impl Log {
    fn new(inputs: &[RcVariable], base: Option<f32>) -> Rc<RefCell<Self>> {
        Rc::new(RefCell::new(Self {
            inputs: inputs.to_vec(),
            output: None,
            base: base,
            generation: 0,
            id: id_generator(),
        }))
    }
}

pub fn log(x: &RcVariable, base: Option<f32>) -> RcVariable {
    let y = log_f(&[x.clone()], base);
    y
}

fn log_f(xs: &[RcVariable], base: Option<f32>) -> RcVariable {
    Log::new(xs, base).borrow_mut().call()
}

底の指定を Sum関数 と同じようにOption型 で渡し、forward,backward で場合分けして処理します。log関数の計算に慣れていればそれほど難しくはないでしょう。

では底で場合分けしてテストします。

#[test]
    fn log_test() {
        use crate::core_new::ArrayDToRcVariable;

        let a = array![3.0, 3.0, 3.0].rv();
        let b = array![3.0, 3.0, 3.0].rv();

        let mut y0 = log(&a, None); //底がe
        let mut y1 = log(&b, Some(2.0)); //底が2.0

        println!("y0 = {}", y0.data()); // 1.098...
        println!("y1 = {}", y1.data()); // 1.584...

        y0.backward(false);
        y1.backward(false);

        println!("a_grad = {:?}", a.grad().unwrap().data()); // 0.3333...
        println!("b_grad = {:?}", b.grad().unwrap().data()); // 0.4808...
    }

クロスエントロピー誤差

では最後にクロスエントロピー誤差を実装していきます。今回は、最初の説明のように softmax 関数もセットにした関数にします。

pub fn softmax_cross_entropy_simple(x: &RcVariable, t: &RcVariable) -> RcVariable {
    if x.data().shape() != t.data().shape() {
        panic!("交差エントロピー誤差でのxとtの形状が異なります。tがone-hotベクトルでない可能性があります。")
    }

    let n = x.data().shape()[0] as f32;

    let p = softmax_simple(&x);

    let clamped_p = clamp(&p, 1.0e-15, 1.0); // pを1.0e-15～1.0に収める

    let log_p = log(&clamped_p, None); // 0の値は渡されない

    let tlog_p = log_p * t.clone();

    let y = (-sum(&tlog_p, None)) / n.rv();
    y
}

引数の\(x\)はモデルの出力、\(t\)は正解ラベルを想定しています。 先ほどの説明のように、正解ラベルのデータは one_hotベクトル に変換することでモデルの出力と同じ形状にし、掛け算をします。なので形状一致を確認します。

次にここで初めて登場した clamp を説明します。この関数は引数に p 1.0e-15, 1.0 の二つの値を渡していますが、最初の変数の値をその二つの値内に収める処理を行います。つまり、p が\(1.0^{-15}\leqq p\leqq 1.0\)に収まるならそのまま値を流し、この範囲からはみ出たら、例えば2.0なら1.0にするというように修正して返します。これは次に処理する log関数 でlogに0の値を渡すことを防ぐためのものです。この関数に関しては 補足 で説明を行います。

TODO:clamp関数 補足で説明

では実装したクロスエントロピー誤差を試してみましょう。データやModelを自由にを用意、構築し、誤差を求めてみましょう。

fn main() {
    let x = array![[0.1f32, -0.4], [0.3, 0.6], [-1.3, -1.2], [3.1, -0.5]].rv();
    let t = array![
        [0.0f32, 0.0, 1.0],
        [1.0, 0.0, 0.0],
        [0.0, 1.0, 0.0],
        [1.0, 0.0, 0.0]
    ]
    .rv(); // one_hotベクトルで渡す

    let mut model = BaseModel::new();
    model.stack(L::Dense::new(5, false, None, Activation::Sigmoid));
    model.stack(L::Dense::new(3, false, None, Activation::Sigmoid));

    let y = model.call(&x);

    let loss = F::softmax_cross_entropy_simple(&y, &t);

    println!("loss = {:?}", loss.clone().data()); //誤差が求まる
}

これらにより、Modelは多クラスにおける予測値と答えとの誤差を求めることができました。次はこの多クラス分類を用いてより実践的な学習を試みます。

データセットの用意

ここではデータセットをモデルが扱えるように処理、管理する構造体を実装していきます。
今までのデータセットと言えば、ニューラルネットワークの構築の際に用いた\(y = sin(2\pi\cdot x) + b\)があります。これは私たちがテストするために用意した簡易的なデータです。しかし、深層学習で学習させるデータセットは莫大なデータ数で、とても複雑なものです。これを毎回人の手でニューラルネットワークにうまく変換することは大変です。なので、それらを一元的に処理する構造体を実装していきます。

こちらの構造体も今までと同様に、Datasetトレイト を実装し、様々な構造体に継承させるので、Datasetトレイトを実装した構造体をDataset構造体 と呼ぶことにします。

先ほどと同じように core.rs ファイルと同じ階層に dataset.rs ファイルを追加します。モジュールとして認識してもらうよう、 lib.rs、mod.rs に dataset.rs の名前を追加しておきます。

pub trait Dataset {
    fn len(&self) -> usize;
    fn data_setup(&mut self);
}

こちらが基本的なトレイトとなります。はじめは比較的簡易的なデータセットを扱うので、初歩的な機能のみを実装させます。fn len() はデータ数を返し、fn data_setup は実際にデータを初期化する処理を行います。

スパイラルデータの用意

では次に学習させるスパイラルデータを例に挙げて実際に実装していきます。

スパイラルデータは「ゼロから作るDeep Learning3　フレームワーク編」で使用されているスパイラルデータを参考にしています。

pub struct Spiral {
    pub data: Array2<f32>,
    pub label: Array1<u32>,
}

impl Dataset for Spiral {
    fn data_setup(&mut self) {}

    /*
    fn get_item(&self, index: i32) -> ArrayViewD<f32> {
        self.data.view().into_dyn()
    }
    */
    fn len(&self) -> usize {
        self.data.shape()[0]
    }
}

impl Spiral {
    pub fn new() -> Self {
        let data_label = get_spiral_data();
        let data = data_label.0;
        let label = data_label.1;
        let spiral = Self {
            data: data,
            label: label,
        };
        spiral
    }
}

fn get_spiral_data() -> (Array2<f32>, Array1<u32>) {
    let data_len = 100;
    let num_class = 3;
    let input_dim = 2;

    let data_size = data_len * num_class;

    let mut x = Array2::zeros((data_size, input_dim));

    let mut t = Array1::zeros(data_size);
    //let normal_dist = Normal::new(0.0f32, 1.0).unwrap();

    for j in 0..num_class {
        for i in 0..data_len {
            let rate = i as f32 / data_len as f32;
            let radius = 1.0 * rate as f32;

            let mut rng = rand::thread_rng();
            let normal: f32 = rng.sample(StandardNormal);

            let theta = j as f32 * 4.0 + 4.0 * rate as f32 + normal * 0.2;

            let ix = data_len * j + i;
            let mut x_row_view = x.row_mut(ix);

            let row_array = array![radius as f32 * theta.sin(), radius as f32 * theta.cos()];

            x_row_view.assign(&row_array);
            t[ix] = j as u32;
        }
    }
    //一つ2次元の行列と一次元の行列の対となる行が同じ位置に来るようにシャッフルして新しい二つの行列を返す
    double_matrix_shuffle_rows_immutable(x.view(), t.view())
}

pub fn double_matrix_shuffle_rows_immutable(
    arr1: ArrayView2<f32>,
    arr2: ArrayView1<u32>,
) -> (Array2<f32>, Array1<u32>) {
    if arr1.nrows() != arr2.len() {
        panic!("arr1とarr2の行列の行数が異なります")
    }
    // 行のインデックスを作成 (0, 1, 2, ...arr1.nrows())
    let mut indices: Vec<usize> = (0..arr1.nrows()).collect();

    // インデックスをシャッフル
    indices.shuffle(&mut thread_rng());
    
    let new_arr1 = arr1.select(Axis(0), &indices);
    let new_arr2 = arr2.select(Axis(0), &indices);

    (new_arr1.to_owned(), new_arr2.to_owned())
}

スパイラルデータを扱う構造体を spiral 構造体とします。今回のデータの生成は比較的簡単なので、data_setup()を使わず、get_spiral_data() で直接生成してフィールドに保持させます。このデータをrustのグラフを扱うライブラリである plotters で表示させると、渦巻きのような模様が現れます。そして、この三色の色に分類するというわけです。このグラフ表示はplottersを使うで解説します。

そして、生成されたデータを正解ラベルとともにシャッフルしていきます。というのも、学習させるデータをシャッフルしてモデルに渡すことで、順番による影響を最小限にすることができます。ここで、データとラベルがずれることがないようにシャッフルする必要があります。この処理は double_matrix_shuffle_rows_immutable() で行います。今回は二次元の座標データなので、次元は2次元\((x,y)\)です。


$$ data:\begin{pmatrix} x_0 & y_0 \\ x_1 & y_1 \\ x_2 & y_2 \\ x_3 & y_3 \\ x_4 & y_4 \\ x_5 & y_5 \\ \vdots & \vdots \\ x_N & y_N \end{pmatrix} \qquad
T:\begin{pmatrix} t_0 \\ t_1 \\ t_2\\ t_3\\ t_4\\ t_5\\ \vdots\\ t_N \end{pmatrix} \quad \xrightarrow{\text{shuffle}} \quad data’:\begin{pmatrix} x_4 & y_4 \\ x_N & y_N \\ x_2 & y_2 \\ x_3 & y_3 \\ x_1 & y_1 \\ x_0 & y_0 \\ \vdots & \vdots \\ x_5 & y_5 \end{pmatrix} \qquad
T’:\begin{pmatrix} t_4 \\ t_N \\ t_2\\ t_3\\ t_1\\ t_0\\ \vdots\\ t_5 \end{pmatrix} $$


これはあくまでシャッフルの一例ですが、大事なのは、データとラベルの順番がしっかり対応しているということです。これにより、シャッフルされた状態でも、同じデータとラベルの組を正しく取り出すことができます。

この学習データは次の学習で用いるので、実装しておきましょう。

accuracy関数

ここではモデルがどれくらい精度が高いかを評価する accuracy関数 を実装していきます。これはただモデルの性能を測るだけの関数なのでモデルの学習には使用しません。つまり、バックプロパゲーション対象外の微分をしない関数ということです。

モデルの精度を測る関数はいくつか存在しますが、この accuracy関数 は一番オーソドックスな関数で、処理としてはただ単に正解数をデータ数で割り正解率を求めるというシンプルな計算ですが、この正解数を求める処理を考える必要があります。


そもそもモデルが正しく答えを出せたとはどういうことでしょうか。まずはそこから考えます。モデルの学習とはいかに誤差を小さくするか、そして、もっと言えば正解ラベルに近い出力データを出せるかということです。そのような中で、私たちはモデルの予測値が答えにどれほど近ければ正解したとするかを定義しなければなりません。先ほど softmax関数 の中で確率に変換すると説明しましたが、つまりモデルは予測値をあくまで確率として出力するのです。例えば、クラス数が3なら、0の確率は0.1、1は0.3、2は0.6という感じです。なので、私たちは、モデルが正解の値を一番高い確率で出力できたら正解したとします。この場合、モデルは確率の値として一番高い2だと予測したとらえ、それが正解ラベルと一致しているかを確かめればよいのです。


モデルの正解の定義を考えましたが、次はそれを実装していきます。重要なことは正解ラベルの値はインデックスを表しているということです。

pub fn accuracy(y: ArrayView2<f32>, t: ArrayView2<f32>) -> f32 {
    if y.shape() != t.shape() {
        panic!("交差エントロピー誤差でのxとtの形状が異なります。tがone-hotベクトルでない可能性があります。")
    }
    let data_size = y.shape()[0] as f32;
    let num_class = t.shape()[1];
    let argmax_vec: Vec<u32> = y
        .outer_iter()
        .map(|row: ArrayBase<ViewRepr<&f32>, Dim<[usize; 1]>>| row.argmax().unwrap() as u32)
        .collect();
    let max_index = Array::from_vec(argmax_vec);
    let one_hot_y = arr1d_to_one_hot(max_index.view(), num_class);

    assert_eq!(one_hot_y.shape(), t.shape());

    let acc_matrix = &one_hot_y * &t;

    let accuracy = acc_matrix.sum() / data_size; // 正解数 / データ数

    accuracy
}

//二つの行列をかけることで正解数がわかる。
    //　例
    // one_hot_y = [[0.0,0.0,1.0],
    //            [0.0,1.0,0.0],
    //            [1.0,0.0,0.0],]
    //
    //    t  =   [[0.0,0.0,1.0],
    //            [1.0,0.0,0.0],
    //            [1.0,0.0,0.0],]
    //
    //one_hot_y* t = [[0.0,0.0,1.0],　　　　2行目は1.0の位置が違うのでかけてすべてゼロ。
    //            [0.0,0.0,0.0],
    //            [1.0,0.0,0.0],]      行列を要素ごとにかけることで、正しい答えだった場合のみ1.0になり、その他はすべて0.0になる。
    //                                 この行列の合計値sum関数で正解数がわかる。この場合、sumの値は2なので正解数は2。

コメントでも説明していますが、予測値で一番値が高いラベルのところを1にし、他を0にするというone_hotベクトル化をします。これにより、正解ラベルの行列と掛けることで正解したところだけ1が残ります。あとは行列でSum を取れば、正解した個数がわかります。あとはそれをデータ数(正確にはバッチ数)で割れば、そのデータ数(バッチ数)での正解率が求まります。


ここで一つ気になる点として行列の次元を静的で扱っていることが挙げられます。この理由は主に二つあります。一つ目は、最大値を取るインデックスを求める、numpyでいうargmax() のような処理(軸を指定した最大値のインデックスを求める関数はrustのndarrayには実装されていないので、自分で処理を書きます)を行う際、静的な次元で行う必要があるからです。二つ目は、基本的に正解ラベルの行列の形状は2次元だからです。というのも、正解のデータは0,1,2,・・というように行列ではなく、数値で表せるからです。しかし、今後、高次元の正解ラベル(行列の正解データはもはやラベルとは呼べないかもしれませんが)を扱う新たなアルゴリズムでが登場した場合は対応できませんので、より柔軟な関数が求めらます。

また、この関数はバックプロパゲーション対象外なので、Function構造体 や RcVariable といったものを考える必要はなく、ただ単にndarrayの行列計算に集中してよいというわけです。

コードで登場した arr1d_to_one_hot() は一次元の行列から2次元のone_hotベクトルの行列を生成します。損失関数の one_hotベクトル の\(T\)から\(T’\)の変換説明を見るとわかりやすいと思います。この arr1d_to_one_hot() については 補足 で説明します。

TODO:arr1d_to_one_hot()補足

では実際に正解率を求められるかテストします。

fn main() {
    let a: Array2<f32> = array![
        [1.0f32, 2.0, 3.0],  // ○
        [6.0, 4.0, 5.0],     // ○
        [1.0, 2.0, 10.0],    // ×
        [12.0, 15.0, 5.0]    // ○
    ];
    let b = array![2, 0, 1, 1];
    let b = to_one_hot(b.view(), 3);

    let acc = accuracy(a.view(), b.view());
    println!("{}", acc);
}

スパイラルデータの学習

ではいよいよ実際に多クラス分類の学習を行います。先ほどの説明でスパイラルデータの準備は整いました。三色に分類するのでクラス数は3です。この学習は今までのものよりも少し規模が大きくなるので、いくつか準備が必要です。

ミニバッチ処理

ここではじめてバッチというものが登場しました。バッチとはデータを分割して、一度にモデルに渡すデータの数を調整することです。ここで一度に渡すデータ数をバッチ数と言います。今までのデータはデータ数が少なかったのですが、実際に使われる場合は数万個に及ぶデータを扱います。全データを一度にモデルに渡すと、メモリ使用量も莫大になり、パラメーターも収束しづらいです。なので、大きいデータ数の場合は、分割して渡すのが主流です。私たちが計算を行列で行ってきたのも、このミニバッチ学習に対応するためなのです。全データからいくつかのデータとラベルを取り出し、モデルに渡します。この行列データからランダムに複数の行を取り出す必要がありますが、その処理をndarrayに搭載される chunks() メソッドで行えます。

バッチサイズが1の時、つまりデータを一つずつモデルに渡す学習をオンライン学習、バッチサイズが全データ数の時、つまり全データをまとめてモデルに渡す学習をバッチ学習、その中間をミニバッチ学習と言います。

use ndarray::*;
use ndarray_stats::QuantileExt;
use std::array;
use std::cell::RefCell;
use std::f32::consts::PI;
use std::rc::Rc;
use std::time::Instant;
use stucrs::config;
use stucrs::core_new::ArrayDToRcVariable;
use stucrs::core_new::{F32ToRcVariable, RcVariable};
use stucrs::datasets::*;
use stucrs::functions_new::{self as F, accuracy, sum};
use stucrs::layers::{self as L, Activation, Dense, Layer, Linear};
use stucrs::models::{BaseModel, Model};
use stucrs::optimizers::{Optimizer, SGD};

use plotters::prelude::*;
use rand::seq::SliceRandom;
use rand::*;

fn main() {
    let max_epoch = 300;
    let lr = 1.0;
    let batch_size = 30;
    
    let train_spiral = Spiral::new(true);
    let test_spiral = Spiral::new(true);

    let x_train = train_spiral.data.view();
    let y_train = train_spiral.label.view();
    let y_train = to_one_hot(y_train, 3);

    let x_test = test_spiral.data.view();
    let y_test = test_spiral.label.view();
    let y_test = to_one_hot(y_test, 3);

    let mut model = BaseModel::new();
    model.stack(L::Dense::new(10, true, None, Activation::Sigmoid));
    model.stack(L::Linear::new(3, true, None));

    let data_size = x_train.shape()[0];
    let mut optimizer = SGD::new(lr);
    optimizer.setup(&model);

    for epoch in 0..max_epoch {
        let mut indices: Vec<usize> = (0..data_size).collect();
        let mut rng = thread_rng();
        indices.shuffle(&mut rng);
        let mut sum_loss = array![0.0f32];
        let mut sum_acc = array![0.0f32];

        // バッチを取り出す
        for chunk_indices in indices.chunks(batch_size) {
            let x_batch = x_train.select(Axis(0), chunk_indices).to_owned().rv();
            let y_batch = y_train.select(Axis(0), chunk_indices).to_owned().rv();

            //println!("x_batch = {:?}, t_batch = {:?}", x_batch, t_batch);

            let y = model.call(&x_batch);
            let mut loss = F::softmax_cross_entropy_simple(&y, &y_batch);
            let acc = accuracy(
                y.data().into_dimensionality().unwrap().view(),
                y_batch.data().into_dimensionality().unwrap().view(),
            );
            model.cleargrad();
            loss.backward(false);
            optimizer.update();

            //ここでt_batch.lenはu32からf32に変換、さらに暗黙的にndarray型に変換されて、計算される。
            //また、sum_lossは静的次元なので、epoch_lossを動的次元から静的次元に変換して足せるようにする。

            let epoch_loss: Array1<f32> = (&loss.data() * (y_batch.len() as f32))
                .into_dimensionality()
                .unwrap();

            sum_loss = &sum_loss + &epoch_loss; //全てのバッチで合計を取る
            sum_acc = sum_acc + acc * (y_batch.len() as f32); //全てのバッチで合計を取る
        }

        let average_loss = &sum_loss / (data_size as f32); //平均をとる
        let average_acc = sum_acc / (data_size as f32);    //平均をとる

        println!(
            "epoch = {:?}, train_loss = {:?}, accuracy = {}",
            epoch + 1,
            average_loss,
            average_acc
        );

        //推論
        config::set_grad_false();

        let mut indices: Vec<usize> = (0..data_size).collect();
        let mut rng = thread_rng();
        indices.shuffle(&mut rng);
        let mut sum_loss = array![0.0f32];
        let mut sum_acc = array![0.0f32];

        for chunk_indices in indices.chunks(batch_size) {
            let x_batch = x_test.select(Axis(0), chunk_indices).to_owned().rv();
            let y_batch = y_test.select(Axis(0), chunk_indices).to_owned().rv();

            //println!("x_batch = {:?}, t_batch = {:?}", x_batch, t_batch);

            let y = model.call(&x_batch);
            let mut loss = F::softmax_cross_entropy_simple(&y, &y_batch);
            let acc = accuracy(
                y.data().into_dimensionality().unwrap().view(),
                y_batch.data().into_dimensionality().unwrap().view(),
            );

            let epoch_loss: Array1<f32> = (&loss.data() * (y_batch.len() as f32))
                .into_dimensionality()
                .unwrap();

            sum_loss = &sum_loss + &epoch_loss;
            sum_acc = sum_acc + acc * (y_batch.len() as f32);
        }

        let average_loss = &sum_loss / (data_size as f32);
        let average_acc = sum_acc / (data_size as f32);

        println!(
            "epoch = {:?}, test_loss = {:?}, test_accuracy = {}",
            epoch + 1,
            average_loss,
            average_acc
        );

        config::set_grad_true();
    }
}

for文でバッチごとに誤差と正解率を算出し、すべてのデータを学習し終わったらそれらの平均を取ります。この値が、1エポックにおける学習の評価となります。今回は最初の設定から、エポック数300、学習率1.0、バッチ数30としています。前半でバックプロパゲーションを行い、パラメーターを更新します。ここが、学習のところです。そして、後半のコードは推論を行います。つまり、テストです。

エポックの学習を繰り返していくうちに、誤差が減少し、正解率が上がるのがわかると思います。

以上により、スパイラルデータで多クラス分類の学習を行うことができました。

では最後に、フレームワークの実験でよく用いられるデータセット、MNIST を学習させて、フレームワークの基礎編を終わりたいと思います。

MNISTの学習

ではいよいよ最後のMNISTの学習です。MNISTはよくニューラルネットワークの学習で試験的に用いられるデータセットで、これを学習できれば基本的な深層学習のフレームワークとして確立することができます。

データの用意

MNISTのデータセットはwebからダウンロードし、解凍して利用する形です。その際、ライブラリクレートである mnistライブラリ を使用します。このデータの用意に関しては補足のMNISTのデータの用意で説明します。内容としては、先ほどの説明した Dataset構造体 として Mnist構造体 を実装します。 使用例としてはこのようになります。

// ライブラリの呼び出しは省略しています

fn main() {
    let mnist = MNIST::new();
    let x_train = mnist.train_img;
    let y_train = mnist.train_label;
    let x_test =mnist.test_img;
    let y_test = mnist.test_label;

    let image_num = 0;

    println!("{:#.1?}\n",x_train.slice(s![image_num, .., ..]));
    println!("{:?}", x_train.shape());
}

xは画像データ、yは正解ラベルを表していて、trainは学習データ、testはテストデータで分割しています。最後のprintln!("{:#.1?}\n",x_train.slice(s![image_num, .., ..])); で0番目の画像の行列を表示します。MNISTの画像データは28×28のデータです。枚数は学習データが6万枚、テストデータが1万枚なので、学習データの行列の形状は(60000,28,28)となります。MNIST構造体の保持するデータは全てndarray型 なので、そのままモデルに渡すことができます。(正確に言うと、MNIST構造体の中の処理でndarray型に変換します。)

use image::{GrayImage, Luma};
use ndarray::*;
use ndarray_stats::QuantileExt;
use plotters::prelude::*;
use rand::seq::SliceRandom;
use rand::*;
use std::array;
use std::cell::RefCell;
use std::f32::consts::PI;
use std::rc::Rc;
use std::time::Instant;
use stucrs::config;
use stucrs::core_new::ArrayDToRcVariable;
use stucrs::core_new::{F32ToRcVariable, RcVariable};
use stucrs::dataloaders::DataLoader;
use stucrs::datasets::*;
use stucrs::functions_new::{self as F, accuracy, sum};
use stucrs::layers::{self as L, Activation, Dense, Layer, Linear};
use stucrs::models::{BaseModel, Model};
use stucrs::optimizers::{Optimizer, SGD};

fn main() {
    let mnist = MNIST::new();
    let x_train = mnist.train_img.view();
    let y_train = mnist.train_label.view();
    let x_test = mnist.test_img.view();
    let y_test = mnist.test_label.view();

    let image_num = 0;

    let x_train = x_train.to_shape((50000, 28 * 28)).unwrap();
    let x_test = x_test.to_shape((10000, 28 * 28)).unwrap();

    // one_hotベクトル化
    let y_train = arr2d_to_one_hot(y_train.mapv(|x| x as u32).view(), 10);
    let y_test = arr2d_to_one_hot(y_test.mapv(|x| x as u32).view(), 10);

    let max_epoch = 5;
    let lr = 0.01;
    let batch_size = 100;

    let data_size = x_train.shape()[0];
    println!("data_size={}", data_size);

    let mut model = BaseModel::new();
    model.stack(L::Dense::new(1000, true, None, Activation::Sigmoid));
    model.stack(L::Linear::new(10, true, None));

    let mut optimizer = SGD::new(lr);
    optimizer.setup(&model);
    let start = Instant::now();
    for epoch in 0..max_epoch {
        let mut indices: Vec<usize> = (0..data_size).collect();
        let mut rng = thread_rng();
        indices.shuffle(&mut rng);

        let mut sum_loss = array![0.0f32];
        let mut sum_acc = 0.0f32;

        for chunk_indices in indices.chunks(batch_size) {
            let x_batch = x_train.select(Axis(0), chunk_indices).to_owned().rv();
            let y_batch = y_train.select(Axis(0), chunk_indices).to_owned().rv();

            let y = model.call(&x_batch);

            let mut loss = F::softmax_cross_entropy_simple(&y, &y_batch);
            let acc = accuracy(
                y.data().into_dimensionality().unwrap().view(),
                y_batch.data().into_dimensionality().unwrap().view(),
            );
            model.cleargrad();
            loss.backward(false);
            optimizer.update();

            let epoch_loss: Array1<f32> = (&loss.data() * (y_batch.len() as f32))
                .into_dimensionality()
                .unwrap();

            sum_loss = &sum_loss + &epoch_loss;
            sum_acc = sum_acc + acc * (y_batch.len() as f32);
        }

        let average_loss = &sum_loss / (data_size as f32);
        let average_acc = sum_acc / (data_size as f32);

        println!(
            "epoch = {:?}, train_loss = {:?}, accuracy = {}",
            epoch + 1,
            average_loss,
            average_acc
        );

        //推論
        config::set_grad_false();
        let test_data_size = x_test.shape()[0];
        let mut indices: Vec<usize> = (0..test_data_size).collect();
        let mut rng = thread_rng();
        indices.shuffle(&mut rng);

        let mut sum_loss = array![0.0f32];
        let mut sum_acc = array![0.0f32];

        for chunk_indices in indices.chunks(batch_size) {
            let x_batch = x_test.select(Axis(0), chunk_indices).to_owned().rv();
            let y_batch = y_test.select(Axis(0), chunk_indices).to_owned().rv();

            let y = model.call(&x_batch);
            let loss = F::softmax_cross_entropy_simple(&y, &y_batch);
            let acc = accuracy(
                y.data().into_dimensionality().unwrap().view(),
                y_batch.data().into_dimensionality().unwrap().view(),
            );

            let epoch_loss: Array1<f32> = (&loss.data() * (y_batch.len() as f32))
                .into_dimensionality()
                .unwrap();

            sum_loss = &sum_loss + &epoch_loss;
            sum_acc = sum_acc + acc * (y_batch.len() as f32);
        }

        let average_loss = &sum_loss / (test_data_size as f32);
        let average_acc = sum_acc / (test_data_size as f32);

        println!(
            "epoch = {:?}, test_loss = {:?}, test_accuracy = {}",
            epoch + 1,
            average_loss,
            average_acc
        );

        config::set_grad_true();
    }
}

MNIST構造体を用いてデータを読み込みます。この際、y_train、y_test のデータは1次元のラベルなので、arr2d_to_one_hot()関数で one_hotベクトル に変換します。また、MNISTデータは今までのデータ数と桁が異なるので、バッチサイズを100とします。あとは基本的に先ほどのSpiral学習 の時と同じです。model のレイヤーの構造を自由に変更して学習を色々試してみましょう。おおよそ80～90%くらいの精度になると思います。補足でも説明しましたが、他の活性化関数、例えば ReLU関数 はMNIST学習でよく使われるので、それに変更して学習するのも面白いです。


以上でStuCrsフレームワーク:基礎編 は終了です。フレームワークの基礎をこのドキュメントでは実装してきました。続いてはフレームワークのさらなる機能的な拡張、 CNN編 です。

TODO:CNN編ドキュメント、url貼り付け

高階微分

plottersを使う

MNISTのデータの用意