MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

スタンフォード大学の学生が学ぶ最小二乗法入門

 

目次

 

はじめに

スタンフォード大学には

機械学習を学ぶ上での第一歩として、

Introduction to Matrix Methods (EE103)という授業があります。

 

今回の記事では、

この授業の教科書である

Introduction to Applied Linear Algebraを

読んだ際の技術メモです。

 

この教科書は下記のリンクのページから

pdfをダウンロードすることができます。

 

本記事では、

上記の教科書の最小二乗法の部分のみのメモです。

他の部分に関しては、下記の記事を参照下さい。

myenigma.hatenablog.com

myenigma.hatenablog.com

 

最小二乗法

最小二乗法は、方程式の近似解を計算する方法です。

方程式の誤差の二乗値の和を最小化するような解を計算するため、

最小二乗法と呼ばれます。

一般的に最小二乗法というと、線形方程式の近似解を求める方法を指し、

線形ではない方程式の場合は非線形最小二乗法と呼ばれます。

 

今回の記事では、線形方程式における最小二乗法について説明します。

ある線形方程式 Ax = bを解きたいとします。

ここで、Aの次元をm x nとすると、解を求めるには

m ≧ nである必要があります。

m=n である場合は、行列の逆行列を計算することで

x = A^-1bから解を計算することができます。

 

ここで Ax = bの解は、

Ax - bの要素ができるだけ0に近ければ

より良い解であるとすると、

r = Ax -bとし、このrの二乗値の和が最小化されれば良いといえます。

このような方法による方程式の解の求め方を最小二乗法といいます。

 

Juliaにおける最小二乗法の解き方

Julia言語で最小二乗法を解く方法はいくつか存在します。

||Ax- b||^2を最小化したい場合を考えると、

下記の4つの方法があります。

xhat = inv(A’*A)*(A’*b) #通常の逆行列invを使う方法
xhat = pinv(A)*b # 疑似逆行列pinvを使う方法
Q,R = qr(A); xhat = inv(R)*(Q’*b) #QR分解を使う方法
xhat = A\b #バックスラッシュ演算子を使う方法 

以上の4つとも同じ結果が得られるはずです。

 

最小二乗法の応用例1: 広告の投資最適化

ある広告媒体が3つあったとします。(例: TV, 雑誌, ネットなど)

それぞれの広告媒体において、

読者の年齢を10等分した各グループ(20歳~70歳の5歳刻みのグループなど)に対する

ある一定の費用を投資した時に、得られる読者の数のデータがあったとします。

例えば、広告1では、100万円投資すると、20-25歳の読者は0.97人得られ、

25-30歳の読者は1.23人得られるなどです。

このデータをそれぞれの広告と、各グループを行と列にすると、

下記のような行列が得られます。

R = [
    0.97 1.89 0.41;
    1.23 2.18 0.53;
    0.80 1.24 0.62;
    1.29 0.98 0.51;
    1.10 1.23 0.69;
    0.67 0.34 0.54;
    0.87 0.26 0.62;
    1.10 0.16 0.48;
    1.92 0.22 0.71;
    1.29 0.12 0.62;
        ]

ここで、各グループにおいて

ある一定数の数の読者を得たい場合、

各広告にどれだけ投資すれば良いでしょうか?

これはシンプルな最小二乗法で解くことができます。

例えば、

すべての年齢のグループから1000人は見てもらい、

また55歳 - 65歳の人には少し多めに見てもらいたいとします。

そのような場合は、下記のような各グループの目標読者のベクトルを作り、

v = [1000 1000 1000 1000 1000 1000 1000 1300 1500 1000]'

||Rs - v||^2が最小化されるような投資ベクトルsを計算すると、

s_hat = inv(R'*R)*(R'*v)

下記のような結果が得られます。

f:id:meison_amsl:20180722164921p:plain

これは上記の値に100万円掛けた値を

各広告に投資すれば良いということを意味しています。

 

最小二乗法の応用例2:ライトの光量の制御

ある部屋に複数のライトが配置されている時、

部屋の各場所で明かりを一定にしたい場合にも、

最小二乗法を利用することができます。

 

まず、下記のように赤いバツ印の位置にライトがある場合、

すべてのライトの光量を一定にすると、

ライトがまとまっている場所が明るくなってしまい、

ライトが無い場所は暗くなってしまいます。

f:id:meison_amsl:20180722224408p:plain

一方、下記のコードのように、

各グリッドを各ライトまでの距離を格納したものを行列Rにし、

各グリッドの目標光量のベクトルlをもとに、

||Rx - l||^2を最小化することで、

各グリッドの光量を目標値に近づける

各ライトの光量を計算することができます。

using PyPlot

function draw_heatmap(data, minx, maxx, miny, maxy, xyreso)
    u = [i for i in (minx-xyreso/2.0):xyreso:(maxx+xyreso/2.0)-xyreso]
    v = [i for i in (miny-xyreso/2.0):xyreso:(maxy+xyreso/2.0)-xyreso]

    x   = [ u1+xyreso for u1 in u, v1 in v]
    y   = [ v1+xyreso for u1 in u, v1 in v]

    maxvalue = 0.0
    for i in data
        if i == Inf
            continue
        elseif i > maxvalue 
            maxvalue = i
        end
    end
    axis("equal")

    pcolor(x, y, data, vmax=maxvalue)
end

function create_matrix(np, n, lights)
    A = []

    for i in 1:np
        Ai = Float64[]
        for ix in 1:n
            for iy in 1:n
                d = hypot(ix-lights[i,1], iy-lights[i,2]) 
                push!(Ai, d^-1)
            end
        end

        if length(A) == 0
            A = Ai
        else
            A = hcat(A, Ai)
        end
    end

    return A
end


function main()
    nmap = 20 # size of grid map
    nlight = 30 # size of grid map

    lights = rand((nlight,2)).*(nmap-1).+1.0

    A = create_matrix(nlight, nmap, lights)

    # w/o optimization    
    p = fill(0.5, nlight)
    data = A * p
    data = reshape(data, (nmap,nmap))'
    draw_heatmap(data, 0, nmap, 0, nmap, 1.0)
    plot(lights[:,1], lights[:,2], "xr")
    show()

    # w/ optimization    
    l = fill(1.0, nmap*nmap)
    p_hat = inv(A'*A)*(A'*l)
    println(p_hat)
    data = A * p_hat
    data = reshape(data, (nmap,nmap))'
    draw_heatmap(data, 0, nmap, 0, nmap, 1.0)
    plot(lights[:,1], lights[:,2], "xr")
    show()
end

main()

下記が計算されたライトの光量から、

計算された各グリッドの光量です。

f:id:meison_amsl:20180722225347p:plain

前述の結果と比べると、

各グリッドの光量のばらつきが

小さくなっていることがわかります。

 

多目的関数における最小二乗法

前述の線形最小二乗法は、

一つの線形方程式の誤差を最小化するものでしたが、

下記の式のように、複数の線形方程式を最小化したい場合もあります。

f:id:meison_amsl:20180805145553p:plain

これはそれぞれのの線形方程式の誤差に重みλを設定し、

それぞれの重み付き足し算の結果を最小化するというものです。

 

このような最小二乗法では、各重みによって

どの線形方程式の誤差を最小化するか

トレードオフとして決定することができます。

大抵の場合、初めの項の重みをλ1 = 1とし、

この線形方程式を主目的関数(primary objective)とし、

残りの線形方程式の重みを調整することになります。

 

上記のような、複数目的関数の最小二乗法を解く場合は、

下記の式のように、各目的関数に重みのルートをかけ合わせて、

それぞれの項を連結して、一つの線形方程式を作成し、

あとは前述の通常の最小二乗法と同じ方法で解くことで、

多目的関数の場合も解くことが可能です。

f:id:meison_amsl:20180805151733p:plain

 

ちなみに前述の各重み係数λは、線形ではなく、

下記のような対数スペースにおけるサンプリングをしたほうが、調整がしやすい場合が多いようです。

f:id:meison_amsl:20180805152325p:plain

f:id:meison_amsl:20180805152336p:plain

 

Tikhonovの正規化による最小二乗法

ある線形方程式の近似解を計算したいが、

できるだけxのノルムを小さくしたい場合、

下記のような評価式を最小化することになります。

f:id:meison_amsl:20180811115510p:plain

 

これをTikhonovの正規化といいます。

これはxが、システムの入力になっており、

線形方程式を解きたいが、

入力を最小化してエネルギーを最小化したい場合などに使われます。

上記の式は、複数目的関数における最小二乗法になるため、

前述の方法で解くことができます。

下記のような式で計算できます。

f:id:meison_amsl:20180811115543p:plain

 

多目的関数最小二乗法による時系列データフィッテング

時系列のフィッテングをする場合にも、

多目的関数の最小二乗法は有効です。

例えば、時系列のデータでかつ、

周期的なデータを線形方程式でフィッテングしたい場合、

以前の記事で紹介した手法では、時系列間の関係はコスト関数に入っていないため、

元のデータに大きなノイズが含まれる場合、

近似された結果がなめらかではなくなってしまう場合があります。

そのような場合、下記のような多目的関数の最小二乗法を使うことが可能です。

左の項は通常の線形最小二乗法ですが、

右側の項は、Dcircという行列を使って、

隣り合うxの値の差が小さくなるようなコスト関数が追加されています。

f:id:meison_amsl:20180811115919p:plain

f:id:meison_amsl:20180811115929p:plain

これにより、λの値が大きい場合、隣り合うxの値の差が小さくなり、

よりスムーズな近似結果が得られます。

 

画像のボケ除去

前述の多目的関数の最小二乗法と同じ考えを利用することにより、

画像のボケの除去も可能です。

例えば、下記のような最小二乗法を考えます。

xは画像のベクトル情報、yは元画像とした時に、

DhとDvが横方向と縦方向の差分行列です。

この二番目の項の効果で画像がなめらかに補正されます。

f:id:meison_amsl:20180811120126p:plain

上記のような最小二乗法で画像を補正しつつ、

λの値を調整すると下記のようなボケを除去した結果が得られます。

f:id:meison_amsl:20180811120152p:plain

CTスキャンの断面復元

医療分野でよく使われるCTスキャンにおいても、

この多目的関数の最小二乗法はよく使用されます。

それぞれの位置で照射されたレーザの反射データを、

多目的関数の最小二乗法を使うことでノイズ除去をし、

そこから人体の断面図などを復元しているようです。

トモグラフィック復元 - Wikipedia

多目関数によるオーバーフィッティング防止

多目的関数による最小二乗法を使うことで、

データフィッティングの過学習を防止することもできます。

これはde-tuningやshrinkage, ridge 回帰と呼ばれる方法です。

詳細は下記の記事を参照ください。

myenigma.hatenablog.com

 

ちなみに、正規化項のラムダの値は、下記の図のように、

横軸のラムダを様々な値でふり、

テストデータのRMS誤差が最小になるようなラムダを設定します。

f:id:meison_amsl:20180812121335p:plain

 

また過学習を防ぐためには、

モデルパラメータの設計(Feature engineering)も重要であり、

おおよそ、学習データの数の10%以下に、

モデルパラメータの数を減らす必要があるようです。

 

GItHubリポジトリ

前述のJuliaのコードは下記のリポジトリで公開しています。

github.com

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

 

MyEnigma Supporters

もしこの記事が参考になり、

ブログをサポートしたいと思われた方は、

こちらからよろしくお願いします。

myenigma.hatenablog.com