行列の最小二乗法について調べてみた

2020.07.15

カフェチームの山本です。

現在、カフェチームでは、機械学習を用いて、画像内に映っている人物の骨格や手を検出することで、どの商品を取り出したかを判定する方法を検討しています。

前回は、実際にカメラの位置や回転を計測することなく、この外部パラメータを推定する方法として、最小二乗法を用いました。

カメラの外部パラメータ行列を最小二乗法で推定する

今回は、上記の記事を作成するにあたって調べた最小二乗法について、計算過程や証明方法をまとめました。結論としては、どのような次元でも、正規方程式を用いて計算できることがわかりました。

お詫び

本ページの数式が正しく表示されていません。現在、修正中です。数式まで見たい方はこちらのページをご覧ください。

行列の最小二乗法について調べてみた

0.最小二乗法の概要

説明

例として、2次元平面において、[(x0, y0), (x1, y1), ...]のような対応関係がある、サンプル(測定点)があるとします。このサンプルを直線で近似することを考えます。すなわち、以下のようなモデルです。

\[y=qx+p\]

この近似において、「各測定点でのy座標の誤差の2乗」の和を最も小さくするような(q, p)を求めるのが最小二乗法です。これによって、測定点の近くを通る直線を求めることができます。

1.最小二乗法の定式化・計算

(1)入力:1元・出力:1元の場合

ステップ0:誤差の定式化

誤差2乗の和を定式化すると、以下のように表すことができます。

\[E = \left\{y_0-(qx_0+p)\right\}^2+\left\{y_1-(qx_1+p)\right\}^2+\cdots\]

ステップ1:行列で書き直す

これは、行列を使って以下のように書き直すことができます。

\[b = \left( \begin{matrix} y_0 \\ y_1 \\ \vdots \end{matrix} \right)\]

\[A = \left( \begin{matrix} x_0 & 1 \\ x_1 & 1 \\ \vdots \end{matrix} \right)\]

\[x = \left( \begin{matrix} q \\ p \end{matrix} \right)\]

\[E = \left\| b - Ax \right\|^2\]

(このEは、下に示すように上の式と同じです)

\[\begin{matrix} E &=& \left\| b - Ax \right\|^2 \\ &=& \left\| \left(\begin{matrix} y_0 - (qx_0+p) \\ y_1 - (qx_1+p) \\ \vdots \end{matrix}\right)\right\|^2 \\ &=& \left\{y_0-(qx_0+p)\right\}^2+\left\{y_1-(qx_1+p)\right\}^2+\cdots \end{matrix}\]

ステップ2:Eを変形させる

Eを変形させると以下のようになります。

\[\begin{matrix} E &=& \left\| b - Ax \right\|^2 \\ &=& \left(b - Ax \right)^T \left(b - Ax \right) \\ &=& \left(b^T - (Ax)^T \right) \left(b - Ax \right) \\ &=& \left(b^T - x^TA^T \right)\left(b - Ax \right) \\ &=& b^Tb- x^TA^Tb - b^TAx + x^TA^TAx \end{matrix}\]

  • 補足1:n行1列ベクトルx, yがあるとき、以下の関係があります。

\[x^Ty = \sum_ix_iy_i=y^Tx\]

そのため、上のEの変形における最後の式の2項目と3項目は同じ値です。よって下のように変形できます。

\[\begin{matrix} E &=& b^Tb- x^TA^Tb - b^TAx + x^TA^TAx \\ &=& b^Tb- 2b^TAx + x^TA^TAx \end{matrix}\]

これを各成分を代入して表すと以下のようになります。

\[\begin{matrix} E &=& b^Tb- x^TA^Tb - b^TAx + x^TA^TAx \\ &=& \sum_i b_i^2 - \sum_i y_i(qx_i+p) - \sum_i y_i(qx_i+p) + \sum_i(qx_i+p)^2 \\ &=& \sum_i b_i^2 - 2\sum_i y_i(qx_i+p) + \sum_i(qx_i+p)^2 \\ &=& \sum_i b_i^2 - 2\sum_i y_ix_iq - 2\sum_i y_ip + \sum_ix_i^2q^2 + 2\sum_i x_i qp + \sum_ip^2 \end{matrix}\]

ステップ3:偏微分・極小点(成分計算の場合)

上の式から、Eはpの2次式であり、かつ、qの2次式であることがわかります。また、pとqの2次の項が係数が正です。そのため、Eをq,pそれぞれで偏微分した値が同時に0になる点があれば、その点でEが最小になることがわかります。

\[\left(\begin{matrix} \frac{\partial E}{\partial q} \\ \frac{\partial E}{\partial p} \end{matrix}\right)=0\]

\[\Leftrightarrow \left(\begin{matrix} - 2\sum_iy_ix_i + \sum_i x_i^2 2q + 2 \sum_i x_i p \\ - 2 \sum_i y_i + 2\sum_i x_i q + \sum_i 2p \end{matrix} \right)=0\]

\[\Leftrightarrow \left(\begin{matrix} - 2 \sum_iy_ix_i + 2 \sum_i x_i^2 q + 2 \sum_i x_i p \\ - 2 \sum_i y_i + 2 \sum_i x_i q + 2 \sum_i p \end{matrix} \right)=0\]

\[\Leftrightarrow \left(\begin{matrix} - \sum_iy_ix_i + \sum_i x_i^2 q + \sum_i x_i p \\ - \sum_i y_i + \sum_i x_i q + \sum_i p \end{matrix} \right)=0\]

\[\Leftrightarrow \left( \begin{matrix} \sum_i x_i^2 q + \sum_i x_i p \\ \sum_i x_i q + \sum_i p \end{matrix} \right) = \left( \begin{matrix} \sum_iy_ix_i \\ \sum_i y_i \end{matrix} \right)\]

\[\Leftrightarrow \left( \begin{matrix} \sum_i x_i^2 & \sum_i x_i \\ \sum_i x_i & \sum_i 1 \end{matrix} \right) \left( \begin{matrix} q \\ p\end{matrix} \right) = \left( \begin{matrix} \sum_iy_ix_i \\ \sum_i y_i \end{matrix} \right)\]

左辺の第一項が正則行列である場合、

\[\Leftrightarrow \left( \begin{matrix} q \\ p\end{matrix} \right) = \left( \begin{matrix} \sum_i x_i^2 & \sum_i x_i \\ \sum_i x_i & \sum_i 1 \end{matrix} \right)^{-1}\left( \begin{matrix} \sum_iy_ix_i \\ \sum_i y_i \end{matrix} \right)\]

上のように式変形することで、(q,p)を求めることができます。

ステップ3:偏微分・極小点(行列計算の場合)

これをxの各成分に関して偏微分すると以下のようになります。

\[\nabla = \left( \begin{matrix} \frac{\partial}{\partial q} \\ \frac{\partial}{\partial p} \end{matrix} \right)\]

\[\begin{matrix} \nabla E &=& - 2 \left( b^T A\right)^T + \left( A^T A + (A^TA)^T \right)x \\ &=& - 2 A^T b + \left( A^T A + A^T A \right)x \\ &=& - 2 A^T b + 2 A^T A x \end{matrix}\]

Eが最小になるには、偏微分の値が0になることが必要条件です。

\[\nabla E = 0\]

\[\Leftrightarrow - 2 A^T b + 2 A^T A x = 0\]

\[\Leftrightarrow 2 A^T A x = 2 A^T b\]

\[\Leftrightarrow A^T A x = A^T b\]

という式が導けます。A^{T} Aが正則行列のとき、

\[x = \left( A^T A \right)^{-1} A^T b\]

として、xを計算することができます。これは、上の成分を代入して計算した場合と同じです。

  • 補足3:A^{T}Aが正則行列がどうかの判定

A^{T}Aが正則行列であることは、以下のように考えられます。

正則行列の性質から、以下が成り立ちます。

\[A^{T}Aが正則である\Leftrightarrow \rm{rank}(A^{T}A)=2\]

ランク(階数)の性質から、以下が成り立ちます

\[\rm{rank}(A^{T}A)=2 \Leftrightarrow \rm{rank}(A)=2\]

また、階数の定義から、以下が成り立ちます。

\[\rm{rank}(A)=2 \Leftrightarrow Aの行ベクトルにおいて、線形独立なものが2つ以上存在する\]

Aは計測したサンプルのx座標を並べた行列です。なので、上と組み合わせると、x座標の異なる点で2つ以上サンプルを計測すれば、A^{T}Aが正則行列になることがわかります。

(つまり、2点以上測定して直線を引けるようにすれば良いということです)

(2)入力:多元・出力:1元の場合

前節は入力が1元(x)、出力が1元(y)でしたが、入力が多元になっても同様に計算できます。

以下では、入力をx,yとして、出力をzとします

誤差2乗の和を定式化すると、以下のようです。

\[E = \left\{z_0-(ry_0 + qx_0+p)\right\}^2+\left\{z_1-(ry_1+qx_1+p)\right\}^2+\cdots\]

これは、行列を使って以下のように書き直すことができます。

\[b = \left( \begin{matrix} z_0 \\ z_1 \\ \vdots \end{matrix} \right)\]

\[A = \left( \begin{matrix} y_0 & x_0 & 1 \\ y_1 & x_1 & 1 \\ \vdots \end{matrix} \right)\]

\[x = \left( \begin{matrix} r \\ q \\ p \end{matrix} \right)\]

とすると、Eは以下のようになり、前節と同じ形になります。

\[E = \left\| b - Ax \right\|^2\]

よって、Eを最小にするxは、同様に計算できます。

\[x = \left( A^T A \right)^{-1} A^T b\]

(3)入力:多元・出力:多元の場合

出力が多元になっても同様の計算ができます。

以下では、入力を(x,y,z)、出力を(x',y',z')とします。

誤差2乗の和を定式化すると、以下のようです。

\[E = \left\{z_0'-(r_{11}z_0 + r_{21}y_0 + r_{31}x_0 +r_{41})\right\}^2 + \left\{y_0'-(r_{12}z_0 + r_{22}y_0 + r_{32}x_0 + r_{42})\right\}^2 + \left\{x_0'-(r_{13}z_0 + r_{23}y_0 + r_{33}x_0 + r_{43})\right\}^2 + \left\{z_1'-(r_{11}z_1 + r_{21}y_1 + r_{31}x_1 +r_{41})\right\}^2 + \left\{y_1'-(r_{12}z_1 + r_{22}y_1 + r_{32}x_1 + r_{42})\right\}^2 + \left\{x_1'-(r_{13}z_1 + r_{23}y_1 + r_{33}x_1 + r_{43})\right\}^2 + \cdots\]

これは、行列を使って以下のように書き直すことができます。

\[B = \left( \begin{matrix} x_0' & y_0' & z_0' \\ x_1' & y_1' & z_1' \\ & \vdots & \end{matrix} \right) = \left( \begin{matrix} b_x & b_y & b_z \end{matrix} \right)\]

\[A = \left( \begin{matrix} z_0 & y_0 & x_0 & 1 \\ z_1 & y_1 & x_1 & 1 \\ & \vdots & & \end{matrix} \right)\]

\[X = \left( \begin{matrix} r_{11} & r_{12} & r_{13} \\ r_{21} & r_{22} & r_{23} \\ r_{31} & r_{32} & r_{33} \\ r_{41} & r_{42} & r_{43} \end{matrix} \right)= \left( \begin{matrix} r_1 & r_2 & r_3 \end{matrix} \right)\]

このとき、Eをx,y,zそれぞれに関する誤差でまとめると、以下のようになります。

\[E = \left\| b_x - A r_1 \right\|^2 + \left\| b_y - A r_2 \right\|^2 + \left\| b_z - A r_3 \right\|^2\]

Eの形をみると、それぞれの項でr1,r2,r3が独立しており、かつ、各項を最小化されればEが最小化されます。よって、(今回の目的である)Eを最小化するr1,r2,r3を求めるには、各項を最小化するr1,r2,r3をそれぞれ求めれば良いことになります。そしてこれは、前節と同じ計算です。

\[r_1 = \left( A^T A \right)^{-1} A^T b_x\]

\[r_2 = \left( A^T A \right)^{-1} A^T b_y\]

\[r_3 = \left( A^T A \right)^{-1} A^T b_z\]

3つの式を横にならべると

\[\left( \begin{matrix} r_1 & r_2 & r_3 \end{matrix} \right) = \left( A^T A \right)^{-1} A^T \left( \begin{matrix} b_x & b_y & b_z \end{matrix} \right)\]

\[\Leftrightarrow X= \left(A^T A \right)^{-1} A^T B\]

となり、前節と同様の形になります。

まとめ

サンプルとなる対応関係が得られたとき、(線形な)モデルでそれらを近似する方法である、最小二乗法についてまとめました。最小二乗法では、入力出力が多次元の場合でも利用することができます。具体的には、対応関係の行列をA,Bとしたとき、以下のように計算できます。

\[x = \left( A^T A \right)^{-1} A^T b\]

\[X= \left(A^T A \right)^{-1} A^T B\]

参考にさせていただいたページ

最小二乗法の行列表現(単回帰,多変数,多項式)

正規方程式の導出と計算例

最小二乗法の基礎をまとめる - エンジニアを目指す浪人のブログ

行列のランクの定義・例・性質/公式 (証明付) - 理数アラカルト -

補足

これを利用して、カメラの外部パラメータの推定を行った様子を、以下の記事で書いています。

カメラの外部パラメータ行列を最小二乗法で推定する