[Linear Analysis] 최소제곱해(Least Square solution)

최소제곱해(Least Square solution)

행렬방정식 Ax = b의 해가 없는 모순 시스템(inconsistent system)에서 그 시스템에 가장 적합한 근사해를 찾는 방법으로 최소제곱방법을 사용합니다. 다시말하면 Ax와 b와의 거리를 최소로 하는 벡터 x를 최소제곱해(least-square solution)라고 하며 식 1과 같이 나타낼 수 있습니다.

$$\tag{식 1}||b-Ax|| ≈0$$

실수 공간에 m×n 차원의 행렬 A와 벡터 b가 존재한다면 Ax = b의 최소제곱해는 식 2를 만족하는 $\hat{x}$가 됩니다. 이 식의 열공간의 차원은 n이므로 $\mathbb{R}^n$ 공간에 존재합니다.

$$\tag{식 2}||b-A \hat{x}|| \leq ||b-Ax||$$

Ax = b가 해를 가진다는 것은 선형독립으로 A는 기저벡터들로 구성됨을 의미합니다. 그러나 이미 가정한 것과 같이 이 식은 모순된 시스템으로 해를 가질 수 없습니다. 대신에 식 3에서 나타낸 것과 같이 결정된 최소제곱해와의 선형독립의 결과인 벡터 $\hat{b}$가 존재합니다.

\begin{align}\tag{식 3}Ax&\approx b\\ A\hat{x}&=\hat{b} \end{align}

식 3의 b는 식 4와 같이 $\hat{b}$와 두 벡터의 차이로 결정할 수 있습니다.

\begin{align}\tag{식 4}b&=\hat{b}+(b-\hat{b})\\ &= \hat{b}+(b-A\hat{x})\end{align}

식 4가 $\mathbb{R}^2$ 차원에서 이루어진다면 다음과 같이 시각화할 수 있습니다.

fig, ax=plt.subplots(figsize=(3,2))
nme=[r"b-A$\hat{x}$",r"$\hat{b}=A\hat{x}$", "b"]
cod=[(0,3), (2,0), (2,3)]
col=["r", "b", "g"]
for i in range(3):
    ax.arrow(0,0,cod[i][0], cod[i][1], color=col[i], lw=2, head_width=0.05)
    ax.text(cod[i][0]+0.05, cod[i][1]+0.05, nme[i], weight="bold", color=col[i])
ax.vlines(2, 0, 3, ls="dashed", color="k", alpha=0.6)
ax.spines['left'].set_position(("data", 0))
ax.spines['bottom'].set_position(("data", 0))
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_xlabel("Col A(A의 열공간)")
ax.set_xticks([])
ax.set_yticks([])
plt.show()

위 그림의 x 축은 파란색 벡터는 $\hat{b} = A\hat{x}$을 나타냅니다. 이 식은 선형독립이므로 x축은 A의 다양한 기저벡터들로 나타낼 수 있습니다. 기저벡터는 벡터의 방향을 유지하는 기본 벡터들을 나타내며 여러 벡터들로 구성되는 행렬의 기저벡터들의 집합을 열공간(Col 행렬)이라 합니다.

위 그림에서 나타낸 것과 같이 행렬 A를 구성하는 기저벡터들의 공간인 Col A와 벡터 b-A$\hat{x}$는 수직을 이루기 때문에 그들의 내적은 0이 됩니다(직교벡터(Orthogonal vectors)참조). 식 5와 같이 나타낼 수 있습니다.

$$\tag{식 5} A\cdot(b-A\hat{x})=0$$

식 5에서 A는 0행렬이 아니므로 식 6과 같이 나타낼 수 있습니다.

\begin{align}\tag{식 6} b &=A\hat{x}\\A^Tb &=A^TA\hat{x}\end{align}

식 6은 Ax = b의 최소제곱해를 계산하기 위한 식으로 Normal equation이라고 합니다.

예 1)

다음 행렬 A에 대한 행렬 방정식 Ax=b는 모순된 시스템입니다. 이 시스템의 최소제곱해?

$$A=\left[\begin{array}{rr}4&0\\0&2\\1&1\end{array}\right] \quad b=\left[\begin{array}{r}2\\0\\11\end{array}\right]$$

import numpy as np
import numpy.linalg as la
from sympy import *

A=np.array([[4, 0],[0,2],[1,1]])
b=np.array([[2],[0],[11]])
Matrix(np.c_[A, b]).rref()

(Matrix([
 [1, 0, 0],
 [0, 1, 0],
 [0, 0, 1]]),
 (0, 1, 2))

위 결과는 Ax=b를 계산한 것으로 위 결과의 3행은 $0x_1+0x_2=1$이므로 모순된 시스템(inconsistent system)입니다. 최소자승해를 계산하기 위해 식 6을 적용합니다(식 7).

\begin{align}\tag{식 7}A^TA\hat{x} &=A^Tb\\ \hat{x} &= (A^TA)^{-1}A^Tb\end{align}

식 7은 y=ax와 같은 선형결합시스템(방정식)에서 해 x를 계산하는 과정입니다. 이 과정은 np.linalg.solve(y, a) 함수를 적용하는 것으로 해를 확인할 수 있습니다.

ATA=A.T@A
ATb=np.dot(A.T, b)
hatx=la.solve(ATA, ATb)
print(hatx)

ATA=A.T@A
ATb=np.dot(A.T, b)
hatx=la.solve(ATA, ATb)
print(hatx)

$A^TA$와 같이 행렬의 전치행렬과 그 행렬의 행렬곱은 정방행렬을 생성하지만 모든 정방행렬이 하나의 해를 가지는 것은 아닙니다. 즉, 하나 이상의 다양한 해를 가질 수 있습니다. 이것은 $A^TA$가 역행렬(Inverse matrix)을 가질 수 없는 비가역행렬임을 의미합니다. 그러므로 위에서 적용한 la.solve() 함수도 에러를 발생합니다. 가역행렬은 그 행렬의 역행렬이 0이 아닌 경우 입니다.

예 2)

행렬 A에 대한 Ax=b의 최소 제곱해를 결정합니다.

$$A=\begin{bmatrix}1& 1& 0& 0\\1& 1& 0& 0\\1& 0& 1& 0\\1& 0& 1& 0\\1& 0& 0& 1\\1& 0& 0& 1 \end{bmatrix} \quad b=\begin{bmatrix}-3\\-1\\0\\2\\5\\1 \end{bmatrix}$$

A=np.array([[1,1,0,0],[1,1,0,0],[1,0,1,0],[1,0,1,0],[1,0,0,1],[1,0,0,1]])
b=np.array([-3,-1,0,2,5,1]).reshape(6,1)
ATA=np.dot(A.T, A)
ATb=np.dot(A.T, b)
la.det(ATA).round(10)

0

위에서 적용한 np.linalg.det(X) 함수는 행렬 X의 행렬식을 계산합니다. 위 결과는 $A^TA$의 행렬식이 0이므로 이 행렬은 비가역행렬입니다. 그러므로 자명하지 않은 해(다양한 해)의 존재를 의미합니다. 이 경우 기약행사다리꼴(rref)을 적용하여 해를 결정할 수 있습니다(기약행 사다리꼴 형태 (Reduced row echelon form, rref) 참조). sympy객체.rref() 함수를 사용합니다.

hatx=Matrix(np.c_[ATA, ATb]).rref()
hatx

(Matrix([
 [1, 0, 0,  1,  3],
 [0, 1, 0, -1, -5],
 [0, 0, 1, -1, -2],
 [0, 0, 0,  0,  0]]),
 (0, 1, 2))

위 결과는 행렬의 0, 1, 2 열이 피봇열이며 3열이 자유변수이므로 3열에 해당하는 변수에 따라 나머지 변수값이 결정됨을 의미합니다. 이 관계는 식 7과 같이 나타낼 수 있습니다.

\begin{align}\begin{bmatrix}1& 0& 0&1\\0& 1& 0&-1 \\0& 0& 1& -1\\0& 0& 0& 0 \end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\\x_4\end{bmatrix} &=\begin{bmatrix}x_1+x_4\\x_2-x_4\\x_3-x_4\\0\end{bmatrix}= \begin{bmatrix}-3\\-5\\-2\\0 \end{bmatrix}\\\hat{x}=\begin{bmatrix}x_1\\x_2\\x_3\\x_4\end{bmatrix}=\begin{bmatrix}3-x_4\\-5+x_4\\-2+x_4\\x_4\end{bmatrix}&=\begin{bmatrix}-3\\-5\\-2\\0 \end{bmatrix}+x_4\begin{bmatrix}-1\\1\\1\\1 \end{bmatrix}\end{align}

최소제곱해를 하나만 가지는 조건은 행렬방정식에서 유일한 해 즉, 자명한 해를 가지는 조건과 같습니다. 즉, $m \times n$의 행렬 A에서 하나의 최소제곱해를 갖는 조건은 다음과 같습니다.

• 행렬 A의 각 열공간은 선형독립입니다.
• $A^TA$는 가역행렬입니다.

모순된 행렬방정식에서 계산할 수 있는 최소제곱해와 b와의 차이 즉, $\Vert{b-A\hat{x}}\Vert$를 최소제곱오차(least square error)라 합니다.

예 3)

Ax=b의 최소제곱오차를 결정합니다.

$$A=\left[\begin{array}{rr}4&0\\0&2\\1&1\end{array}\right] \quad b=\left[\begin{array}{r}2\\0\\11\end{array}\right]$$

A=np.array([[4,0],[0,2],[1,1]])
b=np.array([2,0,11]).reshape(3,1)
ATA=np.dot(A.T, A)
ATb=np.dot(A.T, b)
la.det(ATA).round(10)

84.0

hatx=la.solve(ATA, ATb)
print(hatx)

[[1.]
 [2.]]

lse=la.norm(b-np.dot(A, hatx))
lse.round(3)

9.165

예 4)

다음 행렬방정식 Ax=b의 최소제곱해?

$$A=\begin{bmatrix}1& -6\\1& -2\\1& 1\\1& 7 \end{bmatrix} \quad b=\begin{bmatrix} -1\\2\\1\\6\end{bmatrix}$$

A=np.array([[1,-6],[1,-2],[1,1],[1,7]])
b=np.array([[-1],[2],[1],[6]])
ATA=np.dot(A.T, A)
ATb=np.dot(A.T, b)
la.det(ATA).round(10)

360.0

hatx=la.solve(ATA, ATb)
print(hatx)

[[2. ]
 [0.5]]