# 선형 회귀와 다항 회귀

#### 선형 회귀 기초

회귀(regression)는 주어진 입력과 출력 간의 관계를 추정하는 과정이다. 최소제곱법(Least Squares)은 오차 제곱합을 최소화하도록 회귀 함수를 정하는 대표적인 방법이며, 데이터가 주어졌을 때 함숫값과 실제 관측값의 차이를 나타내는 오차(error)가 가장 작아지도록 계수(coefficients)를 찾는다. 이를 통해 예측에 사용하거나, 데이터의 경향을 파악하고자 한다.

선형 회귀는 종속 변수 $y$와 독립 변수 $x$ 간에 선형 관계를 가정한다. 즉, 입력 $x$가 한 개일 경우 다음과 같이 표현할 수 있다.

$$
\begin{align} y \approx \beta\_0 + \beta\_1 x \end{align}
$$

이때 $\beta\_0$는 상수항(intercept)이고, $\beta\_1$은 기울기(slope)이다. 실제로는 관측 과정이나 모델링 과정에서 오차가 발생하므로, 다음과 같이 오차항 $\epsilon$을 추가한 형태로 표현한다.

$$
\begin{align} y = \beta\_0 + \beta\_1 x + \epsilon \end{align}
$$

하지만 일반적으로 $\epsilon$은 관측 불확실성을 나타내므로 직접 해석하기보다는, 여러 관측 데이터 ${ (x\_i, y\_i)}$를 통해 $\beta\_0, \beta\_1$를 추정한다.

#### 다변량 선형 회귀 일반화

독립 변수가 $x\_1, x\_2, \dots, x\_n$과 같이 여러 개 존재한다면, 이를 벡터 $\mathbf{x} \in \mathbb{R}^n$로 표현하여 다음과 같은 선형 모델을 세울 수 있다.

$$
\begin{align} y \approx \beta\_0 + \beta\_1 x\_1 + \beta\_2 x\_2 + \cdots + \beta\_n x\_n \end{align}
$$

이는 다음과 같은 벡터 표현으로도 나타낼 수 있다. 편의상 상수항을 포함시키기 위해 독립 변수 벡터에 1을 추가한 $\mathbf{x} = (1, x\_1, x\_2, \dots, x\_n)^\top$로 정의하고, 계수 벡터를 $\boldsymbol{\beta} = (\beta\_0, \beta\_1, \dots, \beta\_n)^\top$라 할 때,

$$
\begin{align} y \approx \mathbf{x}^\top \boldsymbol{\beta} \end{align}
$$

이제 관측 데이터가 $m$개 주어졌다고 하자. 각 관측값을 $y\_i$라 하고, 대응하는 독립 변수 벡터를 $\mathbf{x}\_i$라 하면, 전체 시스템을 다음과 같은 형태로 표현할 수 있다.

$$
\begin{align} \mathbf{y} = \begin{pmatrix} y\_1 \ y\_2 \ \vdots \ y\_m \end{pmatrix}, \quad \mathbf{X} = \begin{pmatrix} 1 & x\_{1,1} & x\_{1,2} & \cdots & x\_{1,n} \ 1 & x\_{2,1} & x\_{2,2} & \cdots & x\_{2,n} \ \vdots & \vdots & \vdots & & \vdots \ 1 & x\_{m,1} & x\_{m,2} & \cdots & x\_{m,n} \end{pmatrix}, \quad \boldsymbol{\beta} = \begin{pmatrix} \beta\_0 \ \beta\_1 \ \vdots \ \beta\_n \end{pmatrix} \end{align}
$$

이때 관측값 $\mathbf{y}$와 모델 예측값 $\mathbf{X}\boldsymbol{\beta}$의 차이를 최소화하는 $\boldsymbol{\beta}$를 찾는 것이 최소제곱법의 목표다.

#### 최소제곱 해(정규방정식)

최소제곱법에서 오차의 제곱합을 다음과 같이 정의한다.

$$
\begin{align} S(\boldsymbol{\beta})  = |\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|\_2^2  = (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \end{align}
$$

이 함수를 $\boldsymbol{\beta}$에 관해 최소화하기 위해, $\boldsymbol{\beta}$에 대한 편미분이 0이 되는 조건을 구한다.

$$
\begin{align} \frac{\partial S}{\partial \boldsymbol{\beta}} = \frac{\partial}{\partial \boldsymbol{\beta}}  \Big( (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \Big) = -2 \mathbf{X}^\top (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})  = \mathbf{0} \end{align}
$$

이를 정규방정식(normal equation)이라 부르며, 아래와 같이 정리할 수 있다.

$$
\begin{align} \mathbf{X}^\top \mathbf{X} \boldsymbol{\beta} = \mathbf{X}^\top \mathbf{y} \end{align}
$$

만일 $\mathbf{X}^\top\mathbf{X}$가 역행렬을 갖는다면, 즉 보통 독립 변수가 선형종속이 아니고 표본 크기가 충분할 경우, 최소제곱 해는

$$
\begin{align} \boldsymbol{\beta} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y} \end{align}
$$

가 된다.

#### 선형 모델의 해석

이렇게 구한 $\boldsymbol{\beta}$는 데이터에 가장 잘 들어맞는(오차 제곱합이 최소가 되는) 파라미터 세트이다. 실제 문제에서는 $\mathbf{X}^\top \mathbf{X}$가 직접적으로 역행렬을 갖지 않거나(수치적으로 매우 불안정), 독립 변수들 간의 다중공선성(multicollinearity)이 존재할 수도 있다. 이러한 문제가 심각해지면 모델의 해가 불안정하거나, 예측력이 저하될 수 있다. 따라서 현실에서는 데이터 전처리, 변수 선택, 정규화(regularization) 등 다양한 기법을 병행하여 모델의 안정성을 높인다.

#### 다항 회귀 기초

선형 회귀에서 말하는 "선형성"은 계수 $\boldsymbol{\beta}$에 대해 선형임을 의미한다. 따라서 독립 변수(특히 단일 변수 $x$)를 다항식 형태로 변형한 뒤, 선형 회귀를 적용하면 다항 회귀(polynomial regression)를 얻을 수 있다. 예를 들어 단일 변수 $x$에 대해 2차 다항식을 가정하면,

$$
\begin{align} y \approx \beta\_0 + \beta\_1 x + \beta\_2 x^2 \end{align}
$$

더 나아가 $d$차 다항식을 고려하면,

$$
\begin{align} y \approx \beta\_0 + \beta\_1 x + \beta\_2 x^2 + \cdots + \beta\_d x^d \end{align}
$$

이 상황을 일반적인 선형 회귀 모형으로 해석하기 위해, 입력 변수 $x$를 확장된 벡터 $\mathbf{x} = (1, x, x^2, \dots, x^d)^\top$로 정의한다면, 계수 벡터 $\boldsymbol{\beta} = (\beta\_0, \beta\_1, \dots, \beta\_d)^\top$에 대해

$$
\begin{align} y \approx \mathbf{x}^\top \boldsymbol{\beta} \end{align}
$$

로 볼 수 있다. 즉, $\mathbf{x}$라는 독립 변수를 다항식 형태로 확장만 해주면, 여전히 계수 벡터에 대해 선형적인 모형이 되므로 최소제곱법을 그대로 적용할 수 있다.

#### 다항 회귀에서의 설계행렬

관측 데이터가 $m$개, 그리고 다항 차수가 $d$라고 하자. 각 관측점은 $(x\_i, y\_i)$이며, 이를 통해 다음과 같은 $\mathbf{X}$와 $\mathbf{y}$를 구성할 수 있다.

$$
\begin{align} \mathbf{X} = \begin{pmatrix} 1 & x\_1 & x\_1^2 & \cdots & x\_1^d \ 1 & x\_2 & x\_2^2 & \cdots & x\_2^d \ \vdots & \vdots & \vdots & & \vdots \ 1 & x\_m & x\_m^2 & \cdots & x\_m^d \end{pmatrix}, \quad \mathbf{y} = \begin{pmatrix} y\_1 \ y\_2 \ \vdots \ y\_m \end{pmatrix} \end{align}
$$

이 구성으로부터 앞서 살펴본 정규방정식을 이용해 $\boldsymbol{\beta}$를 추정하면, 결국은 다항 회귀 모델 계수를 얻을 수 있다.

#### 차수 선택과 과적합

다항 회귀에서 가장 중요한 문제 중 하나는 다항 차수 $d$의 선택이다. 차수가 너무 낮으면 실제 데이터의 복잡한 패턴을 반영하지 못하여 과소적합(underfitting)이 발생할 수 있다. 반대로 차수가 지나치게 높으면 모든 관측점을 '맞추려고' 하여 노이즈까지 과하게 학습하는 과적합(overfitting)이 발생할 수 있다. 이를 방지하려면 외부 데이터(검증 데이터)를 통한 성능 평가나, 교차 검증(cross-validation) 기법 등을 활용하여 적절한 차수를 결정한다.

#### 정규화(Regularization)의 개념

차수가 높아지고 모델이 복잡해질수록, 계수 벡터 $\boldsymbol{\beta}$가 급격히 커질 수 있으며 이는 예측의 불안정성으로 이어질 수 있다. 이런 문제를 완화하기 위해서는 정규화(또는 규제) 기법이 활용된다. 예를 들어 다음과 같은 릿지 회귀(ridge regression)는 계수들의 제곱합에 대한 페널티 항을 추가하여 모델이 계수를 지나치게 크게 갖지 않도록 한다.

정규화 항 $\lambda |\boldsymbol{\beta}|\_2^2$을 추가한 비용함수는

$$
\begin{align} S\_{\text{ridge}}(\boldsymbol{\beta}) = |\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|\_2^2 + \lambda |\boldsymbol{\beta}|\_2^2 \end{align}
$$

이 경우 정규방정식은 조금 달라지며,

$$
\begin{align} (\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I}) \boldsymbol{\beta} = \mathbf{X}^\top \mathbf{y} \end{align}
$$

로 나타난다. 여기서 $\mathbf{I}$는 단위 행렬이다. 적절한 $\lambda$값을 선택함으로써, 과적합을 제어하고 모델의 예측 성능을 높일 수 있다.

#### 고급 회귀 기법 및 다항 회귀 심화

회귀 모델을 설계할 때, 단순 선형 회귀와 다항 회귀만으로는 해결하기 어려운 상황이 자주 발생한다. 이를테면 데이터에 잡음이 많거나, 독립 변수 간의 상관관계가 높은 경우, 차수가 높아 모델이 복잡해지는 경우 등이 있다. 이럴 때 정규화(regularization), 가중 최소제곱법(weighted least squares), 직교 다항식(orthogonal polynomials) 등 다양한 기법을 통해 모델의 안정성과 예측 성능을 높일 수 있다.

#### Lasso 회귀

릿지 회귀(ridge regression)는 계수 벡터의 2-노름에 대한 페널티 항을 추가하여 모형이 계수를 크게 갖지 못하도록 규제한다. 반면 Lasso 회귀(least absolute shrinkage and selection operator)는 다음과 같이 계수 벡터의 1-노름에 대한 페널티 항을 추가한다.

Slasso(β)=∥y−Xβ∥22+λ∥β∥1\begin{align} S\_{\text{lasso}}(\boldsymbol{\beta}) = |\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|\_2^2 + \lambda |\boldsymbol{\beta}|\_1 \end{align}

이때 $|\boldsymbol{\beta}|\_1$은 계수들의 절댓값을 모두 더한 것이므로, 일부 계수가 정확히 0이 되도록 하는 희소성(sparsity)을 유도한다. 즉 Lasso 회귀는 변수 선택(variable selection) 효과가 있으며, 고차 다항 회귀 모델에서 중요하지 않은 항들을 제거해 모델 해석을 보다 간결하게 만든다. 단, 1-노름 규제 항으로 인해 미분 가능성이 깨지므로, 경사 하강(gradient descent)이나 좌표 하강법(coordinate descent) 등의 수치 알고리즘을 사용해 최적 해를 구한다.

#### Elastic Net

릿지 회귀와 Lasso 회귀를 동시에 적용하는 방법이 Elastic Net이다. 계수 벡터의 2-노름 항과 1-노름 항을 모두 포함한 페널티를 가짐으로써, 두 방식의 장점을 적절히 절충한다. 비용 함수는 다음과 같다.

SEN(β)=∥y−Xβ∥22+αλ∥β∥1+(1−α)λ∥β∥22\begin{align} S\_{\text{EN}}(\boldsymbol{\beta}) = |\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|\_2^2 + \alpha \lambda |\boldsymbol{\beta}|\_1 + (1-\alpha) \lambda |\boldsymbol{\beta}|\_2^2 \end{align}

여기서 $\lambda$는 전체 규제 강도, $\alpha$는 Lasso와 릿지 간의 균형을 결정하는 하이퍼파라미터다. Elastic Net은 다항 차수가 높고 변수 간 상관관계가 클 때 유용하며, 변수 선택 및 일반화 성능 사이의 균형을 맞추는 데 적합하다.

#### 가중 최소제곱법(Weighted Least Squares)

일부 상황에서는 모든 관측값에 동일한 비중을 두는 것이 합리적이지 않을 수 있다. 예를 들어 측정값 중 특정 구간에서 오차가 커질 가능성이 높다면, 각 관측값에 가중치 $w\_i$를 부여한 다음, 다음과 같은 가중 잔차 제곱합을 최소화하는 방식으로 회귀를 수행할 수 있다.

Sw(β)=(y−Xβ)⊤W(y−Xβ)\begin{align} S\_w(\boldsymbol{\beta}) = (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top \mathbf{W} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \end{align}

이때 $\mathbf{W}$는 대각 행렬이며, 그 대각 원소가 각 관측값의 가중치 $w\_i$로 구성된다. 정규방정식도 다음과 같이 변형된다.

X⊤WXβ=X⊤Wy\begin{align} \mathbf{X}^\top \mathbf{W} \mathbf{X} \boldsymbol{\beta} = \mathbf{X}^\top \mathbf{W} \mathbf{y} \end{align}

이는 다항 회귀에도 동일하게 적용될 수 있으며, 오차 분산이 균등하지 않은 상황(heteroskedasticity)에서 모델 예측 성능과 추정의 효율을 높인다.

#### 직교 다항식(Orthogonal Polynomials)

다항 차수가 커질수록, 단순한 다항 기저 $1, x, x^2, \dots, x^d$로 구성된 행렬 $\mathbf{X}$는 심각한 수치적 불안정(ill-conditioning)에 직면할 수 있다. 예를 들어 $x^k$ 항들이 매우 큰 값과 작은 값을 동시에 가지면, $\mathbf{X}^\top \mathbf{X}$의 조건수가 크게 악화된다. 이를 완화하기 위해 체비쇼프(Chebyshev), 레제느(Legendre) 등의 직교 다항식을 사용하여 기저 벡터를 구성하기도 한다.

직교 다항식 기저 $\phi\_k(x)$는 서로 다른 $k$에 대해 적분 내적이 0이 되므로, 기저 간 상관관계가 낮다. 결국 수치 해석적인 안정성이 좋아지고, 회귀 계수를 추정할 때 역행렬 계산이 덜 민감하게 된다. 다만, 실제 관측 구간에 맞추어 직교화 과정(예: 그라마-슈미트(Gram–Schmidt) 과정)을 거쳐야 할 수도 있다.

#### 다항 피처의 표준화와 다중공선성

다항 회귀는 입력 변수를 $x, x^2, x^3, \dots, x^d$처럼 확장하므로, 차수가 높아질수록 다중공선성(multicollinearity) 문제가 심화된다. $x^2$와 $x^3$ 등의 항들이 서로 강한 상관관계를 가지면, $\mathbf{X}^\top \mathbf{X}$가 거의 특이(singular)에 가까워지거나 수치적으로 불안정해질 수 있다. 이를 개선하기 위해서는 변수들을 표준화(standardization)하거나, 중심화(mean-centering)하는 등의 전처리를 수행한다.

예를 들어 $x$의 평균을 빼고, 필요하다면 분산을 1로 맞추면(표준화), 다항 항들의 크기 스케일이 줄어들어 상호 상관관계가 완화되고, 모델 추정 과정의 안정성이 높아진다.

#### 오차 분석과 잔차 분석(Residual Analysis)

최소제곱법으로 회귀 모델을 적합한 뒤에는 잔차(residual)를 통해 모델의 적합도를 평가하고, 추정 가정(선형성, 등분산성, 정규성 등)이 어느 정도 성립하는지 살펴본다. 일반적으로는 잔차를 그래프로 시각화하여 패턴이나 이상값(outliers)이 존재하는지 확인한다. 특히 다항 회귀의 경우, 잔차가 특정 구간에서 체계적인 패턴을 보인다면 모델 차수가 부족하거나, 특정 항을 추가/제거해야 할 수도 있다.

#### 간단한 Python 예제

다항 회귀 모델을 예시로 보인다. 아래 예제는 Python 언어와 scikit-learn 라이브러리를 이용한다.

```python
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

# 예시 데이터
x = np.linspace(-2, 2, 10).reshape(-1, 1)
y = 1.2 + 2.5 * x + 0.8 * x**2 + np.random.randn(10, 1) * 0.3

# 2차 다항 변환
poly = PolynomialFeatures(degree=2, include_bias=True)
X_poly = poly.fit_transform(x)

# 선형 회귀 적합
model = LinearRegression()
model.fit(X_poly, y)

print("추정된 계수:", model.coef_)
print("추정된 상수항:", model.intercept_)
```

PolynomialFeatures로부터 2차 다항식 기저항을 생성하였고, 이를 선형 회귀에 적용하였다. scikit-learn 내부적으로 최소제곱법(정규방정식 혹은 QR 분해 등)을 사용하여 계수를 추정한다.

#### 다항 회귀와 다항 보간의 비교

다항식을 이용하여 데이터에 근사하는 접근 방식은 크게 회귀(regression)와 보간(interpolation)으로 나눌 수 있다. 다항 회귀는 오차 제곱합을 최소화하는 해를 구하여 데이터를 "근사"하는 데 초점을 맞추는 반면, 다항 보간은 모든 관측점을 정확히 지나는 다항식을 찾는 목적을 갖는다. 예를 들어 $m$개의 점 $(x\_i, y\_i)$에 대해 정확히 통과하는 다항식 $p(x)$를 찾기 위해서는 $m-1$ 차 이하의 다항식으로 충분하며, 이는 이론상 유일하다.

그러나 실제로는 관측값에 노이즈가 존재하기 때문에, 모든 점을 정확히 통과하는 다항식(고차 다항식일수록 더욱)은 과적합을 유발하기 쉽다. 반면 최소제곱법으로 구한 다항 회귀는 모든 점을 정확히 통과할 필요 없이, 전체적으로 관측값과의 제곱 오차가 작아지도록 계수를 추정한다. 따라서 노이즈가 있는 상황에서는 과적합이 일어나기 쉬운 보간보다는, 다항 회귀가 일반적으로 더 유연하고 안정적이다.

#### 스플라인(Spline) 기법

다항 회귀가 단일 전역 다항식을 사용한다면, 스플라인은 구간별로 낮은 차수(일반적으로 3차 이하) 다항식을 연결하되, 경계점에서 매끄럽게 이어지도록 하는 방법이다. 특히 큐빅 스플라인(cubic spline)은 많이 사용된다. 전역 다항식을 이용하면 차수가 높아질수록 수치적 문제나 과적합이 생기기 쉬우나, 스플라인은 구간별로 다항식을 정의하므로 이런 문제를 줄이면서 유연성을 확보한다.

다항 회귀와 직접적인 비교를 하자면, 스플라인은 "특정 구간" 내에서만 다항 함수를 정의하므로 각각의 구간에서 오차 제곱합을 줄이는 최적화를 수행할 수 있다. 이로써 전체적인 복잡도는 상대적으로 낮게 유지하면서도, 자료가 많은 구간 혹은 급변하는 구간에 대해선 더 잘 적응한다. 다만 스플라인을 구성하기 위해서는 구간 분할점(knots)을 어떻게 배치할지 결정하는 과정이 필요하다.

#### 다항 회귀와 커널 방법

선형 회귀로 시작하여 다항(폴리노미얼) 확장을 하면, 결국 $(x, x^2, \dots, x^d)$와 같은 형태의 고차 피처 공간(feature space)에서 선형 모델을 적용하는 셈이 된다. 커널 트릭(kernel trick)을 이용하면, 이러한 고차원 확장을 직접 하지 않고서도 마치 다항 함수를 사용한 것과 같은 효과를 얻을 수 있다. 예를 들어 서포트 벡터 머신(SVM)에서 다항 커널(polynomial kernel)을 사용하면, 계산 효율성 면에서 매우 큰 이점이 있을 수 있다.

다항 커널의 예로서, 차수 $d$에 대해

$$
K(\mathbf{x}, \mathbf{z}) = (\mathbf{x}^\top \mathbf{z} + c)^d
$$

와 같은 식이 있다. 이는 $\mathbf{x}$와 $\mathbf{z}$가 $\mathbb{R}^n$에 있을 때, 사실상 $\mathbf{x}$, $\mathbf{z}$의 모든 단항·이항·고차 항을 포함하는 고차원 매핑과 동등한 효과를 낸다. 다만, 회귀 문제에서 이 커널을 직접 사용하는 것은 주로 커널 회귀(Kernel Ridge Regression, SVR 등)를 통해 이루어지며, 전형적인 다항 회귀와는 구현 방식이 다르다.

#### 로버스트 회귀(Robust Regression)

최소제곱법은 모든 오차의 제곱을 더하여 이를 최소화하므로, 이상값(outlier)이 존재하면 그 영향이 매우 커진다. 이러한 상황에서 다항 회귀 역시 이상값에 쉽게 휘둘릴 수 있다. 예를 들어 극단적인 $x$ 범위에서 잘못된 측정값이 포함되어 있다면, 고차 항의 계수가 크게 왜곡될 수 있다.

이를 완화하기 위해서는 가중 최소제곱법과는 다른 접근인 로버스트 회귀 기법을 사용할 수 있다. 예를 들어 M-추정(M-estimation)은 제곱 오차 대신 적절한 완화 함수(예: Huber 손실, Hampel 손실 등)를 사용하여 이상값의 영향을 줄인다. 간단히 말해, 일정 크기 이상의 오차가 발생하면, 그 오차에 대한 "가중치"를 자동으로 낮춰 이상값이 전체 모델에 큰 영향을 주지 못하게 하는 방식이다.

#### 차수 확장 시 주의사항

다항 회귀에서 차수를 높이면 통계적 자유도가 증가하여, 모델이 훈련 데이터에 과도하게 적합될 위험이 있다. 이 때문에 차수를 무작정 늘리기보다, 교차 검증(cross-validation) 같은 기법으로 모델의 성능 변화를 관찰하면서 최적 차수를 결정하는 것이 일반적이다. 더불어, 분산이 큰 데이터를 다항으로 과도하게 맞추면 곡선이 지나치게 요동하게 되므로, 데이터의 본질적 구조를 파악하는 데에도 어려움이 생길 수 있다.

#### 교차 검증과 모델 평가

최소제곱법으로 학습된 다항 회귀 모델이 과적합에 빠지지 않았는지 확인하려면, 학습 데이터와는 별도로 검증 데이터 혹은 테스트 데이터를 준비하여 예측 정확도를 평가해야 한다. 교차 검증(k-fold cross-validation)은 데이터가 제한적일 때도 모델의 예측 성능을 안정적으로 측정할 수 있는 대표적인 방법이다. 이를 통해 차수 선택, 정규화 파라미터 $\lambda$ 등의 하이퍼파라미터를 조정하고, 실제로 데이터가 추가되었을 때 모델이 얼마나 잘 일반화할 수 있는지 판단한다.

#### Octave 예시

다항 회귀를 Octave로 수행하는 예시이다. 3차 다항식을 가정하고, QR 분해를 사용하여 계수를 추정한다.

```octave
% 예시 데이터
x = linspace(-1, 1, 10)';
y = 1.5 + 2.0*x + 0.7*x.^2 - 0.2*x.^3 + 0.3*randn(size(x));

% 설계행렬 구성
X = [ones(size(x)), x, x.^2, x.^3];

% QR 분해를 이용한 최소제곱 해
[Q,R] = qr(X,0);
beta = R \ (Q' * y);

disp('추정된 계수(상수항, x^1, x^2, x^3 순서):');
disp(beta);
```

여기에서 `X`는 $\[1, x, x^2, x^3]$ 형태로 구성된 설계행렬이며, `qr(X,0)`는 축소 QR 분해(thin QR decomposition)를 수행한다. 이후 $R \ (Q' \* y)$ 연산으로 $\boldsymbol{\beta}$를 구하는 것은, 정규방정식을 $R$과 $Q$를 통해 안정적으로 푸는 과정과 동일하다.