[ML] Constrained Least Squares: Lagrangian 을 활용.

문제설정

흔히 볼 수 있는 least squares problem은 다음과 같음.

$$\mathbf{A'f'}=\mathbf{b}$$

 

이를 다음과 같이 augmented matrix와 homogeneous coordinate를 사용하여 homogeneous equation형태로 정리할 수 있음.

$$ \mathbf{A} = \begin{bmatrix} \mathbf{A'} & - \mathbf{b} \end{bmatrix} \\ \mathbf{f}= \begin{bmatrix} \mathbf{f'} \\ 1 \end{bmatrix}$$

$$ \mathbf{Af}=\mathbf{0} $$

위와 같이 homogeneous equation를 이용하여 푸는 방법은 least square로써 
augmented matrix $\begin{bmatrix} \textbf{A} & \textbf{b} \end{bmatrix}$ 를 사용하는 Total Least Squares 와 구분되므로 주의할 것.

2024.06.22 - [Programming/ML] - [Fitting] Total Least Square Regression

[Fitting] Total Least Square Regression

 

 

위의 문제는 noise 등으로 0이 되지는 못하는 over-determined 이기 때문에,

다음과 같은 형태의 최소화 문제로 풀게 된다. 많은 경우 미지수들의 벡터에 다음과 같은 constraint가 붙는다.

(이런 형태는 computer vision 등에서 많이 다루게 됨.)

 

$$\underset{\mathbf{f}}{\text{argmin }} ||\mathbf{A}\mathbf{f}||^2$$
subject to:

• $||\mathbf{f}||^2 = 1$ 임

여기서 $\mathbf{A}$는 matrix이고, $\mathbf{f}$는 구하고자 하는 미지수들의 벡터임.

$||\mathbf{f}||^2 = 1$는 $\mathbf{f}$의 유클리드 노름이 1이라는 제약 조건을 의미함.

 

예를 들면, Stereo calibration에서 다른 각도로 촬영한 두 이미지의 corresponding features를 통해 fundamental matrix를 구하는 경우, $\mathbf{f}$는 fundamental matrix $\mathbf{F}$를 flatten한 vector로 처리 될 수 있음.

---

Sol.

1. Object function와 Constraint 설정
   • Object function:
     • $\underset{\mathbf{f}}{\text{argmin }}||\mathbf{Af}||^2$
   • Constrints:
     • $||\mathbf{f}||^2 = \mathbf{f}^\top \mathbf{f} = 1$
2. Lagrangian 설정
   Constrained Optimization Problem를 라그랑주 승수법(Lagrange multiplier)을 사용하여 unconstrained optimization problem으로 바꾸어 푼다.
   
   Lagrangian는 다음과 같이 정의됨:
   $$
   \mathcal{L}(f, \lambda) = ||\mathbf{Af}||^2 + \lambda (||\mathbf{f}||^2 - 1)
   $$
   여기서 $\lambda$는 Lagrange Multiplier임.
   2023.06.26 - [.../Math] - [Math] Lagrange Method or Lagrange Multiplier Method
3. Lagrangian 미분
   Optimal solution를 찾기 위해 $\mathcal{L}$을 $\mathbf{f}$에 대해 편미분하고 이를 0으로 설정함:
   $$
   \frac{\partial \mathcal{L}}{\partial \mathbf{f}} = 2\mathbf{A}^\top \mathbf{A} \mathbf{f} + 2\lambda \mathbf{f} = 0
   $$
   이를 정리하면 다음과 같은 형태가 됨:
   $$
   \mathbf{A}^\top \mathbf{Af} + \lambda \mathbf{f} = 0
   $$
4. 위 식을 일반화된 eigen value problem로 변환
   위 식은 다음과 같은 일반화된 eigen value problem으로 해석할 수 있음:
   $$
   \mathbf{A}^\top \mathbf{A f} = -\lambda \mathbf{f}
   $$
   여기서 $\lambda$는 $\mathbf{A}^\top \mathbf{A}$의 고유값(eigenvalue)이고,
   $\mathbf{f}$는 대응하는 고유벡터(eigenvector)임.
5. 고유값 문제 해결
   행렬 $\mathbf{A}^\top \mathbf{A}$의 고유값과 고유벡터를 계산함.
6. $\lambda$가 가장 작은 고유값일 때, 대응하는 고유벡터 $\mathbf{f}$를 선택함.
   이 고유벡터가 주어진 조건을 만족하는 $\mathbf{f}$가 됨.
7. Norm이 1이 되도록 정규화
   고유벡터 $\mathbf{f}$는 일반적으로 노름이 1이 아닐 수 있음. 따라서, 최종 결과로 얻은 $\mathbf{f}$를 정규화하여 노름이 1이 되도록 함:
   $$
   \mathbf{f}_{\text{normalized}} = \frac{\mathbf{f}}{||\mathbf{f}||}
   $$

2024.11.06 - [.../Linear Algebra] - [LA] Eigenvalue and Eigenvector

[LA] Eigenvalue and Eigenvector

 

---

요약

1. 행렬 $\mathbf{A}^\top \mathbf{A}$를 계산함.
2. $\mathbf{A}^\top \mathbf{A}$의 고유값과 고유벡터를 구함.
3. 가장 작은 고유값에 대응하는 고유벡터를 선택함.
4. 선택한 고유벡터를 정규화하여 노름이 1이 되도록 함.

이 과정을 통해 주어진 제약 조건을 만족하는 벡터 $\mathbf{f}$를 구할 수 있음.

---

다음은 Homogeneous Least Squares와 Total Least Squares를 비교한 예임.

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import svd

# Generate 10 samples in 2D with noise
np.random.seed(0)
X = np.random.rand(10, 1) * 10  # X values between 0 and 10
true_slope = 2.0
true_intercept = 1.0
noise = np.random.randn(10, 1)
Y = true_slope * X + true_intercept + noise  # Y values with noise

# Linear regression using homogeneous least squares
A_ls = np.hstack([X, np.ones((10, 1))])
A_homogeneous = np.hstack((A_ls, -Y))

# Perform SVD on the homogeneous augmented matrix
U_hom, S_hom, Vt_hom = svd(A_homogeneous)
f_homogeneous = Vt_hom[-1]
f_homogeneous = f_homogeneous / f_homogeneous[-1]  # normalize so the last entry is 1

# Extract f from f_homogeneous
f_homogeneous = f_homogeneous[:-1]

# Linear regression using Total Least Squares
X_mean = np.mean(X)
Y_mean = np.mean(Y)
X_centered = X - X_mean
Y_centered = Y - Y_mean
A_tls = np.hstack([X_centered, Y_centered])
U_tls, S_tls, Vt_tls = svd(A_tls)
v_tls = Vt_tls[-1]

# Calculate the slope and intercept for TLS
slope_tls = -v_tls[0] / v_tls[1]
intercept_tls = Y_mean - slope_tls * X_mean

# Plotting the results and error directions with x-axis range 0 to 10 and extending the lines over the entire range
x_range = np.linspace(0, 10, 100)

# Line for Homogeneous Least Squares
y_homogeneous_line = f_homogeneous[0] * x_range + f_homogeneous[1]

# Line for Total Least Squares
y_tls_line = slope_tls * x_range + intercept_tls

plt.figure(figsize=(14, 6))

# Plotting Homogeneous Least Squares result
plt.subplot(1, 2, 1)
plt.scatter(X, Y, color='blue', label='Data with noise')
plt.plot(x_range, y_homogeneous_line, color='red', label='Homogeneous Least Squares Fit')
for i in range(len(X)):
    plt.plot([X[i, 0], X[i, 0]], [Y[i, 0], Y_pred_homogeneous[i]], 'k--')
plt.title('Homogeneous Least Squares Fit')
plt.xlabel('X')
plt.ylabel('Y')
plt.xlim(0, 10)
plt.legend()

# Plotting Total Least Squares result
plt.subplot(1, 2, 2)
plt.scatter(X, Y, color='blue', label='Data with noise')
plt.plot(x_range, y_tls_line, color='red', label='Total Least Squares Fit')
for i in range(len(X)):
    # Calculate the projection of the point (X[i], Y[i]) onto the TLS line
    x_proj = (X[i, 0] + slope_tls * (Y[i, 0] - intercept_tls)) / (1 + slope_tls**2)
    y_proj = slope_tls * x_proj + intercept_tls
    plt.plot([X[i, 0], x_proj], [Y[i, 0], y_proj], 'k--')
plt.title('Total Least Squares Fit')
plt.xlabel('X')
plt.ylabel('Y')
plt.xlim(0, 10)
plt.legend()

plt.show()

# Error calculation
error_homogeneous, error_tls

---

같이보면 좋은 자료

2022.04.28 - [Programming/ML] - [Fitting] Ordinary Least Squares : OLS, 최소자승법

[Fitting] Ordinary Least Squares : OLS, 최소자승법