読者です 読者をやめる 読者になる 読者になる

MATHGRAM

主に数学とプログラミング、時々趣味について。

[Rust] 練習で最急降下法やってみた

pythonばっか触ってきた自分からすると, Rustは楽しいけど難しいです.

もっといろんな機能を使いこなせるようになりたい・・・

最急降下法

理論はクソ簡単なので, 特に説明しません. 坂道転がってくだけです.
実装では自分で微分した式を使っちゃってるので, すごいダサく仕上がってます. 将来的にはn次元空間でも使えるように, 微分の計算部分も実装したい.

zの空間は, こちらの記事をそのまま使っちゃってて,

z = 3x^2 + 5y^2 - 6xy

となっています. 一応グラフにしとくとこんな感じ.
f:id:ket-30:20170202231513p:plain:w400:h400
まだgnuplotに慣れてなくて, 汚い図です. いずれgnuplotの使い方をまとめるかも.

そして当たり前ですが

z \geq 0, \,\,z(0,0) = 0


なので, 目的となる点は(0,0)です.
何度か勾配計算 & 更新を繰り返して, 0に辿りつければ成功です.

結果 & コード

初期値は適当に(-7, 4)に落として, そのあと100回更新しました. 更新幅でbreakさせた方が良かったですね. 他の手法を実装する時そうします. f:id:ket-30:20170202233824p:plain:w400:h100
まぁとりあえずオッケー!

以下全コード

extern crate gnuplot;

use gnuplot::*;
use std::vec::Vec;

fn main() {
    // --- 最適化したい空間の設定 & 作図 ---
    let x = linspace(-3.0, 3.0, 1000);
    let y = linspace(-3.0, 3.0, 1000);

    let z = make_zspace(&x, &y);

    let mut fg = Figure::new();
    fg.axes3d()
        .surface(z.iter(), x.len(), y.len(), Some((-4.0, -4.0, 4.0, 4.0)), &[])
        .show_contours(true, true, Linear, Fix(""), Auto)
        .set_title("z space", &[]);
    fg.set_terminal("png", "zspace.png");
    fg.show();

    let mut init_p: Vec<f32> = vec![-7.0, 4.0];
    let mut gd = GradientDescent::new(init_p, 0.1);

    for i in 0..100 {
        gd.update();
        println!("{:?}", gd.p);
    }
}

// 連続なベクトルの生成
fn linspace(start: f32, limit: f32, n: usize) -> Vec<f32> {
    if n <= 1 {
        vec![limit]
    } else {
        let mut v = vec![];
            let s = start;
            let d = (limit - start) / (n - 1) as f32;
            for i in (0..n) {
                v.push(s + i as f32 * d);
            }
        v
    }
}

// 最適化したい空間の生成
fn make_zspace (x: &Vec<f32>, y: &Vec<f32>) -> Vec<f32> {
    let mut z = Vec::with_capacity((x.len()*y.len()) as usize);
    for xi in x.iter() {
        for yi in y.iter() {
            z.push( 3f32*xi.powi(2) + 5f32*yi.powi(2) - 6f32*xi*yi );
        }
    }
    z
}

struct GradientDescent {
    p: Vec<f32>,
    eta: f32,
}

impl GradientDescent {

    fn new(p: Vec<f32>, eta: f32) -> GradientDescent {
        GradientDescent{p: p, eta: eta}
    }

    fn update (&mut self, ) {
        self.p[0] -= self.eta*self.calc_grads()[0];
        self.p[1] -= self.eta*self.calc_grads()[1];
    }

    // 勾配も&で受けれるようにメンバーに加えるべきだったかも.
    // 毎回grad_vを作成しているのはセンスない.
    // 次から気をつける.
    fn calc_grads(&mut self) -> Vec<f32> {
        let mut grad_v = vec![];
        grad_v.push( 6f32 * self.p[0] - 6f32 * self.p[1]);
        grad_v.push(10f32 * self.p[0] - 6f32 * self.p[1]);

        grad_v
    }
}

まとめ

もっと, trait, structの構成方法を勉強しないとなって感じです.
かける人から見たらすごく汚く見えるんだろうなってのが自分でわかります.
もし気になったところがあったらアドバイスください. すごい喜びます.

以上です.

[Rust] zipでfor文 & 速度比較

Rustで2つのイテレータをペアで回す方法です.

二つの方法があるようなので速度比較します.

標準装備を使う

fn main() {
    let x = vec![0; 1000000];
    let y = vec![0; 1000000];
    pair_iter(x, y);
}

fn pair_iter(x: Vec<i32>, y: Vec<i32>) {
    let start = time::now();

    for (i, j) in x.iter().zip(y.iter()) {
    }
    let end = time::now();
    println!("{}", end - start);
}

出力(実行時間)

PT0.033350S

crateのitertoolsを使う

cargo.tomlに

[dependencies]
itertools = "0.4.4"

を書いてから,

extern crate itertools;
use itertools::Zip;

fn main() {
    let x = vec![0; 1000000];
    let y = vec![0; 1000000];
    pair_iter(x, y);
}

fn pair_iter(x: Vec<i32>, y: Vec<i32>) {
    let start = time::now();

    for (i, j) in Zip::new((&x, &y)) {
    }
    let end = time::now();
    println!("{}", end - start);
}

出力(実行時間)

PT0.071162S

まとめ

単にzipで回すだけならcrateを使わない方が2倍ほど早いらしい.

以上です.

[Rust] 正規分布に従う乱数を生成する [Box-Mullerの実装]

正規分布を実装したので, 正規乱数も作りたいなぁと思い実装しました. これ使ってガウスノイズ付きのデータ生成 & 線形回帰とかやっていこうかと思います.

そんなことより, ボックス=ミュラー法, 今日知りました.

うん, 普通に感動したわ.(小並感)

理論の参考はwikipediaと,

ボックス=ミュラー法(正規乱数の生成)の証明 | 高校数学の美しい物語 です.

Rustで実装する.

実装は脳死で数式を打ち込めばオッケーですね.
とりあえず正規乱数を生成する関数の部分だけ. nは個数指定のためです. 一様分布からの乱数はrand使ってんのかい!ってツッコミはなしで.

    // ボックスミュラー法で標準正規分布に従った乱数を生成する.
    fn norm_random(&mut self, n: i32) -> Vec<f64> {
        let mut y = vec![];
        let mut rng = rand::thread_rng();
        for i in 0..n {
            let u1: f64 = rng.gen_range(0f64,1f64);
            let u2: f64 = rng.gen_range(0f64,1f64);
            let a = (- 2f64 * u1.ln()).sqrt();
            let b = (2f64*f64::consts::PI*u2).cos();
            y.push( self.mu + self.sigma*a*b )
        }
        y
    }

実際に100個の乱数生成して, ヒストグラムにしてみました.

f:id:ket-30:20170202043150p:plain:w400:h400

histgramはgnuplotのboxesで書いてるんですけど, そこに渡すベクトルの作り方がくそ汚い. しかもx軸にはboxの中心になる数を与えるのにミスってるし. スマートな方法わかるまで, とりあえずこれで. (match初めてだったし, 多少はね…)

以下全コード.

fn main() {

    //平均 0, 分散 1のNormalDistribution型が生成される.
    let mut norm = NormalDistributionBilder::new().set();

    // 標準正規分布に従う乱数の生成.
    let y: Vec<f64> = norm.norm_random(100);

    // histgram生成のための下準備
    // クソコード
    let mut n: Vec<i32> = vec![0; 6];
    let m: Vec<i32> = vec![−3, -2, -1, 0, 1, 2];

    for e in &y {
        match () {
            _ if (e < &-2f64) => n[0] += 1,
            _ if (&-2f64 <= e && e < &-1f64) => n[1] += 1,
            _ if (&-1f64 <= e && e < &0f64) => n[2] += 1,
            _ if ( &0f64 <= e && e < &1f64) => n[3] += 1,
            _ if ( &1f64 <= e && e < &2f64) => n[4] += 1,
            _ if ( &2f64 <= e) => n[5] += 1,
            _ => println!("error"),
        };
    }

    //プロット
    let mut fg = Figure::new();
    fg.axes2d().boxes(&m, &n, &[Color("blue")]);
    fg.set_terminal("png", "norm_random_hist.png");
    fg.show();
}

struct NormalDistribution {
    mu : f64,
    sigma : f64,
}

impl NormalDistribution {
    fn new(mu: f64, sigma: f64) -> NormalDistribution {
        NormalDistribution{mu: mu, sigma: sigma,}
    }

    fn calc(&mut self, x: &Vec<f64>) -> Vec<f64> {
        let mut y: Vec<f64> = vec![];
        for e in x{
            let a = - (e - self.mu).powi(2) / (2.0_f64 * self.sigma.powi(2));
            let b = (2.0_f64 * f64::consts::PI).sqrt() * self.sigma;
            y.push( a.exp() / b );
        }

        y
    }

    fn show_info(&mut self){
        println!("------- Information -------");
        println!("Average: {}, Variance: {}", self.mu, self.sigma);
        println!("---------------------------");
    }

    // ボックスミュラー法で標準正規分布に従った乱数を生成する.
    fn norm_random(&mut self, n: i32) -> Vec<f64> {
        let mut y = vec![];
        let mut rng = rand::thread_rng();
        for i in 0..n {
            let u1: f64 = rng.gen_range(0f64,1f64);
            let u2: f64 = rng.gen_range(0f64,1f64);
            let a = (- 2f64 * u1.ln()).sqrt();
            let b = (2f64*f64::consts::PI*u2).cos();
            y.push( self.mu + self.sigma*a*b )
        }
        y
    }

    fn plot(&mut self, x: &Vec<f64>) {
        let mut fg = Figure::new();
        let y: Vec<f64> = self.calc(&x);
        fg.axes2d()
        .lines(x, &y, &[Color("blue")]); 
        fg.set_terminal("png", "NormalDistribution.png");
        fg.show();
    }
}

struct NormalDistributionBilder {
    mu: f64,
    sigma: f64,
}

impl NormalDistributionBilder {
    fn new() -> NormalDistributionBilder {
        NormalDistributionBilder{mu: 0f64, sigma: 1f64,}
    }

    fn mu(&mut self, mu: f64) -> &mut NormalDistributionBilder {
        self.mu = mu;
        self
    }

    fn sigma(&mut self, sigma: f64) -> &mut NormalDistributionBilder {
        self.sigma = sigma;
        self
    }

    fn set(&mut self) -> NormalDistribution {
        NormalDistribution{ mu: self.mu, sigma: self.sigma}
    }
}

次はtraitとthreadをちゃんと理解しよう.

以上です.