最小二乗法（三角行列化計算）

こんにちは。
最小二乗法は $J (\boldsymbol{x}) \triangleq \left(\boldsymbol{y} - A \boldsymbol{x} \right)^\textrm{T} \left(\boldsymbol{y} - A \boldsymbol{x} \right)$ を最小とする $\boldsymbol{x}_{\min}$ を計算します。三角行列を経るのが計算効率が良く、常道です。詳しくは下記の関係式ですが、計算を要約しますと、

1. 行列 $\tilde{A} \triangleq \begin{bmatrix} A & \boldsymbol{y} \end{bmatrix}$ （下記例では Ay）を用意し、それから（正方）上三角行列 $\tilde{R} \triangleq \begin{bmatrix} R & \tilde{\boldsymbol{y}} \\ \boldsymbol{0}^\textrm{T} & \alpha \end{bmatrix}$ （下記例では Ry）を求め、
2. これを使い、$R\,\boldsymbol{x} = \tilde{\boldsymbol{y}}$ を解いて、$\boldsymbol{x}_{\min}$ を求める。

これのPython実装例は¹、

def lsqsolv(Ay):
    n = Ay.shape[1] - 1
    Ry = triangulateupper(Ay)
    return numpy.linalg.solve(Ry[0:n,0:n], Ry[0:n,n]) # x_min

def triangulateupper(A):
    return numpy.linalg.cholesky(numpy.dot(A.transpose(),A)).transpose()

---

関係式をまとめますと（$\tilde{Q}$, $Q$ は直交行列）、

\begin{align*}
J (\boldsymbol{x}) &\triangleq \left(\boldsymbol{y} - A \boldsymbol{x} \right)^\textrm{T} \left(\boldsymbol{y} - A \boldsymbol{x} \right) \\
&= \begin{bmatrix} \boldsymbol{x}\\ -1\end{bmatrix}^\textrm{T} \tilde{A}^\textrm{T} \, \tilde{A} \,\begin{bmatrix}  \boldsymbol{x}\\ -1\end{bmatrix} \\
&= \begin{bmatrix} \boldsymbol{x}\\ -1\end{bmatrix}^\textrm{T} \tilde{R}^\textrm{T} \, \tilde{R} \,\begin{bmatrix}  \boldsymbol{x}\\ -1\end{bmatrix} \\
A &= Q R \\
\tilde{A} &= \tilde{Q} \tilde{R} \\
\tilde{R}^\textrm{T} \, \tilde{R} &= \tilde{A}^\textrm{T} \tilde{A} \\
\tilde{A} &\triangleq \begin{bmatrix} A & \boldsymbol{y} \end{bmatrix} \\
\tilde{R} &\triangleq \begin{bmatrix} R  & \tilde{\boldsymbol{y}} \\ \boldsymbol{0}^\textrm{T} & \alpha \end{bmatrix} \\
\tilde{\boldsymbol{y}} &\triangleq Q^\textrm{T} \boldsymbol{y} \\
\alpha &\triangleq \left( \min_{\boldsymbol{x}} J (\boldsymbol{x}) \right)^{1/2} \\
 \\
\end{align*}

---

1. この例の三角行列化 triangulateupper() は計算効率重視の Cholesky 分解を利用していますが、多くの場合では数値安定性重視の QR 分解をお勧めます。 ↩