# QR 분해를 이용한 안정적 최소제곱해

일반적인 선형회귀나 곡선 근사 문제에서 최소제곱해를 구하기 위해 많이 사용하는 고전적 방법은 정상방정식에 의한 접근이다. 행렬을 $\mathbf{A}\in \mathbb{R}^{m\times n}$, 관측된 종속변수를 나타내는 열벡터를 $\mathbf{b}\in \mathbb{R}^{m}$이라 하고, 모수를 담은 미지의 열벡터를 $\mathbf{x}\in \mathbb{R}^{n}$라 할 때, 최소제곱문제는

$$
\min\_{\mathbf{x}\in \mathbb{R}^{n}} |\mathbf{A}\mathbf{x} - \mathbf{b}|\_2
$$

를 만족하는 $\mathbf{x}$를 찾는 것이다.

이 문제를 정상방정식으로 접근하면

$$
\mathbf{A}^\top \mathbf{A}\mathbf{x} = \mathbf{A}^\top \mathbf{b}
$$

를 풀어

$$
\mathbf{x} = (\mathbf{A}^\top \mathbf{A})^{-1}\mathbf{A}^\top \mathbf{b}
$$

를 얻게 된다. 그러나 이 방식은 $\mathbf{A}^\top \mathbf{A}$의 조건수가 좋지 않을 때 수치적으로 매우 불안정해질 수 있으며, 데이터의 수 $m$이 매우 큰 경우에도 효율적이지 않은 문제가 발생한다. 이러한 수치적 불안정을 완화하기 위해 QR 분해를 이용하는 방법이 자주 활용된다.

#### QR 분해의 개념

QR 분해(Orthogonal-Triangular Decomposition)는 실직선 공간에서 직교행렬(orthogonal matrix)과 상삼각행렬(upper triangular matrix)의 곱으로 주어진 행렬을 분해하는 기법이다. 즉 $\mathbf{A}\in \mathbb{R}^{m\times n}$, $m\ge n$인 경우,

$$
\mathbf{A} = \mathbf{Q}\mathbf{R}
$$

로 나타낼 수 있다. 여기에서 $\mathbf{Q}\in \mathbb{R}^{m\times m}$은 직교행렬, $\mathbf{R}\in \mathbb{R}^{m\times n}$은 위가 삼각 형태를 띠고 나머지는 0인 행렬이다. 그러나 실질적으로는 $\mathbf{A}$의 열 개수를 $n$이라 하면, 앞의 $m\times n$ 형태로 $\mathbf{R}$를 자른 뒤 이를 $\mathbf{R}\_1\in \mathbb{R}^{n\times n}$이라 하고, 나머지를 $\mathbf{R}\_2$라고 하면

$$
\mathbf{A} = \mathbf{Q} \begin{bmatrix} \mathbf{R}\_1 \ \mathbf{0} \end{bmatrix}
$$

와 같은 구조로 나타낸다. 즉 실제 계산에서는 직교행렬도 $\mathbf{Q}\in \mathbb{R}^{m\times n}$의 일부 열만 활용하거나, 하우스홀더 반사(Householder reflection) 등을 통해 필요한 만큼만 연산하여 메모리를 아끼는 방법을 쓴다.

일단 $\mathbf{Q}$가 직교행렬이라는 말은

$$
\mathbf{Q}^\top \mathbf{Q} = \mathbf{I}
$$

를 만족함을 의미하고, 이는 $\mathbf{Q}$가 열벡터들로 구성된 정규직교 기저(orthonormal basis)임을 말해준다.

#### QR 분해를 활용한 최소제곱문제 해

$\mathbf{A}\mathbf{x}\approx \mathbf{b}$의 최소제곱해를 구하기 위해 QR 분해를 사용하면 다음의 과정을 거치게 된다. 먼저

$$
\mathbf{A} = \mathbf{Q}\mathbf{R}
$$

로 분해한다. 선형계 $\mathbf{A}\mathbf{x} = \mathbf{b}$는

$$
\mathbf{Q}\mathbf{R}\mathbf{x} = \mathbf{b}
$$

와 같으며, 양변에 $\mathbf{Q}^\top$를 곱해주면

$$
\mathbf{R}\mathbf{x} = \mathbf{Q}^\top \mathbf{b}
$$

가 된다. 여기서 $\mathbf{R}\in \mathbb{R}^{m\times n}$는 사실상 위쪽 $n\times n$ 부분만이 역행렬 계산에 관여하므로, 이를 $\mathbf{R}\_1$이라 하면

$$
\mathbf{R}\_1\mathbf{x} = \mathbf{Q}\_1^\top \mathbf{b}
$$

를 통해 $\mathbf{x}$를 바로 결정할 수 있다. 이때 $\mathbf{Q}\_1$은 $\mathbf{Q}$의 앞쪽 열 혹은 차원을 축소한 형태라고 볼 수 있으며, $\mathbf{R}\_1$은 정방행렬이고 상삼각 형태이므로 뒤에서부터 간단히 되풀이대입(back-substitution)으로 $\mathbf{x}$를 구할 수 있다.

위 과정을 통해 알 수 있듯이, $\mathbf{A}^\top \mathbf{A}$와 같은 연산이 직접적으로 사용되지 않아 $\mathbf{A}^\top \mathbf{A}$의 조건수가 매우 나빠지는 경우에도 수치안정성이 훨씬 향상된다. 실제 구현에서는 클래식 그람-슈미트(Classical Gram-Schmidt)나 수정 그람-슈미트(Modified Gram-Schmidt)보다는 하우스홀더 반사를 통한 QR 분해가 수치안정성이 우수하다고 알려져 있다. 하우스홀더 반사를 이용하면 직교행렬을 명시적으로 저장하지 않고도 연산을 효율적으로 수행할 수 있다.

#### Householder를 이용한 QR 분해의 개요

하우스홀더 반사는 벡터를 특정 목표방향에 대해 대칭 반사(reflection)시키는 선형변환이다. 어떤 단위벡터 $\mathbf{u}$가 주어지면, 하우스홀더 행렬 $\mathbf{H}$는

$$
\mathbf{H} = \mathbf{I} - 2, \mathbf{u}\mathbf{u}^\top
$$

의 형태로 정의된다. $\mathbf{u}$는 정규화되었으므로 $|\mathbf{u}|\_2 = 1$이며, $\mathbf{H}$는 직교행렬임을 확인할 수 있다.

특정 열벡터 $\mathbf{v}$를 $\mathbf{e}\_1$ 방향으로 보내려 할 때, $\mathbf{u}$를 적절히 잡으면 $\mathbf{v}$를 $|\mathbf{v}|\_2 \mathbf{e}\_1$로 바꿀 수 있다. 이 과정을 반복하여 $\mathbf{A}$의 각 열을 삼각 형태로 만들어가면 상삼각 형태의 $\mathbf{R}$가 완성되며, 그 과정을 추적해보면 곧 $\mathbf{Q} = \mathbf{H}\_1 \mathbf{H}\_2 \dots \mathbf{H}\_n$으로 주어진다는 사실을 알게 된다.

이러한 방법이 안정적인 이유는, 각 단계에서 소멸시키려는 원소들을 대규모의 곱셈 없이 직교변환으로 처리하기 때문이다. 그람-슈미트는 내적 연산에 의존하기에, 벡터 크기가 매우 커지거나 작아질 때 발생하는 부동소수점 연산 오류에 민감하지만, 하우스홀더 방식은 상대적으로 수치오차를 더 잘 억제한다.

#### 최소제곱문제에 대한 QR 분해의 활용

최소제곱문제 $|\mathbf{A}\mathbf{x}-\mathbf{b}|\_2$를 풀 때, 하우스홀더 QR 분해를 구하면 결국

$$
\mathbf{Q}^\top \mathbf{A} =  \begin{bmatrix} \mathbf{R}\_1 \ \mathbf{0} \end{bmatrix}
$$

와

$$
\mathbf{Q}^\top \mathbf{b} =  \begin{bmatrix} \mathbf{c}\_1 \ \mathbf{c}\_2 \end{bmatrix}
$$

를 얻게 된다. 여기서 $\mathbf{R}\_1\in \mathbb{R}^{n\times n}$은 가역행렬이고 $\mathbf{c}\_1\in \mathbb{R}^{n}$, $\mathbf{c}\_2\in \mathbb{R}^{m-n}$이다. 최소제곱해 $\mathbf{x}$는

$$
\begin{align} \mathbf{R}\_1 \mathbf{x}  &= \mathbf{c}\_1 \ \mathbf{x}  &= \mathbf{R}\_1^{-1} \mathbf{c}\_1 \end{align}
$$

로 얻어지며, 이때 $\mathbf{c}\_2$는 $\mathbf{A}\mathbf{x}$로도 설명할 수 없는 잔차(residual)에 해당한다. 이 방식은 $\mathbf{A}^\top \mathbf{A}$를 전혀 구성하지 않으므로 수치적으로 훨씬 더 안정적이다.

#### Python 예시

아래는 간단히 Python에서 Numpy를 사용하여 QR 분해를 통해 최소제곱해를 구하는 예시 코드이다.

```python
import numpy as np

# 예제 데이터 구성
# m >= n 인 상황에서 최소제곱문제
m = 6
n = 3
np.random.seed(0)
A = np.random.randn(m, n)
b = np.random.randn(m)

# QR 분해를 수행 (numpy.linalg.qr는 기본적으로 경제적 분해 모드='reduced' 사용)
Q, R = np.linalg.qr(A, mode='reduced')

# Q^T b를 계산
b_tilde = Q.T @ b

# 상삼각 행렬 R을 이용하여 뒤에서부터 x 구하기
x = np.linalg.solve(R, b_tilde)

print("A:", A)
print("b:", b)
print("Q:", Q)
print("R:", R)
print("x (Least Squares Solution):", x)
```

이 코드에서는 $A$를 QR 분해한 후, $Q^\top b$를 계산하고 나서 $R x = Q^\top b$를 푸는 과정을 수행한다. 실제 대규모 데이터에서는 하우스홀더 QR 분해 알고리즘이나, 더 효율적이고 안정적인 방식(예: LAPACK 루틴)이 내부적으로 수행된다. 이 과정을 통해 최소제곱해가 구해진다.

#### Givens 회전을 통한 QR 분해

하우스홀더 반사 외에도 QR 분해를 구현하는 또 다른 대표적 방법으로 Givens 회전(Givens rotation)이 있다. Givens 회전은 어떤 2차원 평면 상에서의 회전을 이용해 특정 원소를 소멸시키는 방식으로, 행렬의 한 원소를 0으로 만들기 위해 필요한 회전각을 구한 뒤 회전에 해당하는 직교행렬을 왼쪽에서 곱하는 과정을 반복한다. 하우스홀더 반사가 하나의 열을 한 번에 처리하는 반면, Givens 회전은 각각의 원소를 개별적으로 제거하므로, 대규모 희소(sparse) 행렬일 때 유리한 경우가 있다.

어떤 행렬 $\mathbf{A}\in \mathbb{R}^{m\times n}$에 대해, $i$번째 행과 $j$번째 행에서 한 원소를 제거하려 할 때, Givens 회전 행렬 $\mathbf{G}\in \mathbb{R}^{m\times m}$는 항등행렬 $\mathbf{I}\in \mathbb{R}^{m\times m}$의 부분에 다음과 같은 2차원 회전행렬을 끼워 넣은 형태로 구성된다. $\mathbf{G}$는 1부터 $m$까지의 행과 열 인덱스 중 $(i,i)$, $(j,j)$, $(i,j)$, $(j,i)$ 위치에만 회전 요소가 들어가고, 나머지는 모두 항등 요소를 유지한다. 이를 수식으로 나타내면,

$$
\mathbf{G}(i,i) = c,\quad \mathbf{G}(i,j) = s,\quad \mathbf{G}(j,i) = -s,\quad \mathbf{G}(j,j) = c
$$

가 되는데, $c$와 $s$는 $\mathbf{A}$의 해당 원소를 0으로 만들기 위해 적절히 설정된 코사인, 사인에 해당한다. 실제 계산에서 $c$와 $s$는

$$
c = \frac{a}{\sqrt{a^2 + b^2}},\quad s = \frac{b}{\sqrt{a^2 + b^2}}
$$

와 같은 방식으로 구한다. 여기서 $a, b$는 목표로 하는 두 행($i$행과 $j$행)에 있는 원소들이다. 이렇게 구성된 $\mathbf{G}$를 $\mathbf{A}$에 왼쪽에서 곱해주면, 특정 위치의 원소가 0으로 만들어지며, 이것을 반복해 나가면 결국 상삼각 행렬 $\mathbf{R}$가 형성된다. 모든 Givens 회전 행렬의 곱을 역순으로 취합하면 $\mathbf{Q}$를 얻을 수 있다.

하우스홀더 반사와 달리, Givens 회전 방식은 다량의 0 원소가 존재하는 희소 행렬이나 실시간으로 행렬의 특정 위치만 수정해야 하는 상황에서 유리한 편이다. 예를 들어 선형회귀 분석을 실시간으로 업데이트해야 하는 경우, 하우스홀더 방식을 매번 전면 재계산하기보다 Givens 회전을 사용하면 국소적 업데이트가 가능하다. 반면, 전체 행렬이 빽빽한(dense) 경우에는 하우스홀더가 더 효율적인 연산량을 보일 수 있다.

#### 컬럼 피벗팅을 통한 순위부족 문제 대처

실제 데이터나 응용 상황에서 $\mathbf{A}$가 완전한 열 순위(full column rank)를 가지지 않을 수도 있다. 즉 $\mathbf{A}$의 열들 중 선형종속이 존재하거나, 매우 높은 상관관계로 인해 열벡터들이 거의 선형의존 상태에 가까울 때가 있다. 이런 상황에서는 최소제곱문제를 풀 때 난해하거나 수치적으로 불안정해질 수 있다. 이를 방지하기 위해 QR 분해를 수행할 때 컬럼 피벗팅(column pivoting)을 적용하기도 한다.

컬럼 피벗팅은 $\mathbf{A}=\mathbf{Q}\mathbf{R}$ 대신

$$
\mathbf{A}\mathbf{P} = \mathbf{Q}\mathbf{R}
$$

와 같은 형태로 $\mathbf{P}$라는 열 교환 행렬(permutation matrix)을 도입하여 수행한다. $\mathbf{P}$는 항등행렬에서 열 순서를 교환한 형태이므로 직교행렬(정확히는 정방교환행렬)이며, 결국 열들의 순서를 재배치한 행렬 $\mathbf{A}\mathbf{P}$에 대해 QR 분해를 진행하는 것과 같다. 그 결과,

$$
\mathbf{R} =  \begin{bmatrix} \mathbf{R}*{11} & \mathbf{R}*{12} \ \mathbf{0} & \mathbf{R}\_{22} \end{bmatrix}
$$

형태로 부블록이 뚜렷해지면, $\mathbf{R}\_{22}$가 근사적으로 0에 가까워지는 지점을 감지해 $\mathrm{rank}(\mathbf{A})$를 추정할 수 있다. 다시 말해, $\mathbf{A}$가 열 순위가 낮거나, 열들 간의 상관성이 높은 상황에서도, 컬럼 피벗팅을 통해 열 순서를 적절히 선택하여 수치 안정도를 상당히 개선할 수 있다.

#### 경제적 QR(Economy-size QR) 분해

$\mathbf{A}\in \mathbb{R}^{m\times n}$에서 $m\ge n$이라고 할 때, 완전한 형태의 $\mathbf{Q}\in \mathbb{R}^{m\times m}$을 전부 저장하는 것은 비효율적일 수 있다. 최소제곱문제에서 실제로 필요한 것은 $\mathbf{Q}$의 앞쪽 $n$개 열만이므로, 그것을 $\mathbf{Q}\_1\in \mathbb{R}^{m\times n}$이라 하고, $\mathbf{R}$ 또한 앞쪽 $n\times n$ 부분을 $\mathbf{R}\_1$이라 한다. 그 결과

$$
\mathbf{A} = \mathbf{Q}\_1 \mathbf{R}\_1
$$

와 같이 쓸 수 있고, 이 경우 $\mathbf{Q}\_1$는 직교행렬의 부분 형태로서 $\mathbf{Q}\_1^\top \mathbf{Q}*1 = \mathbf{I}*{n\times n}$를 만족한다. 이를 경제적 QR 분해(Economy-size QR decomposition)라고 부르며, 실제 구현(예: NumPy, LAPACK 등)에서는 이 경제적 QR 분해가 주로 사용된다. 최소제곱문제를 풀 때 $\mathbf{Q}\_1$와 $\mathbf{R}\_1$만 있으면 충분하므로 연산량과 메모리를 절감할 수 있다.

#### C++ 예시

C++ 언어에서 Eigen 라이브러리를 사용하면 QR 분해를 매우 간단하게 구현할 수 있다. 아래는 간략한 예시이다.

```cpp
#include <iostream>
#include <Eigen/Dense>

int main() {
    using namespace Eigen;
    int m = 6;
    int n = 3;
    
    // 예제 행렬 A와 벡터 b 생성
    MatrixXd A = MatrixXd::Random(m, n);
    VectorXd b = VectorXd::Random(m);

    // Householder QR 분해
    HouseholderQR<MatrixXd> qr(A);
    MatrixXd Q = qr.householderQ(); 
    MatrixXd R = qr.matrixQR().triangularView<Upper>();

    // 최소제곱해 x 계산
    // 경제적 QR 사용 시에는 R과 Q의 일부분만 이용해도 무방
    VectorXd x = qr.solve(b);

    std::cout << "A:\n" << A << std::endl;
    std::cout << "b:\n" << b << std::endl;
    std::cout << "Q:\n" << Q << std::endl;
    std::cout << "R:\n" << R << std::endl;
    std::cout << "x (Least Squares Solution):\n" << x << std::endl;

    return 0;
}
```

이 예시는 Eigen 라이브러리의 HouseholderQR 클래스를 이용하여 QR 분해를 계산한 뒤, `.solve(b)`를 통해 최소제곱해를 구한다. 내부적으로는 $R x = Q^\top b$ 과정을 수행하며, $Q$와 $R$를 직접 확인할 수도 있다. 실무 환경에서는 굳이 $Q$와 $R$를 명시적으로 모두 저장하지 않고, 필요한 연산만을 수행하여 해를 효율적으로 구하기도 한다.

#### 수치 안정성 분석

QR 분해를 통한 최소제곱해 구하기는 $\mathbf{A}^\top \mathbf{A}$를 직접 다루지 않으므로, 정상방정식 접근에 비해 훨씬 양호한 수치 안정성을 제공한다. 이 점을 좀 더 구체적으로 살펴보면, 우선 $\mathbf{Q}$가 직교행렬이라는 사실이 매우 중요하다. 직교행렬은 다음과 같은 성질이 있다:

$$
|\mathbf{Q}\mathbf{x}|\_2 = |\mathbf{x}|\_2
$$

직교행렬의 곱셈은 벡터의 길이를 유지하기 때문에, 불필요한 스케일링으로 인한 수치적 오차 증폭이 최소화된다. 반면에 정상방정식을 풀 때는 $\mathbf{A}^\top \mathbf{A}$가 포함되는데, 이는 $\mathbf{A}$의 조건수 $\kappa(\mathbf{A})$를 제곱으로 악화시키는 효과가 있다:

$$
\kappa(\mathbf{A}^\top \mathbf{A}) = \[\kappa(\mathbf{A})]^2
$$

조건수(condition number)가 커질수록 미세한 오차가 해 전체에 크게 반영될 가능성이 높아지므로, $\mathbf{A}^\top \mathbf{A}$가 직접 등장하는 방법보다는 QR, SVD 같은 직교(혹은 유니터리) 변환을 이용한 방식이 일반적으로 더 안정적이다.

하우스홀더 반사를 이용한 방식은 이론적으로도 매우 튼튼한 안정성을 가진다고 알려져 있다. Givens 회전 역시 안정적이지만, 각 원소를 개별적으로 소거해 나가기 때문에 전체 연산 횟수 면에서 하우스홀더보다 다소 비효율적일 수 있다. 그러나 희소 행렬이나 실시간 업데이트 상황에서는 Givens 회전이 유리할 수 있다.

#### 최소제곱해와 잔차의 관찰

$\mathbf{x}$가 최소제곱해일 때, 즉

$$
\mathbf{x} = \underset{\mathbf{x}\in \mathbb{R}^n}{\mathrm{argmin}}|\mathbf{A}\mathbf{x}-\mathbf{b}|\_2
$$

를 만족하는 $\mathbf{x}$를 찾았다고 하면, 그 때의 잔차(residual) 벡터 $\mathbf{r}$는

$$
\mathbf{r} = \mathbf{b} - \mathbf{A}\mathbf{x}
$$

가 된다. QR 분해로 해를 구하면

$$
\mathbf{Q}^\top \mathbf{A} =  \begin{bmatrix} \mathbf{R}\_1 \ \mathbf{0} \end{bmatrix}, \quad \mathbf{Q}^\top \mathbf{b} =  \begin{bmatrix} \mathbf{c}\_1 \ \mathbf{c}\_2 \end{bmatrix}
$$

에 따라

$$
\mathbf{R}\_1 \mathbf{x} = \mathbf{c}\_1
$$

를 만족하므로,

$$
\mathbf{x} = \mathbf{R}\_1^{-1} \mathbf{c}\_1
$$

로 해를 구할 수 있다. 이때

$$
\mathbf{c}\_2 = \mathbf{Q}\_2^\top \mathbf{b}
$$

(여기서 $\mathbf{Q} = \begin{bmatrix}\mathbf{Q}\_1 & \mathbf{Q}\_2\end{bmatrix}$)는 실제로 $\mathbf{A}\mathbf{x}$로는 설명되지 않는 오차 부분이다. 따라서

$$
\mathbf{r} = \mathbf{b} - \mathbf{A}\mathbf{x}  = \mathbf{b} - \mathbf{Q}\mathbf{R}\mathbf{x} = \mathbf{b} - \mathbf{Q} \begin{bmatrix}\mathbf{R}\_1 \ 0\end{bmatrix} \mathbf{x} = \mathbf{b} - \mathbf{Q} \begin{bmatrix}\mathbf{R}\_1\mathbf{x} \ \mathbf{0}\end{bmatrix} = \mathbf{b} - \mathbf{Q} \begin{bmatrix}\mathbf{c}\_1 \ \mathbf{0}\end{bmatrix}
$$

이고, 이를 $\mathbf{Q}$의 직교성으로 풀면

$$
\mathbf{r} = \mathbf{Q} \Big( \mathbf{Q}^\top \mathbf{b} -  \begin{bmatrix}\mathbf{c}\_1 \ \mathbf{0}\end{bmatrix} \Big) = \mathbf{Q} \begin{bmatrix}\mathbf{c}\_1 - \mathbf{c}\_1 \ \mathbf{c}\_2\end{bmatrix} = \mathbf{Q} \begin{bmatrix}\mathbf{0}\ \mathbf{c}\_2\end{bmatrix} = \mathbf{Q}\_2 \mathbf{c}\_2
$$

가 된다. 즉 잔차 벡터는 $\mathbf{Q}\_2$ (정규직교 기저의 나머지 부분) 방향으로만 존재한다. 따라서 QR 분해를 통해 구해진 잔차의 구조와 크기를 이해하고, 필요한 경우 $\mathbf{c}\_2$를 분석함으로써 데이터가 얼마나 잘 근사되고 있는지를 확인할 수 있다.

#### QR 분해와 SVD(특이값 분해) 비교

더 나아가, 최소제곱해 문제에서 가장 안정적인 방법 중 하나로 꼽히는 것은 SVD(특이값 분해, Singular Value Decomposition) 접근이다. SVD는 $\mathbf{A}$를

$$
\mathbf{A} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^\top
$$

로 분해하되, $\mathbf{U}$와 $\mathbf{V}$가 각각 직교행렬이고, $\boldsymbol{\Sigma}$는 대각선에 특이값(singular value)이 정렬된 형태이다. 이때 최소제곱해는

$$
\mathbf{x} = \mathbf{V},\boldsymbol{\Sigma}^{\dagger},\mathbf{U}^\top,\mathbf{b}
$$

로 구할 수 있다. 여기서 $\boldsymbol{\Sigma}^{\dagger}$는 $\boldsymbol{\Sigma}$의 무어-펜로즈 역행렬(Moore-Penrose inverse)이며, 특이값 중 매우 작은 값들은 안정성 측면에서 제거(또는 감쇠)할 수도 있어, 랭크 결핍(rank-deficient) 상황이나 노이즈가 큰 상황에서의 해석에도 유리하다.

SVD는 이론적으로는 수치 안정성 면에서 가장 탁월한 방법이지만, 일반적으로 QR 분해에 비해 계산비용이 크다. 선형회귀나 대규모 데이터의 경우에는 QR 분해가 더 빠르고 실용적이므로, 보통은 QR 분해를 쓰되 수치 문제가 심각하거나, 정확도가 매우 요구되는 상황에서는 SVD를 사용하기도 한다.

#### 투영 행렬의 해석

최소제곱해 $\mathbf{x}$를 통해 $\mathbf{A}\mathbf{x}$는 $\mathbf{A}$의 열공간(column space)에 속하는 벡터 중 $\mathbf{b}$와 가장 가까운 벡터가 된다. 즉 정사영(projection)을 생각하면,

$$
\mathbf{\hat{b}} = \mathbf{A}\mathbf{x}
$$

는 $\mathbf{b}$의 $\mathbf{A}$의 열공간으로의 직교정사영에 해당한다. 이를 투영행렬(projection matrix) $\mathbf{P}$라 할 때,

$$
\mathbf{\hat{b}} = \mathbf{P}\mathbf{b}
$$

로 표현할 수 있다. 정상방정식 관점에서 보면,

$$
\mathbf{P} = \mathbf{A}(\mathbf{A}^\top \mathbf{A})^{-1}\mathbf{A}^\top
$$

형태이지만, QR 분해 관점에서는

$$
\mathbf{P} = \mathbf{Q}\_1 \mathbf{Q}\_1^\top
$$

가 된다. 여기서 $\mathbf{Q}\_1\in \mathbb{R}^{m\times n}$은 열벡터들이 정규직교를 이루므로, $\mathbf{Q}\_1 \mathbf{Q}\_1^\top$는 차원 $m\times m$의 투영행렬이 된다. 투영행렬의 대각원소 합은 투영되는 공간의 차원과 같으므로, 여기서는 $n$에 해당한다.

#### 가중 최소제곱(Weighted Least Squares)의 QR 분해

일부 응용에서는 모든 관측값에 동일한 중요도를 두는 것이 합리적이지 않을 수 있다. 예를 들어 특정 구간에서 측정값의 신뢰도가 낮거나 오차가 큰 반면, 다른 구간에서는 상대적으로 정밀한 관측값을 얻을 수 있는 경우가 있다. 이런 상황에서 보다 정확한 근사를 위해, 각 데이터 점에 대해 서로 다른 가중치(weight)를 부여하는 가중 최소제곱문제(Weighted Least Squares)를 고려한다.

가중행렬(대각 행렬) $\mathbf{W}\in \mathbb{R}^{m\times m}$가 주어졌다고 하자. 예컨대 $\mathbf{W}$가

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

형태인 대각행렬이라고 할 때, 우리의 목표는 다음을 최소화하는 것이다.

$$
\min\_{\mathbf{x}\in \mathbb{R}^{n}} | \mathbf{W}^{1/2}(\mathbf{A}\mathbf{x} - \mathbf{b}) |\_2
$$

여기서 $\mathbf{W}^{1/2}$는 $\mathbf{W}$의 대각원소에 각각 루트를 취한 것으로,

$$
\mathbf{W}^{1/2} =  \begin{bmatrix} \sqrt{w\_1} & 0 & \cdots & 0 \ 0 & \sqrt{w\_2} & \cdots & 0 \ \vdots & \vdots & \ddots & \vdots \ 0 & 0 & \cdots & \sqrt{w\_m} \end{bmatrix}
$$

이다. 이렇게 정의하면,

$$
\mathbf{W}^{1/2}(\mathbf{A}\mathbf{x} - \mathbf{b}) |\_2^2  = (\mathbf{A}\mathbf{x} - \mathbf{b})^\top \mathbf{W} (\mathbf{A}\mathbf{x} - \mathbf{b})
$$

가 된다.

이제 $\mathbf{\tilde{A}} \equiv \mathbf{W}^{1/2}\mathbf{A}$, $\mathbf{\tilde{b}} \equiv \mathbf{W}^{1/2}\mathbf{b}$ 라고 놓으면, 문제는

$$
\min\_{\mathbf{x}\in \mathbb{R}^{n}} |\mathbf{\tilde{A}}\mathbf{x} - \mathbf{\tilde{b}}|\_2
$$

가 된다. 이는 결국 일반적인 최소제곱문제와 동일한 형태이므로, QR 분해를 이용하여 안정적으로 해를 구할 수 있다. 즉 다음과 같은 단계를 거친다.

먼저, $\mathbf{\tilde{A}}$에 대해 QR 분해를 수행하면

$$
\mathbf{\tilde{A}} = \mathbf{Q}\mathbf{R}
$$

가 된다. 여기서 $\mathbf{Q}$는 $m\times m$ 직교행렬(또는 필요에 따라 경제적 분해 시 $m\times n$)이고, $\mathbf{R}$는 $m\times n$의 상삼각 형태다(상단 $n\times n$ 부분이 유효). 양변에 $\mathbf{Q}^\top$를 곱해주면

$$
\mathbf{Q}^\top \mathbf{\tilde{A}}, \mathbf{x}  = \mathbf{Q}^\top \mathbf{\tilde{b}},
$$

즉

$$
\mathbf{R}\mathbf{x} = \mathbf{Q}^\top \mathbf{\tilde{b}}
$$

가 되고, $\mathbf{R}$의 상단 정방 부분을 $\mathbf{R}\_1$라 할 때,

$$
\mathbf{R}\_1 \mathbf{x} = \mathbf{Q}\_1^\top \mathbf{\tilde{b}}
$$

가 성립한다. 따라서

$$
\mathbf{x} = \mathbf{R}\_1^{-1} \mathbf{Q}\_1^\top \mathbf{\tilde{b}}
$$

를 통해 가중 최소제곱해를 얻는다. 여기서 $\mathbf{\tilde{b}} = \mathbf{W}^{1/2}\mathbf{b}$임을 기억하면, 실제 해석에선

$$
\mathbf{x} = \mathbf{R}\_1^{-1} \mathbf{Q}\_1^\top (\mathbf{W}^{1/2},\mathbf{b})
$$

가 된다.

가중값 $\mathbf{W}$의 대각원소들이 클수록 해당 관측값의 오차를 더 크게 벌(penalize) 주는 역할을 하기 때문에, 그 관측값에 대한 근사를 더 정확히 하도록 해를 유도한다. 반대로 작은 가중이 부여된 관측값은 오차가 다소 커져도 전체 목적함수에 큰 영향을 주지 않는다. 수치해석적으로는 마찬가지로 QR 분해를 통해 안정적으로 해를 구할 수 있으므로, 원 자료에 일괄 변환을 가해 $\mathbf{\tilde{A}}$, $\mathbf{\tilde{b}}$를 형성한 뒤, 일반적인 방식대로 최소제곱해를 구하면 된다.

#### 랭크 결핍(Deficiency)이 있는 경우의 QR 접근

실제 데이터는 종종 랭크가 결핍된 상태, 즉 $\mathrm{rank}(\mathbf{A}) < n$일 수 있다. 이런 경우 $\mathbf{A}$의 열들이 선형종속적이라서 일반적인 최소제곱문제 해 $\mathbf{x}$가 유일하지 않다. 또한 $\mathbf{A}^\top \mathbf{A}$가 가역이 아닐 수 있다(또는 수치적으로 사실상 단수(singular)에 가깝다). QR 분해를 이용하면 컬럼 피벗팅(column pivoting)을 통해 유효랭크를 판별해볼 수 있다.

컬럼 피벗팅을 동반한 QR 분해는

$$
\mathbf{A}\mathbf{P} = \mathbf{Q}\mathbf{R}
$$

로 나타나며, $\mathbf{P}$는 교환행렬(permutation matrix)이다. $\mathbf{R}$를 상위 블록과 하위 블록으로 나누면,

$$
\mathbf{R} =  \begin{bmatrix} \mathbf{R}*{11} & \mathbf{R}*{12} \ \mathbf{0} & \mathbf{R}\_{22} \end{bmatrix},
$$

여기서 $\mathbf{R}*{11}$의 크기는 $r\times r$, $r=\mathrm{rank}(\mathbf{A})$라고 할 때, $\mathbf{R}*{22}$는 거의 0에 가까운 값이 있는 블록으로 분리될 수 있다. 그렇게 되면 사실상

$$
\mathbf{Q}^\top \mathbf{b} =  \begin{bmatrix} \mathbf{b}\_1 \ \mathbf{b}\_2 \end{bmatrix}
$$

에 대해, 상단부

$$
\mathbf{R}\_{11},\mathbf{x}*1 + \mathbf{R}*{12},\mathbf{x}\_2 = \mathbf{b}\_1
$$

만이 유효하게 동작하고, 하단부는

$$
\mathbf{0},\mathbf{x}*1 + \mathbf{R}*{22},\mathbf{x}\_2 \approx \mathbf{b}\_2
$$

에서 $\mathbf{R}\_{22}$가 매우 작으므로, $\mathbf{x}*2$에 대한 구속이 사실상 느슨해지거나 무효가 된다(즉 무수히 많은 해가 존재). 따라서 랭크 결핍 상태임을 식별하고, $\mathbf{R}*{22}$가 충분히 작으면 해당 부분을 0으로 처리하거나, SVD 접근으로 특이값이 매우 작은 방향을 제거함으로써 해를 특정 형태로 선택할 수도 있다.

#### 해석적 예시 (Octave)

아래는 간단한 Octave 스크립트 예시로, 가중 최소제곱과 컬럼 피벗팅 QR을 함께 시연해볼 수 있다.

```octave
% Octave 예시
% 가중 최소제곱 + 컬럼 피벗팅 QR

% 랜덤 시드 고정
rand("seed", 0);

% 간단한 데이터 생성
m = 6;
n = 4;
A = randn(m, n);
b = randn(m, 1);

% 특정 열을 인위적으로 거의 선형 결합상태로 만들기
A(:,4) = A(:,1) + 0.0001 * A(:,2);  

% 가중 벡터를 임의로 설정
Wdiag = [10; 1; 5; 1; 1; 2];
W = diag(Wdiag);

% 가중 최소제곱문제에서의 A, b 변환
Atilde = sqrt(W) * A;
btilde = sqrt(W) * b;

% 컬럼 피벗팅 QR 분해
[Q, R, P] = qr(Atilde, 0);   % economy-size QR with pivot
% Q: m x n, R: n x n, P: n x n

% 최소제곱해 x 계산
btilde_tilde = Q' * btilde;    % Q' btilde
x_tilde = R \ btilde_tilde;    % R x = btilde_tilde
x_wls = P * x_tilde;           % 원래 열 순서 복원

disp("A = "), disp(A);
disp("b = "), disp(b);
disp("W = "), disp(W);
disp("Pivot P = "), disp(P);
disp("R = "), disp(R);
disp("x_wls (Weighted LS) = "), disp(x_wls);
```

이 예시에서

1. `A`와 `b`에 가중 `W`를 반영하여 `Atilde`, `btilde`를 만든다.
2. `[Q,R,P]=qr(Atilde,0)` 명령을 통해 컬럼 피벗팅 QR 분해를 수행한다. (Octave에서 두 번째 인자로 0을 주면 경제적 QR 분해를 의미한다.)
3. 해는 $R \backslash (Q^\top \texttt{btilde})$로 구한 뒤, `P`를 통해 원래 열 순서를 되돌려놓아야 한다. (즉 $x = P \cdot \texttt{x\_tilde}$)

이 과정을 거치면, 만약 열 간 선형종속성이 심각하다면 $R$의 특정 부분이 거의 0으로 나타나 수치 안정도가 위험 신호를 내보낼 것이다. 랭크가 결핍된 경우에 대해선, 실제로는 SVD를 함께 사용하여 작은 특이값들을 제거하거나, 정칙화(regularization)를 도입하기도 한다.

#### 고차원 문제 및 대규모 데이터에서의 QR

실제 응용에서는 $m$과 $n$이 수천, 수만, 혹은 그 이상으로 큰 경우가 많다. 이때 전통적인 QR 분해(하우스홀더 또는 Givens)도 연산량이 매우 커질 수 있다. 대규모 문제에서는 다음과 같은 방안이 고려된다.

* 블록 QR(Block QR): 행렬을 일부분씩 나누어 QR 분해를 수행하고, 결과를 합치는 기법.
* 반복적 방법(Iterative methods): GMRES(Generalized Minimal Residual), LSQR(Least Squares QR) 같은 반복 알고리즘을 통해 축소된 시스템을 순차적으로 근사하면서 해를 갱신.
* 희소행렬을 이용한 알고리즘: Givens나 하우스홀더 반사를 희소 구조에 맞춰 최적화.
* 랜덤화 기법(Randomized algorithms): 입력 행렬을 무작위로 스케치(sketch)하여 저차원 근사를 구한 뒤, 효율적으로 해를 추정.

이와 같은 방법들을 QR 관점에서 해석하면, 결국 직교 투영의 개념이 중심에 있다는 사실은 변하지 않는다. 다만, 연산량과 메모리 사용량을 줄이거나, 랭크 결핍 및 정칙화 이슈를 극복하기 위해 다양한 변형 기법이 함께 연구되고 있다.

#### Rank-Revealing QR(RRQR)과 랭크 추정

컬럼 피벗팅을 동반한 QR 분해는 행렬의 랭크 추정에 어느 정도 도움을 준다. 그러나 좀 더 효과적으로 “랭크를 밝혀내는” 분해 기법으로 Rank-Revealing QR(RRQR) 알고리즘이 제안되었다. RRQR은 $\mathbf{A}\in \mathbb{R}^{m\times n}$ ($m\ge n$)에 대해

$$
\mathbf{A},\mathbf{P} = \mathbf{Q} \begin{bmatrix} \mathbf{R}*{11} & \mathbf{R}*{12} \ \mathbf{0} & \mathbf{R}\_{22} \end{bmatrix}
$$

와 같이 분해하되, 추가적으로 $\mathbf{R}*{11}$과 $\mathbf{R}*{22}$가 특이값 분해와 유사한 형태의 경계(bound) 조건을 만족하도록 열들을 재배치한다. 즉 $\mathbf{R}*{11}$ 부분이 가능한 한 비정칙성(singular)을 피하고, $\mathbf{R}*{22}$ 부분이 매우 작아지는 형태가 되도록 피벗팅 전략을 최적(또는 준최적)에 가깝게 수행한다.

이 과정을 통해 특정 임계값 $\epsilon$ 이하로 작아진 부분을 0이라 보고 랭크를 추정하거나, 필요 시 그 부분에 대응하는 열들을 제외하는 방식으로 저차원 근사를 얻을 수 있다. 예컨대 대규모 데이터에서 열 간 중복도가 심하거나 잡음이 많은 상태라면, RRQR을 통해 ‘중요도’가 떨어지는 열들을 제거하거나 축소하여 계산량을 현저히 줄일 수 있다. 이러한 랭크저감(Rank Reduction) 기법은 고차원 통계나 머신러닝, 이미지 처리 분야에서 자주 쓰인다.

#### 수정 그람-슈미트(Modified Gram-Schmidt)와 안정성

QR 분해를 구현할 때, 가장 직관적인 방법은 그람-슈미트(Gram-Schmidt) 직교화 과정이다. 이 과정에서 $\mathbf{A} = \[\mathbf{a}\_1 ,\mathbf{a}\_2 ,\dots, \mathbf{a}\_n]$의 열벡터들을 순차적으로 직교화하는 식으로 $\mathbf{Q}$와 $\mathbf{R}$를 형성할 수 있다. 그러나 고전적 그람-슈미트(Classical Gram-Schmidt)는 수치오차에 매우 취약하여, 벡터 간 내적 연산이 커지거나(벡터 크기가 큰 경우) 부동소수점 정밀도가 한계에 달하면, 중간 계산이 틀어지기 쉽다.

이를 보완하기 위해 수정 그람-슈미트(Modified Gram-Schmidt, MGS)가 제안되었다. MGS는 열 하나를 정규화한 뒤, 그 결과를 즉시 다음 열에 반영하면서 직교화 오차가 누적되지 않도록 순차적으로 보정해 나간다. 즉,

$$
\mathbf{q}\_1 = \frac{\mathbf{a}\_1}{|\mathbf{a}\_1|\_2}
$$

를 먼저 구하고, 이어서

$$
\mathbf{a}\_2 \leftarrow \mathbf{a}\_2 - (\mathbf{q}\_1^\top \mathbf{a}\_2),\mathbf{q}\_1
$$

로 제거하고, 그 다음에 $\mathbf{a}\_2$를 정규화하여 $\mathbf{q}\_2$를 얻는다. 그 뒤,

$$
\mathbf{a}\_3 \leftarrow \mathbf{a}\_3 - (\mathbf{q}\_1^\top \mathbf{a}\_3),\mathbf{q}\_1 - (\mathbf{q}\_2^\top \mathbf{a}\_3),\mathbf{q}\_2
$$

와 같이 진행한다. 이때 각 단계마다 직교화를 수행하므로, 고전적 그람-슈미트에 비해 훨씬 더 안정적인 결과를 얻을 수 있다.

그럼에도 불구하고, 하우스홀더 반사(Householder Reflection)에 비하면 MGS조차 수치안정성이 열악할 수 있다. 특히 $m\gg n$인 문제나, 열벡터 간 상관관계가 높아 $\mathrm{rank}(\mathbf{A})$가 애매한 경우(수치적 랭크 부족)에는 하우스홀더가 상대적으로 더 안전하다.

#### 연산 복잡도와 메모리 사용량

* 하우스홀더 QR 분해: $\mathbf{A}\in \mathbb{R}^{m\times n}$에서 $m\ge n$일 때, 일반적으로 $2mn^2 - \tfrac{2}{3}n^3$ 정도의 부동소수점 연산이 필요하다고 분석된다. 메모리 사용 측면에서는 각 단계별로 계수(반사벡터)를 저장해야 하므로, $O(mn)$ 수준의 저장공간이 요구된다.
* Givens QR 분해: Givens 회전을 개별 원소마다 적용하면, 각 소멸연산마다 $O(1)$의 연산이지만 전체적으로 $O(mn^2)$ 정도가 걸린다. 희소행렬에서 0이 많은 위치를 건너뛸 수 있다면 상당히 이점을 얻지만, 밀집(dense) 행렬에는 오히려 비효율적일 수 있다.
* 그람-슈미트(고전 / 수정) 방식: $m\gg n$인 일반적 상황에서는 $O(mn^2)$ 정도의 연산이 든다. 다만 하우스홀더에 비해서는 내적을 많이 수행하기 때문에, 누적 오차가 커질 가능성이 있다. MGS를 이용하면 어느 정도 개선되지만, 여전히 하우스홀더보다 수치안정성이 떨어질 수 있다.
* SVD(특이값 분해): 보통 $O(\min(mn^2,,nm^2))$ 정도의 연산이 필요하므로, 큰 행렬에는 부담이 될 수 있다. 그러나 랭크 결핍 문제나 노이즈가 큰 상황에서는 가장 안전한 해석을 제공하므로, 정밀도를 최우선으로 요구하는 과제에서는 SVD 접근이 각광받는다.

#### 정칙화(Regularization)와 QR

랭크 결핍 혹은 측정 노이즈가 큰 상황에서 직접 최소제곱문제를 풀면, 해가 불안정하게 요동치거나, 계수가 비정상적으로 커져서(과적합 현상) 원하는 근사를 얻지 못하는 경우가 있다. 이를 방지하기 위해 정칙화(regularization) 방법을 사용할 수 있다. 대표적으로 리지 회귀(Ridge Regression)에 해당하는 Tikhonov 정칙화는 다음과 같은 문제를 푼다.

$$
\min\_{\mathbf{x}\in \mathbb{R}^{n}}  |\mathbf{A}\mathbf{x} - \mathbf{b}|\_2^2 + \lambda |\mathbf{x}|\_2^2
$$

여기서 $\lambda>0$는 정칙화 파라미터이다. 정상방정식으로 접근하면

$$
(\mathbf{A}^\top \mathbf{A} + \lambda \mathbf{I}) \mathbf{x} = \mathbf{A}^\top \mathbf{b}
$$

와 같으며, QR 분해 관점에서도 마찬가지로 $\mathbf{A}$를 $\mathbf{Q}\mathbf{R}$로 분해한 뒤, 문제를 적절히 변형하면 해결이 가능하다. 다만, 실제로는 SVD 관점에서

$$
\mathbf{x} = \sum\_{i=1}^r \frac{\sigma\_i}{\sigma\_i^2 + \lambda},\mathbf{u}\_i^\top \mathbf{b},\mathbf{v}\_i
$$

로 명시적인 해석이 가능하여, 정칙화의 효과(작은 특이값 방향을 억제)를 직관적으로 확인할 수 있다. QR 접근에서는 정칙화 항이 $\mathbf{R}$에 직접 가산되지 않는 이상, 추가 단계를 거쳐야 한다. LAPACK이나 ScaLAPACK 같은 고성능 선형대수 패키지에는 이러한 정칙화 문제를 직접 풀어주는 루틴도 존재한다.

#### 곡선 근사(Curve Fitting)에서의 QR 활용

다항회귀(Polynomial Regression)나 지수함수, 로그함수 등 특정한 모형의 곡선 근사를 위해, 일반적으로는 모형함수의 기저를 구성하여 계수추정을 최소제곱으로 푼다. 예를 들어, $f(x)\approx \beta\_0 + \beta\_1 x + \beta\_2 x^2 + \dots + \beta\_{n-1} x^{n-1}$ 꼴의 근사식을 생각하면, $m$개의 측정점 $(x\_i, y\_i)$를 통해

$$
\mathbf{A} =  \begin{bmatrix} 1 & x\_1 & x\_1^2 & \dots & x\_1^{n-1} \ 1 & x\_2 & x\_2^2 & \dots & x\_2^{n-1} \ \vdots & \vdots & \vdots & \ddots & \vdots \ 1 & x\_m & x\_m^2 & \dots & x\_m^{n-1} \end{bmatrix}, \quad \mathbf{b} =  \begin{bmatrix} y\_1 \ y\_2 \ \vdots \ y\_m \end{bmatrix}
$$

형태로 시스템을 만든 뒤, 이를 QR 분해로 푼다. 하우스홀더 반사나 Givens 회전을 이용하면 내적의 증폭에 의한 오차 문제가 줄어들어, 수치적으로 안정적인 계수를 얻을 수 있다.

또한, $x\_i$가 큰 값들이거나 매우 작은 값들일 때, $x\_i^k$ 꼴이 매우 큰 스케일차이를 일으켜 수치오차를 유발할 수 있는데, QR 분해 접근은 정상방정식 대비 상대적으로 오차에 잘 대응한다. 그러나 $n$이 커질수록 다항식 기저가 서로 선형의존에 가까워지는 문제가 발생할 수 있으므로, 적절히 정규직교 다항식을 이용하거나(예: 체비쇼프 다항식), 다른 기법(Splines 등)과 병행할 수도 있다.

#### 전체 정리

QR 분해는 정상방정식 접근보다 훨씬 안정적이면서도, SVD만큼의 계산량은 들이지 않고 최소제곱문제를 풀 수 있는 효율적 방법이다. 하우스홀더, Givens, (수정) 그람-슈미트 등 다양한 구현 방식이 있으며, 실제로는 하우스홀더가 가장 폭넓게 쓰인다. 랭크 결핍이 우려되면 컬럼 피벗팅 또는 Rank-Revealing QR(RRQR) 기법을 사용하여 안정도를 높일 수 있고, 가중 최소제곱이나 정칙화에도 무리 없이 응용 가능하다.