6.101 재귀적 최소 제곱법과 온라인 추정

6.101 재귀적 최소 제곱법과 온라인 추정

1. 개요

로봇 시스템의 파라미터 추정에서는 새로운 데이터가 실시간으로 수집될 때마다 추정값을 갱신해야 하는 경우가 빈번하다. 배치(batch) 최소 제곱법은 전체 데이터를 한꺼번에 처리하므로, 데이터가 추가될 때마다 전체 계산을 반복해야 한다. 재귀적 최소 제곱법(Recursive Least Squares, RLS)은 이전 시점의 추정 결과를 활용하여 새로운 데이터를 효율적으로 반영하는 온라인 추정 기법이다. 이 절에서는 RLS의 수학적 유도, 수렴 성질, 그리고 로봇공학에서의 적용 사례를 다룬다.


2. 배치 최소 제곱법의 복습

k개의 관측이 주어진 시점에서 선형 모델

y_i = \mathbf{a}_i^\top \boldsymbol{\theta} + e_i, \quad i = 1, 2, \dots, k

를 행렬 형태로 쓰면

\mathbf{y}_k = \mathbf{A}_k \boldsymbol{\theta} + \mathbf{e}_k

이다. 여기서 \mathbf{A}_k \in \mathbb{R}^{k \times n}는 회귀 행렬, \mathbf{y}_k \in \mathbb{R}^k는 출력 벡터, \boldsymbol{\theta} \in \mathbb{R}^n는 추정할 파라미터 벡터, \mathbf{e}_k는 오차 벡터이다. 배치 최소 제곱 추정은 다음과 같다.

\hat{\boldsymbol{\theta}}_k = (\mathbf{A}_k^\top \mathbf{A}_k)^{-1} \mathbf{A}_k^\top \mathbf{y}_k

정보 행렬(information matrix)을 \mathbf{P}_k^{-1} = \mathbf{A}_k^\top \mathbf{A}_k로 정의하면

\hat{\boldsymbol{\theta}}_k = \mathbf{P}_k \mathbf{A}_k^\top \mathbf{y}_k

이다.


3. 재귀적 최소 제곱법의 유도

3.1 기본 아이디어

k번째 시점까지의 추정값 \hat{\boldsymbol{\theta}}_k가 주어졌을 때, (k+1)번째 관측 (y_{k+1}, \mathbf{a}_{k+1})이 추가되면 추정값을 다음과 같이 갱신한다.

\hat{\boldsymbol{\theta}}_{k+1} = \hat{\boldsymbol{\theta}}_k + \mathbf{K}_{k+1} (y_{k+1} - \mathbf{a}_{k+1}^\top \hat{\boldsymbol{\theta}}_k)

여기서 \mathbf{K}_{k+1} \in \mathbb{R}^n는 이득 벡터(gain vector)이며, y_{k+1} - \mathbf{a}_{k+1}^\top \hat{\boldsymbol{\theta}}_k는 혁신(innovation) 또는 예측 오차이다.

3.2 정보 행렬의 갱신

(k+1)번째 시점의 정보 행렬은 다음과 같다.

\mathbf{P}_{k+1}^{-1} = \mathbf{A}_{k+1}^\top \mathbf{A}_{k+1} = \mathbf{A}_k^\top \mathbf{A}_k + \mathbf{a}_{k+1} \mathbf{a}_{k+1}^\top = \mathbf{P}_k^{-1} + \mathbf{a}_{k+1} \mathbf{a}_{k+1}^\top

3.3 공분산 행렬의 갱신: 행렬 역보조 정리

행렬 역보조 정리(matrix inversion lemma, Woodbury identity)를 적용하면

(\mathbf{M} + \mathbf{u}\mathbf{v}^\top)^{-1} = \mathbf{M}^{-1} - \frac{\mathbf{M}^{-1} \mathbf{u} \mathbf{v}^\top \mathbf{M}^{-1}}{1 + \mathbf{v}^\top \mathbf{M}^{-1} \mathbf{u}}

여기서 \mathbf{M} = \mathbf{P}_k^{-1}, \mathbf{u} = \mathbf{v} = \mathbf{a}_{k+1}로 두면

\mathbf{P}_{k+1} = \mathbf{P}_k - \frac{\mathbf{P}_k \mathbf{a}_{k+1} \mathbf{a}_{k+1}^\top \mathbf{P}_k}{1 + \mathbf{a}_{k+1}^\top \mathbf{P}_k \mathbf{a}_{k+1}}

3.4 이득 벡터의 유도

이득 벡터를 다음과 같이 정의한다.

\mathbf{K}_{k+1} = \frac{\mathbf{P}_k \mathbf{a}_{k+1}}{1 + \mathbf{a}_{k+1}^\top \mathbf{P}_k \mathbf{a}_{k+1}}

이를 사용하면 공분산 갱신은 다음과 같이 간결하게 표현된다.

\mathbf{P}_{k+1} = (\mathbf{I} - \mathbf{K}_{k+1} \mathbf{a}_{k+1}^\top) \mathbf{P}_k

3.5 파라미터 갱신의 유도

\hat{\boldsymbol{\theta}}_{k+1} = \mathbf{P}_{k+1} \mathbf{A}_{k+1}^\top \mathbf{y}_{k+1}에서 출발하여 전개하면

\begin{aligned} \hat{\boldsymbol{\theta}}_{k+1} &= \mathbf{P}_{k+1} (\mathbf{A}_k^\top \mathbf{y}_k + \mathbf{a}_{k+1} y_{k+1}) \\ &= \mathbf{P}_{k+1} (\mathbf{P}_k^{-1} \hat{\boldsymbol{\theta}}_k + \mathbf{a}_{k+1} y_{k+1}) \\ &= \mathbf{P}_{k+1} \mathbf{P}_k^{-1} \hat{\boldsymbol{\theta}}_k + \mathbf{P}_{k+1} \mathbf{a}_{k+1} y_{k+1} \end{aligned}

\mathbf{P}_{k+1} \mathbf{P}_k^{-1} = \mathbf{I} - \mathbf{K}_{k+1} \mathbf{a}_{k+1}^\top임을 이용하고, \mathbf{P}_{k+1} \mathbf{a}_{k+1} = \mathbf{K}_{k+1}임을 확인하면

\hat{\boldsymbol{\theta}}_{k+1} = \hat{\boldsymbol{\theta}}_k + \mathbf{K}_{k+1}(y_{k+1} - \mathbf{a}_{k+1}^\top \hat{\boldsymbol{\theta}}_k)

이 성립한다.


4. RLS 알고리즘 요약

초기화:

  • \hat{\boldsymbol{\theta}}_0: 초기 파라미터 추정값 (보통 \mathbf{0})
  • \mathbf{P}_0 = \alpha \mathbf{I} (\alpha \gg 1): 초기 공분산 행렬

매 시간 단계 k = 0, 1, 2, \dots에서:

  1. 새로운 관측 (y_{k+1}, \mathbf{a}_{k+1})을 수신한다.
  2. 이득 벡터를 계산한다.

\mathbf{K}_{k+1} = \frac{\mathbf{P}_k \mathbf{a}_{k+1}}{1 + \mathbf{a}_{k+1}^\top \mathbf{P}_k \mathbf{a}_{k+1}}

  1. 파라미터를 갱신한다.

\hat{\boldsymbol{\theta}}_{k+1} = \hat{\boldsymbol{\theta}}_k + \mathbf{K}_{k+1}(y_{k+1} - \mathbf{a}_{k+1}^\top \hat{\boldsymbol{\theta}}_k)

  1. 공분산 행렬을 갱신한다.

\mathbf{P}_{k+1} = (\mathbf{I} - \mathbf{K}_{k+1} \mathbf{a}_{k+1}^\top) \mathbf{P}_k


5. 가중 재귀적 최소 제곱법

관측의 신뢰도가 서로 다를 때, 가중 계수 w_{k+1} > 0을 도입한 가중 RLS는 다음과 같다.

\mathbf{K}_{k+1} = \frac{\mathbf{P}_k \mathbf{a}_{k+1}}{w_{k+1}^{-1} + \mathbf{a}_{k+1}^\top \mathbf{P}_k \mathbf{a}_{k+1}}

\hat{\boldsymbol{\theta}}_{k+1} = \hat{\boldsymbol{\theta}}_k + \mathbf{K}_{k+1}(y_{k+1} - \mathbf{a}_{k+1}^\top \hat{\boldsymbol{\theta}}_k)

\mathbf{P}_{k+1} = (\mathbf{I} - \mathbf{K}_{k+1} \mathbf{a}_{k+1}^\top) \mathbf{P}_k


6. 망각 인자를 사용한 RLS

시변(time-varying) 파라미터를 추적하기 위해 망각 인자(forgetting factor) 0 < \mu \leq 1을 도입한다. 이 경우 비용 함수는 다음과 같다.

J_k(\boldsymbol{\theta}) = \sum_{i=1}^{k} \mu^{k-i} (y_i - \mathbf{a}_i^\top \boldsymbol{\theta})^2

\mu < 1이면 오래된 데이터에 지수적으로 감소하는 가중치를 부여하여 최근 데이터를 더 중시한다. 갱신식은 다음과 같이 수정된다.

\mathbf{K}_{k+1} = \frac{\mathbf{P}_k \mathbf{a}_{k+1}}{\mu + \mathbf{a}_{k+1}^\top \mathbf{P}_k \mathbf{a}_{k+1}}

\hat{\boldsymbol{\theta}}_{k+1} = \hat{\boldsymbol{\theta}}_k + \mathbf{K}_{k+1}(y_{k+1} - \mathbf{a}_{k+1}^\top \hat{\boldsymbol{\theta}}_k)

\mathbf{P}_{k+1} = \frac{1}{\mu} (\mathbf{I} - \mathbf{K}_{k+1} \mathbf{a}_{k+1}^\top) \mathbf{P}_k

\mu의 전형적인 값은 0.95 \leq \mu \leq 0.999 범위이다. \mu = 1이면 표준 RLS와 동일하다.


7. 수렴 성질

7.1 일관성 (Consistency)

관측 잡음 e_i가 평균 0, 분산 \sigma^2인 i.i.d. 확률 변수이고 회귀 벡터 \mathbf{a}_i가 충분히 풍부하게(persistently exciting) 변화하면, RLS 추정값은 참값에 확률적으로 수렴한다.

\hat{\boldsymbol{\theta}}_k \xrightarrow{k \to \infty} \boldsymbol{\theta}^* \quad \text{(확률 1)}

7.2 지속적 가진 조건 (Persistent Excitation)

RLS의 수렴을 보장하기 위해서는 다음 지속적 가진 조건이 필요하다. 양의 상수 \alpha_1, \alpha_2 > 0과 정수 N > 0이 존재하여 모든 k에 대해

\alpha_1 \mathbf{I} \leq \sum_{i=k+1}^{k+N} \mathbf{a}_i \mathbf{a}_i^\top \leq \alpha_2 \mathbf{I}

가 성립하여야 한다. 이 조건은 입력 신호가 충분히 다양한 방향으로 변화함을 의미한다.


8. 로봇공학에서의 응용

8.1 동역학 파라미터 추정

로봇 매니퓰레이터의 역동역학 모델은 관절 토크 \boldsymbol{\tau}에 대해 다음과 같이 선형 파라미터 형태로 표현된다.

\boldsymbol{\tau} = \mathbf{Y}(\mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}}) \boldsymbol{\theta}

여기서 \mathbf{Y}는 회귀 행렬(regressor matrix), \boldsymbol{\theta}는 관성 파라미터(질량, 관성 모멘트, 질량 중심 위치 등)를 포함하는 벡터이다. RLS를 적용하면 로봇 운용 중에 실시간으로 파라미터를 추정할 수 있다.

매 제어 주기 k에서 관절 각도 \mathbf{q}_k, 각속도 \dot{\mathbf{q}}_k, 각가속도 \ddot{\mathbf{q}}_k, 토크 \boldsymbol{\tau}_k를 측정하여

\hat{\boldsymbol{\theta}}_{k+1} = \hat{\boldsymbol{\theta}}_k + \mathbf{K}_{k+1}(\boldsymbol{\tau}_{k+1} - \mathbf{Y}_{k+1} \hat{\boldsymbol{\theta}}_k)

로 갱신한다.

8.2 적응 제어와의 연계

RLS에 의한 온라인 파라미터 추정은 적응 제어(adaptive control)의 핵심 구성 요소이다. 추정된 파라미터를 제어기에 반영하면 모델 불확실성에 강인한 제어 성능을 달성할 수 있다.

8.3 마찰 모델 추정

관절 마찰은 쿨롱 마찰과 점성 마찰의 합으로 모델링하는 경우가 많다.

\tau_f = f_c \text{sgn}(\dot{q}) + f_v \dot{q}

이를 \boldsymbol{\theta}_f = [f_c, f_v]^\top에 대한 선형 모델로 표현하고, RLS로 온라인 추정하면 로봇의 운용 시간에 따른 마찰 특성 변화를 추적할 수 있다.


9. 계산 복잡도

방법시간 복잡도 (관측당)공간 복잡도
배치 최소 제곱법O(kn^2)O(kn)
RLSO(n^2)O(n^2)

RLS는 관측 횟수 k에 무관하게 일정한 계산량으로 갱신할 수 있어 실시간 응용에 적합하다.


10. 수치적 안정성

RLS의 공분산 갱신에서 \mathbf{P}_{k+1} = (\mathbf{I} - \mathbf{K}_{k+1} \mathbf{a}_{k+1}^\top) \mathbf{P}_k는 수치 오차의 누적에 의해 \mathbf{P}의 대칭성과 양의 정부호 성질이 손실될 수 있다. 이를 방지하기 위해 Joseph 안정화 형태를 사용한다.

\mathbf{P}_{k+1} = (\mathbf{I} - \mathbf{K}_{k+1} \mathbf{a}_{k+1}^\top) \mathbf{P}_k (\mathbf{I} - \mathbf{K}_{k+1} \mathbf{a}_{k+1}^\top)^\top + \mathbf{K}_{k+1} \sigma^2 \mathbf{K}_{k+1}^\top

여기서 \sigma^2는 관측 잡음의 분산이다. 이 형태는 수치적으로 대칭성과 양의 정부호 성질을 보장한다.


11. 참고 문헌

  • Ljung, L. (1999). System Identification: Theory for the User (2nd ed.). Prentice Hall.
  • Åström, K. J., & Wittenmark, B. (2008). Adaptive Control (2nd ed.). Dover Publications.
  • Slotine, J.-J. E., & Li, W. (1991). Applied Nonlinear Control. Prentice Hall.
  • Siciliano, B., Sciavicco, L., Villani, L., & Oriolo, G. (2009). Robotics: Modelling, Planning and Control. Springer.
  • Haykin, S. (2002). Adaptive Filter Theory (4th ed.). Prentice Hall.

v 0.1