多次元の最小二乗法(least squares method)についてまとめます。最小二乗法とは、観測データとの誤差が最小となるような、最も確からしい関係式を求める方法を指します。

また最後にこれらのpythonでの実装方法について記しておきます。

i) 一次関数の場合

一次関数の最小二乗法について。

$$y = ax + b$$

ここで観測データをもとに、最も誤差が小さくなるa, bの値を導く。

二乗和誤差を$E(a, b)$ (Errorの頭文字)と置く。

$$E(a,b) = \sum_{i} {(y_i - a x_i - b) ^ 2} \tag{*}$$

であるから、aとbについてそれぞれを偏微分し、その最小となる値を求めればよい。

aについての偏微分

bを固定し、aについて微分する。

$$\frac{\partial E(a,b)}{\partial a} = -2 \sum_{i} {x_i (y_i - a x_i - b)} \tag{1}$$

bについての偏微分

以下同様。

$$\frac{\partial E(a,b)}{\partial b} = -2 \sum_{i} {(y_i - a x_i - b)} \tag{2}$$

(*)の右辺はaについて見ればaの二次関数であり、また$a^2$の係数は正の実数（$x ≠ 0$）であるから、これは下に凸の放物線と言える。bについても同じことが言える。

つまりEの最小値は(1), (2)の導関数が0のときとなる。よって下記の連立方程式を解いてやれば良い。

$$\begin{equation*} \begin{cases} -2 \sum_{i} {x_i (y_i - a x_i - b)} = 0 \\ -2 \sum_{i} {(y_i - a x_i - b)} = 0 \end{cases} \end{equation*}$$

これを解いてやれば最終的に以下のような値になる。

$$a = \frac{\bar{xy} - \bar{x}\bar{y}}{\bar{x^2} - \bar{x} ^ 2}\\ b = -a \bar{x} + \bar{y}$$

ii) 行列・ベクトルの場合

多項式の場合の重みベクトルの二乗和誤差を考える。

例えば、m種類の果物が売られていてそれぞれの単価は$w_i$であり、これを$x_i$個ずつ購入する。このとき合計金額は $x_i*w_1$ の総和となる。

購入した個数と合計金額のデータから、最小二乗法を用いて果物の値段wの近似値を求める。

合計金額を多項式で表すと下記になる。

$$y = x_1 w_1 + x_2 w_2 + \dots + x_m w_m$$

次にxを変数ベクトル、wを係数ベクトルとして表現すると、

$$x = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_m \end{pmatrix} ,\ w = \begin{pmatrix} w_1 \\ w_2 \\ \vdots \\ w_m \end{pmatrix}$$

と表すことができる。

これを用いると、合計金額yはxの転置行列とwの積であり

$$y = x^{\top}w$$

と書くことができる。

ここで目標値を$t$とおくと、二乗和誤差は下記のように表せる。

$$E(w) = \frac{1}{2} \sum_{i=1}^{n} {(y_i - t_i) ^ 2}$$

さて、ここでn人の購入データ（観測データ）があるとする。

つまり、上記のベクトルは行列$\mathbb{R}^{n \times m}$として置き換えることができる。

$$y = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} ,\ X = \begin{pmatrix} x_{11} & x_{12} & \cdots & x_{1m} \\ x_{21} & x_{22} & \cdots & x_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \cdots & x_{nm} \end{pmatrix}$$

よって下記の式が成りたつ。

$$y = Xw$$

このときの二乗和誤差は、下記のようにノルムの２乗で書くことができる（$t$は目標ベクトル）。

$$E(w) = \frac{1}{2} {\left\| \mathbf{Xw - t} \right\| ^ 2} \tag{3}$$

ここで

$${\left\| \mathbf{x} \right\| ^ 2} = x^{\top} x$$

であるので、(3)は式変形して

$$\begin{align*} E(w) &= \frac{1}{2} (Xw - t)^{\top}(Xw - t) \\ &= \frac{1}{2} (X^{\top}w^{\top} - t^{\top})(Xw - t) \\ &= \frac{1}{2} (X^{\top}w^{\top}Xw - tX^{\top}w^{\top} - t^{\top}Xw + \|t\|^2) \end{align*}$$

となる。

i)一次関数で計算したように、wに関してE(w)を微分し、その導関数が0になるwを求めればよい。

先にベクトルの微分の定義を書いておく。

$$\frac{d}{dx} (x^{\top} A x) = (A + A^{\top})x \\ \frac{d}{dx} (a^{\top} x) = \frac{d}{dx} (x^{\top} a) = a$$

であるから、E(w)の微分は

$$\begin{align*} \frac{dE(w)}{dw} &= \frac{1}{2} \left( X^{\top} X + (X^{\top} X)^{\top} \right) w - X^{\top} t \\ &= X^{\top} X w - X^{\top} t \end{align*}$$

と式変形することができる。

$$\frac{dE(w)}{dw} = 0 \\ w = (X^{\top}X)^{-1}X^{\top}t$$

であるときに、E(w)は誤差が最小となる。

Pythonの実装

ii)の最小二乗法をpythonで実装すると、下記のように書ける。

import numpy as np

def predict_simple_lstsq(X, t):
    """simple least square solution"""
    w = np.linalg.inv(X.T @ X) @ (X.T @ t)
    return w

def main():
    """
    t: 目標値
    w: 重み
    X: 入力データ
    """
    t = np.array([660, 1000, 1670, 770, 890])
    X = np.array([[3, 1, 4, 1, 5],
                  [9, 2, 6, 5, 3],
                  [5, 8, 9, 7, 9],
                  [3, 2, 3, 8, 4],
                  [6, 2, 6, 4, 3]])

    w_estimated = predict_simple_lstsq(X, t)
    print(w_estimated) # [27.76613348 38.3563155  64.69387755 26.70159956 50.57363486]

if __name__ == '__main__':
    main()

補足

pythonにはlinalg.lstsq()を用いた最小二乗法の解法がある。

def predict_lstsq(X, t):
    """least square solution"""
    w = np.linalg.lstsq(X, t, rcond=None)[0]
    return w

また、重みが負の実数にならないようなNon-Negative Least Squares (NNLS)も用意されている。例えば上記の果物の例では、商品の値段が負の実数であることはあり得ないため、上述の単純な最小二乗法を用いることはできない。

from scipy.optimize import nnls

def predict_nnls(X, t):
    """non-negative least square solution"""
    w = nnls(X, t)[0]
    return w

実際に NNLS で重みを推定し、ノイズを加えた場合の誤差を比較するコードは下記。

import numpy as np
from scipy.optimize import nnls

def predict(X, y):
    # non-negative least square solution
    w = np.round(nnls(X, y)[0], 2)
    return w

def add_noise(y, noise_scale=1, noise_mean=0):
    return y + np.random.normal(noise_mean, noise_scale, len(y))

def calculate_error(w_true, w_estimated):
    return np.round(np.mean(np.abs(w_true - w_estimated)), 2)

def main():
    """
    X ∈ R^{n×m}, y ∈ R^n, w ∈ R^m
    """
    X = np.array([[3, 1, 4, 1, 5],
                  [9, 2, 6, 5, 3],
                  [5, 8, 9, 7, 9],
                  [3, 2, 3, 8, 4],
                  [6, 2, 6, 4, 3]])

    w_true = np.array([31, 41, 59, 26, 53])

    y = np.array(np.dot(X, w_true))
    y = add_noise(y)

    w_estimated = predict(X, y)
    E = calculate_error(w_true, w_estimated)

    print(f"""
    True    : {', '.join(map(str, w_true))}
    Estimates: {', '.join(map(str, w_estimated))}
    Error   : {E}
    """)

if __name__ == '__main__':
    main()