# 일반화 최소제곱법(가중치 적용)

#### 개요와 기본 아이디어

최소제곱법은 관측값과 모델로부터 예측된 값 사이의 오차 제곱합을 최소화하는 문제이다. 관측 오차가 모두 동일한 분산을 갖는다고 가정하면, 표준 최소제곱법으로도 충분하다. 그러나 어떤 관측점들은 측정 과정에서 다른 점보다 더 정확하거나, 혹은 관측에 들어간 비용이나 중요도가 다르다면, 모든 관측점을 똑같이 취급하는 것보다는 상대적 중요도를 반영할 필요가 있다. 이러한 가정을 일반화한 것이 가중치 적용 최소제곱법이다.

가중치 적용 최소제곱법에서는 관측값 중 상대적으로 신뢰도가 높은(오차 분산이 작은) 점에 더 큰 가중치를 부여하고, 신뢰도가 낮은(오차 분산이 큰) 점에는 더 작은 가중치를 부여한다. 이렇게 하면 전체적으로 오차가 클 수 있는 데이터의 영향은 줄어들고, 오차가 작은(또는 잘 측정된) 데이터가 해의 결정에 더 큰 영향을 미치게 된다.

#### 일반화 최소제곱 문제의 정의

표준 선형 최소제곱법과 유사하게, 일반화 최소제곱법에서도 목표는 모델 파라미터를 찾는 것이며, 모델의 형태는 선형 또는 비선형일 수 있다. 선형 모델부터 간단히 살펴보면, 독립변수 벡터를 $\mathbf{x}\_i$라 하고, 종속변수(관측값)를 $y\_i$라 하면, 어떤 파라미터 벡터 $\mathbf{a}$에 대해 예측값 $\hat{y}\_i = \mathbf{x}\_i^T \mathbf{a}$가 주어진다고 하자. (비선형 모델인 경우 $\hat{y}\_i = f(\mathbf{x}\_i, \mathbf{a})$의 형태가 될 것이다.)

관측 오차로 인한 잔차(residual)는

$$
r\_i = y\_i - \hat{y}\_i
$$

로 정의할 수 있다. 가중치 적용 최소제곱법에서는 오차의 제곱합 대신, 가중치 행렬 $\mathbf{W}$를 고려한 다음 식을 최소화한다.

$$
\begin{align} S(\mathbf{a})  &= \mathbf{r}^T \mathbf{W} \mathbf{r}\\
&= \sum\_{i=1}^n w\_i \bigl(y\_i - \hat{y}\_i \bigr)^2 \end{align}
$$

여기서 $\mathbf{r}$은 잔차 벡터

$$
\mathbf{r} = \begin{bmatrix} r\_1 \ r\_2 \ \vdots \ r\_n \end{bmatrix}
$$

을 의미하고, $\mathbf{W}$는 대각 성분에 가중치가 들어 있는 (또는 좀 더 일반적으로 분산 공분산 행렬의 역행렬) 대칭 양정치 행렬이다. 가중치가 단순히 점마다 다른 불확실도의 역수라면

$$
\mathbf{W} = \begin{bmatrix} w\_1 & 0 & \cdots & 0 \ 0 & w\_2 & \cdots & 0 \ \vdots & & \ddots & \vdots \ 0 & 0 & \cdots & w\_n \end{bmatrix}
$$

처럼 대각 행렬로 놓을 수 있다. $w\_i = \frac{1}{\sigma\_i^2}$ 같은 형태로, 관측값의 오차 분산을 $\sigma\_i^2$라 할 때 $\sigma\_i^2$가 작은 관측값에 더 큰 가중 $w\_i$가 부여된다.

#### 선형 모델에서의 해법

가중치가 주어진 선형 최소제곱법은 행렬 $\mathbf{X}$를 사용하면 더욱 간단히 표현된다. 독립변수 벡터 $\mathbf{x}\_i$를 행으로 가지도록 구성한 행렬을

$$
\mathbf{X} = \begin{bmatrix} - (\mathbf{x}\_1)^T - \ - (\mathbf{x}\_2)^T - \ \vdots \ - (\mathbf{x}\_n)^T -  \end{bmatrix}
$$

라 하고, 종속변수 벡터를

$$
\mathbf{y} = \begin{bmatrix} y\_1 \ y\_2 \ \vdots \ y\_n \end{bmatrix}
$$

라 하면, 잔차 벡터는

$$
\mathbf{r} = \mathbf{y} - \mathbf{X}\mathbf{a}
$$

가 된다. 그러면 최소화해야 할 목적함수는

$$
S(\mathbf{a}) = (\mathbf{y} - \mathbf{X}\mathbf{a})^T \mathbf{W} (\mathbf{y} - \mathbf{X}\mathbf{a})
$$

가 된다. 이 식을 $\mathbf{a}$에 대해 미분하고 0으로 놓으면 정상방정식이 얻어진다.

$$
\begin{align} \frac{\partial}{\partial \mathbf{a}} S(\mathbf{a}) &=  \frac{\partial}{\partial \mathbf{a}} \bigl\[(\mathbf{y} - \mathbf{X}\mathbf{a})^T \mathbf{W} (\mathbf{y} - \mathbf{X}\mathbf{a})\bigr]\\
&=  -2 \mathbf{X}^T \mathbf{W} (\mathbf{y} - \mathbf{X}\mathbf{a})  = \mathbf{0} \end{align}
$$

이를 정리하면 다음의 (가중치가 반영된) 정상방정식을 얻는다.

$$
\mathbf{X}^T \mathbf{W} \mathbf{X} ,\mathbf{a} = \mathbf{X}^T \mathbf{W} ,\mathbf{y}
$$

이 식이 비가역이 아닌 경우(즉 $\mathbf{X}^T \mathbf{W} \mathbf{X}$가 역행렬을 갖는 경우), 해는

$$
\mathbf{a} =  \bigl(\mathbf{X}^T \mathbf{W} \mathbf{X}\bigr)^{-1} \mathbf{X}^T \mathbf{W} \mathbf{y}
$$

가 된다. 이를 통해 데이터의 상대적 신뢰도를 반영하면서 선형 모델 파라미터를 추정할 수 있다.

#### 가중치 행렬과 공분산 행렬

가중치 행렬 $\mathbf{W}$는 단순히 각 점의 상수 가중만을 의미하지 않고, 더 일반적으로는 오차의 공분산 행렬의 역행렬로 간주할 수도 있다. 만약 오차 항들이 서로 독립적이고 각 점마다 분산이 $\sigma\_i^2$라면,

$$
\mathbf{V} = \begin{bmatrix} \sigma\_1^2 & 0 & \cdots & 0 \ 0 & \sigma\_2^2 & \cdots & 0 \ \vdots & & \ddots & \vdots \ 0 & 0 & \cdots & \sigma\_n^2 \end{bmatrix}
$$

가 오차 공분산 행렬이 된다. 이때

$$
\mathbf{W}  =  \mathbf{V}^{-1}
$$

가 되므로,

$$
w\_i = \frac{1}{\sigma\_i^2}
$$

의 형태를 갖는다.

만약 오차들이 서로 상관되어 있는 경우(즉 공분산 행렬이 대각이 아닌 항도 갖는 경우), $\mathbf{W}$는 일반적인 대칭 양정치 행렬이 되며, 이 행렬을 통해 상관된 오차까지 반영하여 파라미터를 추정할 수 있다. 이 경우에도 동일한 방식으로 정상방정식을 세울 수 있으며, 해를 구할 때

$$
(\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1}
$$

가 존재해야 한다.

#### 비선형 모델에서의 가중치 적용

모델이 비선형인 경우, 예측값 $\hat{y}\_i = f(\mathbf{x}\_i,\mathbf{a})$라면, 최소화해야 할 가중 잔차 제곱합은

$$
S(\mathbf{a}) = \sum\_{i=1}^n w\_i \bigl\[y\_i - f(\mathbf{x}\_i,\mathbf{a})\bigr]^2
$$

이 된다. 일반화하면

$$
S(\mathbf{a}) = (\mathbf{y} - \mathbf{f}(\mathbf{a}))^T \mathbf{W} (\mathbf{y} - \mathbf{f}(\mathbf{a}))
$$

로 표현할 수 있다. 여기서 $\mathbf{f}(\mathbf{a})$는 각 관측점에서의 모델 예측값을 모은 벡터다.

비선형 문제의 경우, 미분 방정식을 통해 직접 해를 구하기보다는, 반복 알고리즘(예: 가중 Gauss-Newton법, Levenberg-Marquardt 등)을 통해 해를 찾게 된다. 가중치가 들어가는 경우에도 잔차 벡터와 야코비 행렬을 적절히 정의하고, 각 반복 단계에서 해를 갱신하는 방식으로 접근한다. 이 과정에서 잔차에 곱해지는 $\mathbf{W}$가 미분 과정에 함께 들어가게 되어, 일반 Gauss-Newton 방정식도 조금 변형된다. 예를 들어 $k$번째 단계에서의 잔차를 $\mathbf{r}(\mathbf{a}^{(k)}) = \mathbf{y} - \mathbf{f}(\mathbf{a}^{(k)})$라 하고, 야코비 행렬을 $\mathbf{J}(\mathbf{a}^{(k)})$라 하면, 다음 방향 벡터를 얻어 해를 갱신한다.

#### 고급 주제: 오차 추정과 공분산

일반화 최소제곱법에서 얻은 해 $\mathbf{a}$는 표준편차와 같은 불확실성을 함께 추정할 수 있다. 선형인 경우 $\mathbf{a}$의 공분산 행렬은

$$
(\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1}
$$

에 분산 추정값을 곱하여 구할 수 있다. 비선형의 경우에도, 수치적으로 선형화(1차 근사)하여 비슷한 형태의 공분산 추정값을 얻는다. 이는

$$
\hat{\mathbf{V}}\_{\mathbf{a}}  \approx \bigl(\mathbf{J}^T \mathbf{W} \mathbf{J}\bigr)^{-1} \sigma^2
$$

의 형태를 갖게 되며, $\sigma^2$는 잔차로부터 추정된 오차 분산(Scale factor)이다.

***

가중치 적용 최소제곱법은 데이터의 신뢰도나 오차 특성이 균일하지 않을 때, 혹은 오차가 서로 상관되어 있을 때, 보다 합리적인 추정을 가능케 해준다. 선형 모델일 경우에는 가중치 행렬을 사용한 정상방정식을 직접 풀어서 해를 얻을 수 있고, 비선형 모델에서는 반복 알고리즘을 통해 잔차를 점진적으로 줄여나가면서 해를 구한다.

#### 가중 Gauss-Newton 알고리즘

비선형 가중치 최소제곱법에서 자주 활용되는 방법은 가중 Gauss-Newton 알고리즘이다. 모델이 $f(\mathbf{x}\_i,\mathbf{a})$로 주어졌을 때, 각 관측점에서의 잔차를

$$
r\_i(\mathbf{a}) = y\_i - f(\mathbf{x}\_i,\mathbf{a})
$$

라 하고, 이들을 벡터 형태로 모으면

$$
\mathbf{r}(\mathbf{a}) = \begin{bmatrix} r\_1(\mathbf{a}) \ r\_2(\mathbf{a}) \ \vdots \ r\_n(\mathbf{a}) \end{bmatrix} = \mathbf{y} - \mathbf{f}(\mathbf{a})
$$

가 된다. 여기에서 $\mathbf{f}(\mathbf{a})$는

$$
\begin{bmatrix} f(\mathbf{x}\_1,\mathbf{a}) \ f(\mathbf{x}\_2,\mathbf{a}) \ \vdots \ f(\mathbf{x}\_n,\mathbf{a}) \end{bmatrix}
$$

로 구성된다.

가중치를 고려하면 최소화해야 할 함수는

$$
S(\mathbf{a}) = \mathbf{r}(\mathbf{a})^T \mathbf{W} , \mathbf{r}(\mathbf{a})
$$

가 된다. 이때 $\mathbf{W}$는 (대칭 양정치) 가중치 행렬로, 흔히 오차 공분산 행렬의 역행렬 형태를 가진다.

가중 Gauss-Newton 알고리즘은 잔차 벡터의 1차 테일러 전개를 이용하여, 현재 추정값 $\mathbf{a}^{(k)}$ 근방에서

$$
\mathbf{r}(\mathbf{a}) \approx \mathbf{r}(\mathbf{a}^{(k)}) + \mathbf{J}(\mathbf{a}^{(k)})  \bigl(\mathbf{a} - \mathbf{a}^{(k)}\bigr)
$$

로 근사한다. 여기서 야코비(편미분) 행렬 $\mathbf{J}(\mathbf{a}^{(k)})$는

$$
\mathbf{J}(\mathbf{a}^{(k)}) = \begin{bmatrix} \frac{\partial r\_1}{\partial a\_1} & \frac{\partial r\_1}{\partial a\_2} & \cdots & \frac{\partial r\_1}{\partial a\_m} \ \frac{\partial r\_2}{\partial a\_1} & \frac{\partial r\_2}{\partial a\_2} & \cdots & \frac{\partial r\_2}{\partial a\_m} \ \vdots & & \ddots & \vdots \ \frac{\partial r\_n}{\partial a\_1} & \frac{\partial r\_n}{\partial a\_2} & \cdots & \frac{\partial r\_n}{\partial a\_m} \end{bmatrix}
$$

로 구성된다.

이 근사를 이용해 $S(\mathbf{a})$를 선형화하면

$$
S(\mathbf{a}) \approx \bigl( \mathbf{r}(\mathbf{a}^{(k)}) + \mathbf{J}(\mathbf{a}^{(k)}) ,\Delta \bigr)^T \mathbf{W} \bigl( \mathbf{r}(\mathbf{a}^{(k)}) + \mathbf{J}(\mathbf{a}^{(k)}) ,\Delta \bigr)
$$

가 된다. 여기서 $\Delta = \mathbf{a} - \mathbf{a}^{(k)}$는 파라미터의 변화량이다. 이 식을 $\Delta$에 대해 최소화하면, 다음과 같은 방정식이 나온다.

$$
\mathbf{J}(\mathbf{a}^{(k)})^T \mathbf{W} , \mathbf{J}(\mathbf{a}^{(k)}) , \Delta = \mathbf{J}(\mathbf{a}^{(k)})^T \mathbf{W} , \mathbf{r}(\mathbf{a}^{(k)})
$$

여기서 만약 $\mathbf{J}(\mathbf{a}^{(k)})^T \mathbf{W} , \mathbf{J}(\mathbf{a}^{(k)})$가 역행렬을 가진다면,

$$
\Delta = \bigl\[ \mathbf{J}(\mathbf{a}^{(k)})^T \mathbf{W} , \mathbf{J}(\mathbf{a}^{(k)})  \bigr]^{-1} \mathbf{J}(\mathbf{a}^{(k)})^T \mathbf{W} , \mathbf{r}(\mathbf{a}^{(k)})
$$

의 형태로 $\Delta$를 구할 수 있다. 그리고

$$
\mathbf{a}^{(k+1)} = \mathbf{a}^{(k)} + \Delta
$$

로 갱신한다. 이 과정을 잔차가 충분히 작아지거나, 변화량 $\Delta$가 매우 작아질 때까지 반복한다.

이 알고리즘은 잔차가 비교적 작고, 모델이 충분히 선형 근사 가능한 영역에서 수렴이 빠르다. 가중치를 반영함으로써, 잔차 제곱합을 실제 오차 분산의 역수로 스케일링하여 더 정확한 해를 추정할 수 있다.

#### 가중 Levenberg-Marquardt 알고리즘

비선형 문제가 매우 민감하거나, 잔차가 큰 초깃값 부근에서 Gauss-Newton법이 불안정하게 움직일 경우, Levenberg-Marquardt(이하 LM) 알고리즘을 활용할 수 있다. 이 또한 가중치가 있는 형태로 쉽게 일반화할 수 있다. 일반적 LM 알고리즘은 Gauss-Newton 업데이트 식에 견인 계수 $\lambda$를 추가하여 다음과 같이 수정한다.

$$
\bigl\[ \mathbf{J}(\mathbf{a}^{(k)})^T \mathbf{W} , \mathbf{J}(\mathbf{a}^{(k)})  + \lambda ,\mathbf{I} \bigr] ,\Delta = \mathbf{J}(\mathbf{a}^{(k)})^T \mathbf{W} , \mathbf{r}(\mathbf{a}^{(k)})
$$

이 식을 풀어 $\Delta$를 구하고,

$$
\mathbf{a}^{(k+1)}  =  \mathbf{a}^{(k)} + \Delta
$$

로 갱신한다. $\lambda$는 음이 아닌 실수로, 작을 때는 Gauss-Newton에 가깝고, 클 때는 Gradient Descent(경사하강)에 가까워진다. 이 방법은 Gauss-Newton에 비해 안정적으로 수렴 경로를 탐색하는 장점이 있다. 가중치가 있는 경우에도 동일한 형태의 방정식에 $\mathbf{W}$가 포함된다.

#### 가중 최소제곱법에서의 정규화(Regularization)

오차가 상관되어 있거나, 모델의 민감도가 높아 해가 불안정해질 우려가 있다면, Tikhonov(티코노프) 정규화나 Ridge 회귀(선형의 경우)를 추가하는 방법도 고려된다. 예를 들어 선형 모델에서, 가중 최소제곱에 정규화 항 $\alpha |\mathbf{L}\mathbf{a}|^2$를 첨가하면, 최소화해야 할 목적함수는

$$
S(\mathbf{a}) = (\mathbf{y} - \mathbf{X}\mathbf{a})^T \mathbf{W} , (\mathbf{y} - \mathbf{X}\mathbf{a}) + \alpha , |\mathbf{L}\mathbf{a}|^2
$$

와 같다. 이때 $\mathbf{L}$은 원하는 형태의 규제(예: 단순 L2, 미분 연산자 등)를 반영할 수 있다. 이를 해석하면

$$
\begin{align} \frac{\partial}{\partial \mathbf{a}} S(\mathbf{a}) &= -2 \mathbf{X}^T \mathbf{W} (\mathbf{y} - \mathbf{X}\mathbf{a}) + 2,\alpha,\mathbf{L}^T \mathbf{L},\mathbf{a} = \mathbf{0}
\Longrightarrow (\mathbf{X}^T \mathbf{W},\mathbf{X} + \alpha,\mathbf{L}^T \mathbf{L}),\mathbf{a} = \mathbf{X}^T\mathbf{W},\mathbf{y}. \end{align}
$$

여기서

$$
\mathbf{a} = (\mathbf{X}^T \mathbf{W} \mathbf{X} + \alpha,\mathbf{L}^T \mathbf{L})^{-1} \mathbf{X}^T \mathbf{W} ,\mathbf{y}
$$

의 형태로 정규화된 해를 얻을 수 있다. 비선형 문제에서도 마찬가지 아이디어를 이용해, 반복 단계별로 정규화 항을 포함한 Gauss-Newton 혹은 LM 알고리즘을 적용한다.

#### 간단한 Python 예시

비선형 함수 $f(x,a) = a\_1 e^{a\_2 x}$가 있다고 가정하고, 어떤 관측값들이 주어졌을 때 이를 가중 최소제곱으로 피팅한다고 해보자. 예시 코드를 간단히 제시하면 다음과 같다.

```python
import numpy as np

# 독립변수 x, 종속변수 y, 오차분산 sigma^2가 각각 주어졌다고 가정
x_data = np.array([...])
y_data = np.array([...])
sigma_data = np.array([...])  # 각 점마다 분산

# 가중치 행렬을 대각 행렬로 가정
W = np.diag(1.0 / sigma_data**2)

# 초기 추정값
a = np.array([1.0, 0.1])

def model(x, params):
    return params[0] * np.exp(params[1]*x)

def residuals(a):
    return y_data - model(x_data, a)

def jacobian(a):
    # 야코비: d/d(a_1), d/d(a_2)
    J = np.zeros((len(x_data), 2))
    J[:,0] = np.exp(a[1]*x_data)       # df/da_1
    J[:,1] = a[0]*x_data*np.exp(a[1]*x_data)  # df/da_2
    return J

for iteration in range(100):
    r = residuals(a)
    J = jacobian(a)
    grad = J.T @ W @ r
    H = J.T @ W @ J
    # Gauss-Newton update
    delta = np.linalg.inv(H) @ grad
    a_new = a + delta
    if np.linalg.norm(a_new - a) < 1e-12:
        break
    a = a_new

print("Estimated parameters:", a)
```

이 예시는 Gauss-Newton 형태의 업데이트를 직접 구현한 것이며, 가중치가 있는 경우 $W = \mathrm{diag}(1/\sigma\_i^2)$와 같이 각 점의 분산 정보에 기반을 두어 잔차를 평가한다.

#### 가중 Total Least Squares

표준 최소제곱법(Ordinary Least Squares, OLS)에서는 독립변수(또는 디자인 행렬) $\mathbf{X}$에 오차가 없고, 종속변수(관측값) $\mathbf{y}$에만 오차가 있다고 가정한다. 그러나 실제로는 독립변수도 측정 과정에서 오차가 발생할 수 있다. 이처럼 $\mathbf{X}$와 $\mathbf{y}$ 모두가 오차를 포함할 수 있는 경우, 총체적 최소제곱법(Total Least Squares, TLS)이 유용하다.

기본 아이디어는

$$
\min\_{\Delta \mathbf{X}, ,\Delta \mathbf{y}, ,\mathbf{a}} \bigl|\Delta \mathbf{X}\bigr|\_F^2 + \bigl|\Delta \mathbf{y}\bigr|^2
$$

의 형태로, 행렬 $\mathbf{X}$와 벡터 $\mathbf{y}$에 대한 수정량($\Delta \mathbf{X}$, $\Delta \mathbf{y}$) 및 해 $\mathbf{a}$를 동시에 찾아서,

$$
(\mathbf{X} + \Delta \mathbf{X})\mathbf{a}  \approx  (\mathbf{y} + \Delta \mathbf{y})
$$

를 만족시키게 한다. $|\cdot|\_F$는 Frobenius 노름이고, 여기서는 $\mathbf{X}$의 행별 혹은 열별 오차 모두를 고려하게 된다.

만약 자료에 따라 특정 행(혹은 관측점)마다 서로 다른 신뢰도가 있다면, 가중 행렬을 활용한 Total Least Squares로 확장할 수 있다. 이를 위해

$$
\mathbf{W} = \mathrm{diag}(w\_1,\dots,w\_n)
$$

이 주어진다면,

$$
\min\_{\Delta \mathbf{X}, ,\Delta \mathbf{y}, ,\mathbf{a}} \sum\_{i=1}^n w\_i \Bigl| \begin{bmatrix} \Delta \mathbf{X}\_{i,\cdot} & \Delta y\_i \end{bmatrix} \Bigr|^2
$$

와 같은 형태로 문제를 정의할 수 있다. 이 문제는 SVD(특잉값 분해)나 사영 연산 등을 응용하여 해를 구할 수 있다. 계산 과정이 일반적인 최소제곱법에 비해 복잡하지만, 독립변수와 종속변수 모두의 오차를 균형 있게 반영하기 때문에 고정밀 분석이 필요한 상황에서 사용된다.

#### 가중치가 있는 로버스트(robust) 회귀

최소제곱법은 제곱 오차를 최소화하는 데 초점을 맞춘다. 이로 인해 이상값(outlier)이나 노이즈 분포의 꼬리가 두꺼운 경우, 해가 크게 왜곡될 수 있다. 이러한 문제를 완화하기 위해 로버스트 회귀가 도입되며, Huber, Cauchy, Lorentzian 등 다양한 로버스트 손실 함수를 사용한다.

가중치가 있는 로버스트 회귀에서는 보통 잔차에 대한 가중뿐 아니라, 로버스트 가중(이상값일 가능성이 높은 점에 대한 작은 가중)을 추가적으로 적용한다. 예를 들어 관측점 $i$에 대한 초기 가중 $w\_i$가 있고, 로버스트 오차 함수에서 추가로 구한 내부 가중 $\psi\_i$가 있다면, 최종적으로

$$
w\_i^\prime = w\_i \cdot \psi\_i
$$

의 형태로 결합 가중을 구성할 수 있다. 이 결합 가중을 행렬 형태로 모아 $\mathbf{W}^\prime$라 하면, 앞서 살펴본 가중치 적용 최소제곱 문제의 해법을 그대로 적용할 수 있다.

#### 가중 다항식 근사

데이터가 다항식 형태로 잘 설명된다고 가정하면, 예를 들어 $p$차 다항식

$$
f(x)  = a\_0  +  a\_1 x  +  a\_2 x^2  + \cdots + a\_p x^p
$$

로 모델링할 수 있다. 관측값 $(x\_i, y\_i)$, 가중치 $w\_i$가 주어졌다고 할 때, 최소화해야 할 함수는

$$
S(a\_0,a\_1,\dots,a\_p) = \sum\_{i=1}^n  w\_i \bigl\[ y\_i  - f(x\_i) \bigr]^2.
$$

이를 행렬로 표현하면

$$
\mathbf{X} = \begin{bmatrix} 1 & x\_1 & x\_1^2 & \cdots & x\_1^p \ 1 & x\_2 & x\_2^2 & \cdots & x\_2^p \ \vdots & \vdots & & \ddots & \vdots \ 1 & x\_n & x\_n^2 & \cdots & x\_n^p \end{bmatrix},  \quad \mathbf{a} = \begin{bmatrix} a\_0 \ a\_1 \ \vdots \ a\_p \end{bmatrix}, \quad \mathbf{y} = \begin{bmatrix} y\_1 \ y\_2 \ \vdots \ y\_n \end{bmatrix}.
$$

또한, 가중 행렬을

$$
\mathbf{W} = \begin{bmatrix} w\_1 & 0 & \cdots & 0 \ 0 & w\_2 & \cdots & 0 \ \vdots & & \ddots & \vdots \ 0 & 0 & \cdots & w\_n \end{bmatrix}
$$

로 놓으면, 가중 다항식 근사의 해는

$$
\mathbf{a} = (\mathbf{X}^T  \mathbf{W} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W} \mathbf{y}
$$

가 된다. (물론 $\mathbf{X}^T \mathbf{W} \mathbf{X}$가 역행렬을 가져야 한다는 전제가 필요하다.) 다항식 차수가 높아지면 과적합(overfitting) 문제가 발생할 수 있으므로, 필요한 경우 정규화(regularization)를 적용하거나 적절한 차수를 선택하여 해를 안정화한다.

#### 가중 최소제곱과 직교다항식

다항식 근사를 할 때 기저 함수를 단순히 ${1, x, x^2, \dots }$로 놓으면 수치적으로 불안정해질 수 있다. 차수가 커지면 항들 사이의 상관성이 높아지고, $\mathbf{X}^T \mathbf{W} \mathbf{X}$가 ill-conditioned 상태에 놓이기 쉽다. 이를 개선하기 위해 직교다항식(예: Legendre, Chebyshev)을 기저로 활용하면 수치 안정성을 높일 수 있다.

가중이 포함된 상황에서도 직교다항식을 구성할 수 있는데, 예를 들어 구간 $\[a,b]$ 위에 정의된 $p$차 직교다항식 집합 ${P\_k(x)}$와 특정 가중함수 $\omega(x)$를 이용해,

$$
\int\_a^b P\_k(x)  P\_l(x)  \omega(x) , dx = \delta\_{kl},
$$

(크론네커 델타) 형태를 만족하도록 설계한다. 이 경우, 데이터가 $x\_i \in \[a,b]$에 분포해 있고, 점별로 측정된 가중 $w\_i$가 있다면, 이를 근사적으로 변형하거나 샘플링 기반의 내적 정의에 녹여서 직교성을 유지하는 방법도 가능하다. 고차 다항 근사 시, 수치 정밀도나 계산 효율 면에서 유리한 기법이다.

#### 대규모 문제에서의 가중 최소제곱(스파스 구조)

데이터의 개수가 매우 많거나, 설계 행렬 $\mathbf{X}$가 희소(sparse) 구조를 가진 대규모 문제에서, 가중 최소제곱을 직접 풀기 위해서는

$$
(\mathbf{X}^T \mathbf{W} \mathbf{X}) ,\mathbf{a} = \mathbf{X}^T \mathbf{W} ,\mathbf{y}
$$

를 푸는 과정이 병목이 될 수 있다. 이때 스파스 선형 해법(CG, GMRES 등)이나 스파스 LU 분해, 혹은 특잉값 분해(SVD) 알고리즘의 희소 버전을 사용하여 연산량을 줄일 수 있다.

비선형 문제인 경우, 가중 Gauss-Newton이나 LM 알고리즘 단계마다 등장하는 야코비 행렬이 큰 차원을 가질 수 있으므로, 스파스 구조를 활용하면 연산 효율을 높일 수 있다. 실제로 수많은 관측점(예: 3D 스캐닝, 대규모 센서 네트워크)에서, 측정값 간의 상호 상관(오차 공분산)까지 고려하면, 대각이 아닌 위치에도 요소가 발생할 수 있으나 대부분의 항이 0인 경우가 많다. 이를 적절히 관리하면 일반 PC에서도 대규모 문제를 실용적으로 다룰 수 있다.

#### 상태공간, 칼만필터(Kalman Filter)와의 관련성

동적 시스템에서 시간에 따라 관측이 이루어지고, 각 시점마다의 측정에 대해 가중치(또는 공분산)를 적용한 최소제곱 문제를 반복적으로 푼다고 생각해보면, 실제로 이것은 선형 칼만필터(Kalman Filter)와 동일한 수식 구조를 갖는다. 칼만필터는 상태추정 문제에서 예측-갱신을 반복하는 알고리즘으로, 가중치 행렬(오차 공분산의 역행렬)을 고려한 선형화된 최소제곱 문제를 시계열 방향으로 푸는 것으로 해석할 수 있다.

따라서 동적 시스템의 파라미터 추정을 다룰 때, 가중 최소제곱법을 수시로 업데이트하는 접근이 곧 칼만필터(또는 확장 칼만필터, EKF)의 핵심 아이디어와 직결된다. 비선형 모델에 대해서도 확장 칼만필터나 Unscented 칼만필터(UKF) 등을 적용하는 방식으로 연결된다.

#### 가중치 선택의 실제적인 고려사항

이론적으로는 관측값마다 추정된 분산(또는 표준편차)의 역수를 가중치로 사용하는 것이 가장 직관적이다. 그러나 실제 데이터 분석에서 각 점의 오차 분산을 미리 알고 있는 경우는 드물다. 따라서 다음과 같은 상황별 접근이 가능하다.

오차 분산을 사전에 어느 정도 추정할 수 있는 경우, 예를 들어 센서 장비의 스펙이나 반복 측정 등을 통해 $y\_i$의 표준편차 $\sigma\_i$를 측정했다면,

$$
w\_i  = \frac{1}{\sigma\_i^2}
$$

로 설정한다. 만약 다중 센서가 서로 다른 정밀도를 가진다면, 센서 종류별로 $\sigma\_i$를 달리 설정하여 가중치를 배분할 수 있다. 측정 불확실도에 대한 사전 정보가 전혀 없다면, 먼저 표준(비가중) 최소제곱법으로 추정한 뒤 잔차들의 통계를 바탕으로 각 점의 가중치를 갱신하는 과정(Iterative Reweighted Least Squares, IRLS)을 수행하기도 한다. IRLS는 로버스트 회귀 방정식을 풀 때도 자주 이용되며, 각 단계에서 새로운 가중치를 잔차 크기에 대한 함수를 통해 업데이트한다.

예컨대, (가중치 없는) 최소제곱으로 추정한 해 $\mathbf{a}^{(0)}$를 가지고 잔차 벡터 $\mathbf{r}^{(0)}$를 구한다. 어떤 규칙(예: 잔차의 크기가 크면 상대적으로 작은 가중)을 통해

$$
w\_i^{(1)} = g\bigl( r\_i^{(0)} \bigr)
$$

의 형태로 가중치를 다시 설정한 다음, 이를 사용하여 가중 최소제곱을 풀어 새로운 해 $\mathbf{a}^{(1)}$를 구한다. 그 후 다시 $\mathbf{r}^{(1)}$을 구해 $w\_i^{(2)}$를 결정하는 식으로 반복한다. 이 과정을 적절히 수렴할 때까지 진행하면, 이상값(outlier)의 영향을 줄이면서도 각 점의 측정 오차 규모를 고려하는 해를 얻을 수 있다.

#### 베이지안적 관점과 가중 최소제곱

가중 최소제곱법은 베이지안 추정 관점에서도 자연스럽게 해석된다. 관측 오차가 정규분포를 따르며, 각각 평균 0, 분산 $\sigma\_i^2$라고 가정하면, 사후 확률을 최대화(MAP 추정)하는 문제는

$$
p(\mathbf{a} \mid \mathbf{y}) \propto p(\mathbf{y} \mid \mathbf{a}) ,p(\mathbf{a})
$$

의 최대화가 된다. $p(\mathbf{y} \mid \mathbf{a})$가

$$
\prod\_{i=1}^n \exp\Bigl( -\frac{1}{2\sigma\_i^2} \bigl\[ y\_i - \hat{y}\_i(\mathbf{a}) \bigr]^2 \Bigr)
$$

꼴이면, 그 음로그는

$$
\sum\_{i=1}^n \frac{1}{2\sigma\_i^2} \bigl\[ y\_i - \hat{y}\_i(\mathbf{a}) \bigr]^2
$$

와 같아서, 바로 가중 최소제곱의 목적함수와 동일해진다(사전분포 $p(\mathbf{a})$가 상수라고 가정할 경우). 정규화(regularization)는 사전분포를 가우시안 형태(예: $p(\mathbf{a}) = \exp(-\alpha |\mathbf{a}|^2)$)로 두는 베이지안 접근과 동일하며, 이는 릿지(Ridge) 또는 티코노프(Tikhonov) 정규화로 이어진다.

#### 미분동형(Robust) 오차함수와 IRLS

로버스트 추정을 위해 Huber, Lorentzian 등과 같은 오차함수를 사용하는 경우, 직접적으로는 비선형 최적화문제를 풀어야 한다. 하지만 IRLS(Iteratively Reweighted Least Squares) 알고리즘은 이 문제를 단계적으로 가중 최소제곱 형태로 풀 수 있게 만들어 준다.

예시로, Huber 오차함수를 정의하면, 잔차 $r\_i$의 크기가 작을 땐 제곱 오차처럼, 잔차가 큰 이상값 구간에선 절댓값 오차 형태로 동작한다. 이를 미분하면 각 잔차에 대해 $\psi(r\_i)$ 형태의 가중 계수를 얻게 되고, 그 계수를 통해 가중 최소제곱의 가중치가 업데이트된다. 이 때

$$
w\_i = \frac{\psi(r\_i)}{r\_i}
$$

와 비슷한 규칙이 반복적으로 적용되어 이상값의 영향력을 억제한다.

#### 오차모형 검증과 잔차 진단

가중치 적용 최소제곱 결과가 나왔을 때, 각 점의 잔차가 가정된 오차 분포(예: 독립 정규분포)에 어느 정도 부합하는지 확인해야 한다. 가중치가 정말로 분산의 역수가 맞는지, 오차 상관(Autocorrelation)이 없는지 등을 살펴보면서, 필요하면 추가적인 모형 수정이나 가중치 재설정이 이루어진다. 예를 들어, 잔차가 특정 구간에서 체계적으로 커지는 패턴이 보이면, 모델의 비선형성을 더 적절히 반영해야 하거나, 측정 과정에서 시스템적 오차가 발생했을 가능성이 있다.

오차분석에 있어서도, 가중치가 있다면 각 잔차의 기여도(가중 $w\_i$)가 달라지므로, 단순히 잔차의 크기만 보는 것이 아니라 $w\_i , r\_i^2$ 값 등을 통해 영향도를 평가한다. 가중화된 핫elling(Hotelling) 구, 혹은 Q-통계량 등을 사용해 극단값을 찾기도 한다.

#### 데이터 가공에서의 가중치 정책

실험이나 관측 데이터를 전처리할 때, 다양한 이유(신뢰도, 장비 상태, 환경 변화 등)로 인해 측정값이 서로 다른 품질을 지닐 수 있다. 어떤 데이터를 완전히 배제하기보다, "낮은 신뢰도"를 의미하는 작은 가중치를 부여하는 방식이 유연하다. 반대로, 중요한 구간이나 고정밀 계측기의 데이터에 대해서는 높은 가중치를 줄 수 있다. 이렇게 가중치를 부여하면, 전체 피팅 결과가 중요 구간에서 더 잘 맞도록 유도할 수도 있다.

가령 스펙트럼 분석에서 특정 파장대의 측정값이 노이즈가 심하다면, 그 구간에 대해 작은 가중치를 부여하여 전체 모델에 미치는 영향을 줄일 수 있다. 반대로, 특정 구간을 정밀하게 맞추고 싶다면(예: 화학 분석에서 특정 스펙트럼 피크 부근), 그 부근의 데이터를 더 크게 가중하여 모델 파라미터가 그 영역에 최적화되도록 할 수 있다.

#### 수치 안정성과 사후 검증

가중 최소제곱법에서 가장 많이 등장하는 연산은

$$
\mathbf{X}^T \mathbf{W} \mathbf{X}
$$

나

$$
\mathbf{J}(\mathbf{a}^{(k)})^T \mathbf{W} \mathbf{J}(\mathbf{a}^{(k)})
$$

의 역행렬(혹은 분해)을 구하는 과정이다. $\mathbf{W}$가 매우 큰 값이나 매우 작은 값을 지닌 대각원소로 구성되어 있으면, 수치적으로 ill-conditioned(컨디션 수가 매우 큰) 상태를 야기할 수 있다. 따라서 가중치의 크기가 극단적으로 달라지는 경우, 스케일링이나 정규화 기법을 통해 수치적 안정성을 개선해야 한다.

추가적으로, 해 $\mathbf{a}$가 구해진 뒤에 공분산 행렬

$$
(\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1}
$$

또는

$$
(\mathbf{J}^T \mathbf{W} \mathbf{J})^{-1}
$$

를 통해 파라미터 추정값의 신뢰구간을 검토한다. 가중치가 잘못 설정되어 있으면, 실제 오차 수준과 모델에서 추정된 오차 공분산이 크게 어긋나는 문제가 생길 수 있으므로, 이를 진단하여 재조정하는 과정을 거치기도 한다.

#### 추가 예시: C++로 구현한 선형 가중 최소제곱

선형 모델의 가중 최소제곱(Weighted Least Squares)을 C++ 코드로 구현해 볼 수 있다. 아래 예시는 다항회귀(직교 다항식을 쓰지 않은 일반적인 다항식) 문제를 가중하여 푸는 간단한 예시다. 관측값 $(x\_i, y\_i)$와 각 점별 분산 $\sigma\_i^2$가 주어졌다고 하자.

```c++
#include <iostream>
#include <vector>
#include <cmath>
#include <Eigen/Dense> // for matrix operations

// g++ -I /usr/include/eigen3 wls_demo.cpp -o wls_demo

int main() {
    // 예시 데이터 (x, y, sigma)
    std::vector<double> x_data = {0.1, 0.2, 0.4, 0.5, 0.7};
    std::vector<double> y_data = {1.05, 1.15, 1.48, 1.62, 2.04};
    std::vector<double> sigma_data = {0.01, 0.02, 0.02, 0.015, 0.03};
    
    // 다항식 차수
    int p = 2; // 예: 2차 다항식 a0 + a1 x + a2 x^2
    
    int n = x_data.size();
    // 디자인 행렬 X, 크기 n x (p+1)
    Eigen::MatrixXd X(n, p+1);
    Eigen::VectorXd y(n);
    Eigen::VectorXd w(n); // 가중치(= 1/sigma^2)
    
    for(int i = 0; i < n; i++){
        double xi = x_data[i];
        // [1, x, x^2, ...] 형태
        double power = 1.0;
        for(int j = 0; j <= p; j++){
            X(i, j) = power;
            power *= xi;
        }
        y(i) = y_data[i];
        // 분산이 sigma^2이므로 w = 1/sigma^2
        w(i) = 1.0 / (sigma_data[i]*sigma_data[i]);
    }
    
    // W 행렬은 diag(w), 크기 n x n
    Eigen::MatrixXd W = Eigen::MatrixXd::Zero(n, n);
    for(int i = 0; i < n; i++){
        W(i, i) = w(i);
    }
    
    // 정규방정식: (X^T W X) a = (X^T W y)
    Eigen::MatrixXd A = X.transpose() * W * X;
    Eigen::VectorXd B = X.transpose() * W * y;
    
    // 해 a = A^-1 B
    Eigen::VectorXd a = A.inverse() * B;
    
    // 결과 출력
    std::cout << "Estimated coefficients (a0, a1, a2): " << a.transpose() << std::endl;
    
    // 잔차 제곱합과 가중치 적용 오차 확인
    double S = 0.0;
    for(int i = 0; i < n; i++){
        double fit = 0.0;
        double power = 1.0;
        for(int j = 0; j <= p; j++){
            fit += a(j)*power;
            power *= x_data[i];
        }
        double r = y_data[i] - fit;
        S += w(i)*r*r; // 가중치 적용 잔차 제곱합
    }
    std::cout << "Weighted residual sum of squares S = " << S << std::endl;
    
    return 0;
}
```

위 코드는

1. `Eigen` 라이브러리를 사용하여 행렬 연산을 수행한다.
2. `(X^T W X)`와 `(X^T W y)`를 계산하여 해 $\mathbf{a}$를 구한다.
3. 결과로 얻은 계수 `a0, a1, a2`가, 가중 최소제곱법으로 적합된 2차 다항식의 파라미터다.
4. 잔차 $\mathbf{r}$에 대해 가중치가 적용된 제곱합 $S = \sum\_i w\_i,r\_i^2$를 확인할 수 있다.

#### 예측 구간 추정과 불확실성

가중 최소제곱으로 구한 파라미터 $\mathbf{a}$에 대해서, 특정 $x^\ast$에서의 예측값 $y^\ast = \mathbf{x}^\ast \mathbf{a}$(선형 모델의 경우)나 $f(x^\ast,\mathbf{a})$(비선형 모델의 경우)에 대한 오차 분산을 추정할 수도 있다. 예컨대 선형 다항 모델에서, $x^\ast$에서의 예측값 분산은

$$
\mathrm{Var}\bigl\[y^\ast\bigr] = \mathbf{x}^\ast (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} (\mathbf{x}^\ast)^T \sigma^2
$$

와 비슷한 형태를 갖는다(단, $\sigma^2$는 추정된 잔차분산). 일반화하면,

$$
\mathbf{V}\_{\mathbf{a}} = (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} \hat{\sigma}^2
$$

이므로, $\mathbf{x}^\ast (\mathbf{V}\_{\mathbf{a}}) (\mathbf{x}^\ast)^T$가 예측치의 분산을 제공한다. 비선형의 경우엔 야코비과 1차 근사를 이용한 유사한 식으로 추정 가능하다.

#### 부분 최소제곱(Partial Least Squares)과의 비교

데이터 차원이 매우 크고, 다중공선성(독립변수들이 서로 강하게 상관) 문제가 심각할 때, PLS(Partial Least Squares)나 PCR(Principal Component Regression)이 활용된다. 이 기법들은 디자인 행렬 $\mathbf{X}$의 열들 간 상관성을 줄이기 위해 차원 축소 및 특잉값분해 등의 방법을 결합하여 유효한 특성을 추출한다. 가중치가 필요한 경우, 각 관측점이나 변수를 스케일링할 때 가중치를 결합해 나가는 방식으로 확장 가능하다.

#### 정밀도 높은 수치 적분과 계수 추정

수치해석 관점에서, 가중 최소제곱법은 다른 수치 기법(예: 적분, 미분 방정식 해석 등)과 결합하여 “임의의 함수 근사”를 구하기도 한다. 예를 들어, 어떤 함수 $g(x)$를 개략적으로 알고 있고, 실험 데이터 $(x\_i, y\_i)$가 주어졌을 때, $y\_i \approx \int K(x\_i, t),\phi(t),dt$ 같은 형태를 최소제곱으로 피팅하려면, 적분 연산자 $K$에 대한 근사 및 가중 매트릭스를 동시에 고려해야 한다. 이 경우, 수치 적분의 오차를 반영한 가중치를 설계할 수도 있고, $\mathbf{X}$가 커널(핵) 함수로 구성된 “디자인 행렬”이 된다.

#### 심층학습(Deep Learning)과 가중 최소제곱

심층 신경망 기반의 모델들은 보통 확률적 경사하강법(SGD)으로 학습한다. 하지만 특정 계층(특히 마지막 레이어)을 최소제곱으로 푸는 방식과 결합할 수도 있다. 예: 오토인코더, 신경망 회귀 마지막 단계에서 오차분산을 가중치로 반영하여 “가중 Ridge 회귀”를 수행하면, 큰 잔차를 완화하면서 오버피팅을 방지할 수 있다. 이때 역전파(Backpropagation) 대신 닫힌 해를 찾거나, 적어도 부분적으로 해석적 해를 결합하여 학습 효율을 높이기도 한다.

\---은 말하지 않음

지금까지 다룬 가중 최소제곱법(가중 선형/비선형, 로버스트, 정규화, 대규모 문제 처리 등)은 실제 수치해석과 통계 데이터 분석의 여러 분야에서 핵심적 역할을 한다.