機械学習とCS : 所信表明と、手始めの線形回帰のRustによる実装

本ブログの所信表明

こんにちは、名古屋でデータサイエンティストとして働いているボブと申します。日々、機械学習モデルの業務実装に向けて、地道に頑張ってます。一方、趣味でコンピュータサイエンス(CS)を学んでおり、現在はRustとコンピュータアーキテクチャに興味を持ってます。

今後、機械学習モデルの開発・適用の観点でCSを学び、発信していきたいと思います。

本記事では、ブログ執筆の所信表明と、手始めに線形回帰のRustによる実装を共有します。

当面の目論見 : レイヤーを降りていく学習

コンピュータサイエンスの領域の裾野は途方もなく広く、僕には全体像をイメージすることすらできません。実務に活きるコンピュータサイエンスを学びたいのですが、あまりにも領域が広すぎます。
そこで、現時点の目論見では、機械学習モデル、しかも今僕が使っているTransformersを頂点として、そこからトップダウンで降りていく形で勉強・実験・共有をしたいと考えています。大雑把に高レイヤーから順に低レイヤーまでを分けると以下の具合かなと

• 機械学習モデル
• アルゴリズムとデータ構造
• 線形代数等の数学
• プログラミング言語
• コンパイラ、OS
• ハードウェア

当然、これらトピックそれぞれが非常に重厚なものなので、各レイヤーについてまずはインデックスをはることを目指し、その後徐々に掘り下げていきたいです。

また、実際にCSを機械学習の社会実装に活かすために、応用を睨んだ以下のトピックにも取り組みたいと思います

• 数値実験
• プロファイリング（ボトルネックの特定）
• 並列計算・GPU利用
• AIチップの活用

なお、これらテーマは現時点のとても浅い理解によるものであり、今後学習していくにつれて大きく見直す予定です。

なぜRust?

本ブログでは、実装にあたってはRustを主として取り上げたいと考えています。なぜRsutかというと

• 低レイヤーまでにらめるプログラミング言語を題材にしたいが、システムプログラミングまで担える言語は、実質C++とRustの2択と言われている

• 安全に並列プログラミングを実装できると言われている

• これ↓

  If you can use rust, ignore carbon

ただ、そのうち既存のライブラリをみたり触る必要が出てくる可能性もあり、C++を拒絶するつもりはないです（風呂敷広げすぎたくはないですが。。。）

線形回帰をRustで実装してみた

上述した目論見のためには、とにかく実用に近いシチュエーションで色々試して遊ぶことが不可欠です。そのため、今後いくつかの機械学習モデルをRustで実装し、実験やプロファイリング等を行うオモチャにしたいと思います。

本記事では以下、手始めに「はむかずさん」の著書『機械学習のエッセンス』を参考に、線形重回帰を実装したいと思います。

線形回帰とは

詳細な解説は他書に譲りますが、ここでは実装にあたっての最低限の振り返りをしておきます。

線形回帰は、目的変数$y$を説明変数$x$の線型結合

$$\hat{y}_{i} = \sum_{j=1}^{p} w_j x_{ij} \tag{1}$$

により予測するモデルです。なお上記の記法では、$x_0 = 1$であり、$\beta_0$は定数項を示しています。$\beta_j$はモデルの推定（学習）対象のパラメータであり、

$$J(\bold{w}) = \sum_{i=1}^{n} (y_i - w_0 - \sum_{j=1}^{p} w_j x_{ij})^2 \tag{2}$$

の最小化により推定されます。

上記はベクトルおよび行列表記により、以下で現されます

$$\hat{\bold{y}} = \bold{\~{X}}\bold{w} \tag{3}$$

ここで$\bold{\~{X}}$は特徴量行列$\bold{X}$（1行が一つのレコードの特徴量ベクトルで、縦にデータ数$n$個だけ特徴量ベクトルがならんだもの）にたいして、定数項を表現するために左側にすべて値が$1$の列ベクトルを追加したもので、デザイン行列と呼ばれます。

通常、機械学習モデルは勾配降下法を用いて反復的に学習されますが、線形回帰モデルについては以下の式により解析的に最適パラメータが算出されます（導出は他書に譲ります）

$$\bold{w} = \left( \bold{\~{X}}^T \bold{\~{X}} \right)^{-1} \bold{\~{X}}^T \bold{y} \tag{4}$$

実装する処理

模擬的な正解データに対し、2変数の重回帰分析を以下の通り実装

• 2列（$x_{1}$, $x_{2}$）×100行の0〜1の乱数行列$\bold{X}$を生成し、すべての値をscale倍（本件ではscale=10）する（uniform_random_matrix()）

• $\bold{X}$に対して、左側にすべて値が$1$の列ベクトルをconcatしてデザイン行列$\bold{\~{X}}$を定義（design_matrix()）

• 予測すべき回帰係数$\bold{w}$を$(1, 2, 3)$として、以下から正解データを作成

$$y_{正解値} = 1 + 2 * x_{1} + 3 * x_{3} \tag{5}$$

• 正解データにノイズを加える（add_noise()、$\epsilon$は標準正規分布に従う乱数）

$$y_{正解値} = y_{正解値} + \epsilon \tag{6}$$

• 一般的に逆行列を解くより一次方程式を解く方が高速なので、$\bold{w}$の解析解を上記の式(4)で求めるかわりに、以下の一次方程式を解くことで$\bold{w}$を求める

$$\bold{\~{X}}^T \bold{\~{X}} \bold{w} = \bold{\~{X}}^T \bold{y} \tag{7}$$

• 最後に、$\bold{w}$を表示し、意図したとおり$(1, 2, 3)$に近い値となっているか確認

実装してみた

今回、まずは超ナイーブに実装しました。今後リファクタリング・拡張して徐々に成長させます。

まず依存関係です。

[package]
name = "ridge_naive"
version = "0.1.0"
edition = "2021"

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

[dependencies]
nalgebra = "0"
rand ="0"
rand_distr = "0"

次にコードです。
なお、ベクトル・行列・線形代数まわりについてはnalgebraを使ってますが、他にもより良い実装がありうると思っており、今後調査したいです。

main.rs

use nalgebra::{DMatrix, DVector, SVD};
use rand::distributions::{Distribution, Uniform};
use rand::thread_rng;
use rand_distr::{StandardNormal};

fn uniform_random_matrix(dim: usize, n: usize, scale: f64) -> DMatrix<f64> {
    // 0から1までの一様乱数生成器
    let mut rng = thread_rng();
    let uniform = Uniform::new(0.0, 1.0);
    // 一様乱数の特徴量行列
    DMatrix::from_fn(
        n, dim,
        |_, _| uniform.sample(&mut rng) * scale
    )
}

fn design_matrix(x: &DMatrix<f64>) -> DMatrix<f64> {
    DMatrix::from_fn(
        x.nrows(), x.ncols() + 1 ,|row, col| {
            if col == 0 { 1.0 }
            else { x[(row, col - 1)] }
        }
    )
}

fn add_noise(vector: &mut DVector<f64>) {
    // 標準正規分布
    let normal = StandardNormal;

    // 乱数生成器
    let mut rng = thread_rng();

    // 各要素にノイズを加える
    for i in 0..vector.nrows() {
        let noise: f64 = normal.sample(&mut rng);
        vector[i] += noise;
    }
}

fn main() {
    //-------------------------------------------------------------------
    // 諸元の設定
    //-------------------------------------------------------------------
    // 説明変数の数
    let dim: usize = 2;
    // データ数
    let n: usize = 100;
    // 特徴量のスケール
    let scale = 10.0;
    
    //-------------------------------------------------------------------
    // データの生成
    //-------------------------------------------------------------------
    // 一様乱数の特徴量行列
    let x = uniform_random_matrix(dim, n, scale);
    // デザイン行列
    let x_til = design_matrix(&x);

    // 回帰係数の正解値
    let w = DVector::from_vec(
        vec![1.0, 2.0, 3.0]
    );
    // 目的変数（行列積で1行のDMatrixができるので、DVectorに変換）
    let mut y = DVector::from_column_slice(
        ((&x_til * &w)).column(0).as_slice()
    );
    // 目的変数にノイズを加える
    add_noise(&mut y);

    //-------------------------------------------------------------------
    // 回帰分析
    //-------------------------------------------------------------------
    let svd = SVD::new(&x_til.transpose() * &x_til, true, true);
    y = DVector::from_column_slice(
        (x_til.transpose() * y).column(0).as_slice()
    );
    let w_pred = svd.solve(&y, 1.0e-8).unwrap();
    println!("w_pred: {:?}", w_pred);

}

出力は以下のとおりで、意図したとおり$(1, 2, 3)$に近い値となっています

w_pred: VecStorage { data: [1.3246357365512231, 1.995832938791781, 2.9915401004188014], nrows: Dyn(3), ncols: Const }

今後の予定

上記は超ナイーブな実装なので、以下のような観点で機能追加・リファクタリングします。

• 機能追加
  • csvデータを読みとってデータセットとして用いる
  • Ridge回帰、Lasso回帰追加
• リファクタリング
  • 学習器を構造体にし、学習と予測をメソッド化
  • エラー処理の実装
  • テストプログラムの追加

上記の機能追加の後、実験やプロファイリングのオモチャとして今後使っていきたいと思います。