7.25 헤시안의 로봇 최적화 문제 응용

1. 헤시안 행렬의 역할

헤시안 행렬(Hessian matrix)은 스칼라 함수 f: \mathbb{R}^n \to \mathbb{R}의 이계 편도함수로 구성된 대칭 행렬이다.

\mathbf{H}(\mathbf{x}) = \nabla^2 f(\mathbf{x}) = \begin{bmatrix} \dfrac{\partial^2 f}{\partial x_1^2} & \dfrac{\partial^2 f}{\partial x_1 \partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_1 \partial x_n} \\ \dfrac{\partial^2 f}{\partial x_2 \partial x_1} & \dfrac{\partial^2 f}{\partial x_2^2} & \cdots & \dfrac{\partial^2 f}{\partial x_2 \partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{\partial^2 f}{\partial x_n \partial x_1} & \dfrac{\partial^2 f}{\partial x_n \partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_n^2} \end{bmatrix}

fC^2 급(이계 연속 편도함수를 갖는) 함수이면, 슈바르츠 정리(Schwarz’s theorem)에 의해 \dfrac{\partial^2 f}{\partial x_i \partial x_j} = \dfrac{\partial^2 f}{\partial x_j \partial x_i}가 성립하므로 헤시안은 대칭 행렬이다.

로봇공학의 최적화 문제에서 헤시안은 비용 함수의 국소적 곡률(curvature) 정보를 제공하며, 이를 활용하여 최적화 알고리즘의 수렴 속도와 안정성을 개선할 수 있다.

2. 뉴턴법과 헤시안

2.1 뉴턴법의 유도

비용 함수 f(\mathbf{x})를 현재 점 \mathbf{x}_k 근방에서 이차 테일러 전개하면 다음과 같다.

f(\mathbf{x}_k + \boldsymbol{\delta}) \approx f(\mathbf{x}_k) + \nabla f(\mathbf{x}_k)^\top \boldsymbol{\delta} + \frac{1}{2} \boldsymbol{\delta}^\top \mathbf{H}(\mathbf{x}_k) \boldsymbol{\delta}

이 이차 근사 함수를 \boldsymbol{\delta}에 대해 최소화하면 다음을 얻는다.

\nabla f(\mathbf{x}_k) + \mathbf{H}(\mathbf{x}_k) \boldsymbol{\delta} = \mathbf{0}

따라서 뉴턴 스텝(Newton step)은 다음과 같다.

\boldsymbol{\delta}_k = -\mathbf{H}(\mathbf{x}_k)^{-1} \nabla f(\mathbf{x}_k)

갱신 규칙은 다음과 같다.

\mathbf{x}_{k+1} = \mathbf{x}_k - \mathbf{H}(\mathbf{x}_k)^{-1} \nabla f(\mathbf{x}_k)

2.2 수렴 특성

뉴턴법은 최적해 근방에서 이차 수렴(quadratic convergence) 특성을 가진다. 헤시안이 양정치이고 리프시츠 연속이면 다음이 성립한다.

\lVert \mathbf{x}_{k+1} - \mathbf{x}^* \rVert \leq C \lVert \mathbf{x}_k - \mathbf{x}^* \rVert^2

여기서 C > 0는 상수이고 \mathbf{x}^*는 극소점이다. 이는 경사하강법의 선형 수렴에 비해 매우 빠른 수렴 속도이다.

2.3 뉴턴법의 문제점과 해결

헤시안이 양정치가 아닌 경우 뉴턴 스텝이 하강 방향(descent direction)을 보장하지 못한다. 이를 해결하기 위한 방법은 다음과 같다.

수정 뉴턴법(Modified Newton’s Method):

\boldsymbol{\delta}_k = -(\mathbf{H}_k + \mu_k \mathbf{I})^{-1} \nabla f(\mathbf{x}_k), \quad \mu_k \geq 0

여기서 \mu_k\mathbf{H}_k + \mu_k \mathbf{I} \succ 0이 되도록 선택한다. 구체적으로, \lambda_{\min}이 헤시안의 최소 고유값이면 다음 조건을 만족해야 한다.

\mu_k > \max(0, -\lambda_{\min}(\mathbf{H}_k))

3. 준뉴턴법

헤시안의 직접 계산은 O(n^2)개의 이계 도함수를 요구하므로 고차원 문제에서 계산 비용이 높다. 준뉴턴법(quasi-Newton method)은 그래디언트 정보만을 이용하여 헤시안의 근사 행렬 \mathbf{B}_k \approx \mathbf{H}_k를 점진적으로 갱신한다.

3.1 BFGS 방법

BFGS(Broyden-Fletcher-Goldfarb-Shanno) 갱신 공식은 다음과 같다.

\mathbf{B}_{k+1} = \mathbf{B}_k - \frac{\mathbf{B}_k \mathbf{s}_k \mathbf{s}_k^\top \mathbf{B}_k}{\mathbf{s}_k^\top \mathbf{B}_k \mathbf{s}_k} + \frac{\mathbf{y}_k \mathbf{y}_k^\top}{\mathbf{y}_k^\top \mathbf{s}_k}

여기서 다음과 같이 정의한다.

\mathbf{s}_k = \mathbf{x}_{k+1} - \mathbf{x}_k, \quad \mathbf{y}_k = \nabla f(\mathbf{x}_{k+1}) - \nabla f(\mathbf{x}_k)

시컨트 조건(secant condition) \mathbf{B}_{k+1} \mathbf{s}_k = \mathbf{y}_k가 만족되며, \mathbf{y}_k^\top \mathbf{s}_k > 0이면 \mathbf{B}_{k+1}의 양정치성이 보존된다.

3.2 L-BFGS 방법

대규모 문제에서는 \mathbf{B}_k를 명시적으로 저장하지 않고, 최근 m개의 벡터 쌍 \{(\mathbf{s}_i, \mathbf{y}_i)\}_{i=k-m}^{k-1}만 저장하여 헤시안-벡터 곱을 암묵적으로 계산하는 L-BFGS(Limited-memory BFGS) 방법을 사용한다. 메모리 사용량은 O(mn)으로 O(n^2)에 비해 대폭 절감된다.

4. 역기구학 최적화

로봇 매니퓰레이터의 역기구학 문제를 최소화 문제로 정식화하면 다음과 같다.

\min_{\mathbf{q}} f(\mathbf{q}) = \frac{1}{2} \lVert \mathbf{e}(\mathbf{q}) \rVert^2 = \frac{1}{2} \mathbf{e}(\mathbf{q})^\top \mathbf{e}(\mathbf{q})

여기서 \mathbf{e}(\mathbf{q}) = \mathbf{x}_d - \mathbf{k}(\mathbf{q})는 작업 공간 오차이고, \mathbf{k}(\mathbf{q})는 순기구학 함수이다.

4.1 그래디언트와 헤시안

그래디언트는 다음과 같다.

\nabla f(\mathbf{q}) = -\mathbf{J}(\mathbf{q})^\top \mathbf{e}(\mathbf{q})

여기서 \mathbf{J}(\mathbf{q}) = \dfrac{\partial \mathbf{k}}{\partial \mathbf{q}}는 자코비안 행렬이다. 헤시안은 다음과 같다.

\mathbf{H}(\mathbf{q}) = \mathbf{J}(\mathbf{q})^\top \mathbf{J}(\mathbf{q}) - \sum_{i=1}^{m} e_i(\mathbf{q}) \frac{\partial^2 k_i}{\partial \mathbf{q} \partial \mathbf{q}^\top}

여기서 m은 작업 공간의 차원이다. 첫째 항 \mathbf{J}^\top \mathbf{J}는 가우스-뉴턴 근사(Gauss-Newton approximation)에 해당하며, 잔차가 작을 때 둘째 항을 무시할 수 있다.

4.2 가우스-뉴턴법

가우스-뉴턴법은 헤시안의 근사 \mathbf{H} \approx \mathbf{J}^\top \mathbf{J}를 사용한다. 갱신 스텝은 다음과 같다.

\boldsymbol{\delta}_k = -(\mathbf{J}_k^\top \mathbf{J}_k)^{-1} \mathbf{J}_k^\top \mathbf{e}_k

이는 자코비안의 의사역행렬(pseudoinverse)을 이용한 것으로, \mathbf{J}_k가 열 최대 계수(full column rank)이면 \mathbf{J}_k^\top \mathbf{J}_k는 양정치이므로 하강 방향이 보장된다.

4.3 레벤버그-마쿼트법

레벤버그-마쿼트법(Levenberg-Marquardt method)은 가우스-뉴턴법에 감쇠 인자를 추가한 것이다.

\boldsymbol{\delta}_k = -(\mathbf{J}_k^\top \mathbf{J}_k + \lambda_k \mathbf{I})^{-1} \mathbf{J}_k^\top \mathbf{e}_k

감쇠 인자 \lambda_k > 0가 크면 경사하강법에 가깝고, 작으면 가우스-뉴턴법에 가깝다. 이 방법은 자코비안이 특이(singular)에 가까운 경우에도 안정적으로 동작한다.

5. 궤적 최적화에서의 헤시안

로봇 궤적 최적화 문제에서 비용 함수는 이산화된 궤적 \mathbf{q}_0, \mathbf{q}_1, \ldots, \mathbf{q}_N에 대해 다음과 같이 정의한다.

J = \sum_{k=0}^{N-1} c(\mathbf{q}_k, \mathbf{q}_{k+1})

여기서 c는 스테이지 비용 함수이다. 결정 변수 벡터를 \boldsymbol{\xi} = [\mathbf{q}_1^\top, \mathbf{q}_2^\top, \ldots, \mathbf{q}_{N-1}^\top]^\top로 구성하면, 헤시안 행렬은 블록 삼대각(block tridiagonal) 구조를 가진다.

\mathbf{H} = \begin{bmatrix} \mathbf{D}_1 & \mathbf{E}_1 & & \\ \mathbf{E}_1^\top & \mathbf{D}_2 & \mathbf{E}_2 & \\ & \ddots & \ddots & \ddots \\ & & \mathbf{E}_{N-2}^\top & \mathbf{D}_{N-1} \end{bmatrix}

이 희소 구조를 활용하면 뉴턴 스텝을 O(Nn^3)의 계산량으로 효율적으로 구할 수 있다. 여기서 n은 관절 공간의 차원이다.

6. 제약 조건이 있는 최적화

6.1 등식 제약 최적화

등식 제약 \mathbf{h}(\mathbf{x}) = \mathbf{0}이 있는 최적화 문제의 라그랑주 함수는 다음과 같다.

\mathcal{L}(\mathbf{x}, \boldsymbol{\lambda}) = f(\mathbf{x}) + \boldsymbol{\lambda}^\top \mathbf{h}(\mathbf{x})

라그랑주 함수의 헤시안은 다음과 같다.

\nabla^2_{\mathbf{x}\mathbf{x}} \mathcal{L} = \nabla^2 f(\mathbf{x}) + \sum_{i=1}^{p} \lambda_i \nabla^2 h_i(\mathbf{x})

KKT 시스템의 뉴턴 스텝은 다음 선형계를 풀어 구한다.

\begin{bmatrix} \nabla^2_{\mathbf{x}\mathbf{x}} \mathcal{L} & \mathbf{A}^\top \\ \mathbf{A} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \boldsymbol{\delta}_x \\ \boldsymbol{\delta}_\lambda \end{bmatrix} = -\begin{bmatrix} \nabla_{\mathbf{x}} \mathcal{L} \\ \mathbf{h} \end{bmatrix}

여기서 \mathbf{A} = \dfrac{\partial \mathbf{h}}{\partial \mathbf{x}}는 제약 조건의 자코비안이다.

6.2 부등식 제약과 내점법

부등식 제약 g_j(\mathbf{x}) \leq 0이 있는 경우, 내점법(interior point method)은 장벽 함수(barrier function)를 사용한다.

f_\mu(\mathbf{x}) = f(\mathbf{x}) - \mu \sum_{j=1}^{q} \ln(-g_j(\mathbf{x}))

장벽 함수의 헤시안은 다음과 같다.

\nabla^2 f_\mu = \nabla^2 f + \mu \sum_{j=1}^{q} \left( \frac{\nabla g_j \nabla g_j^\top}{g_j^2} - \frac{\nabla^2 g_j}{g_j} \right)

관절 한계(joint limit), 충돌 회피, 토크 제한 등 로봇의 물리적 제약을 부등식 제약으로 표현하고 내점법을 적용하면, 헤시안 정보를 활용한 효율적인 제약 최적화가 가능하다.

7. 수치적 헤시안 계산

해석적 헤시안 계산이 어려운 경우, 유한 차분(finite difference)을 이용하여 수치적으로 근사할 수 있다.

전방 차분:

\frac{\partial^2 f}{\partial x_i \partial x_j} \approx \frac{f(\mathbf{x} + h\mathbf{e}_i + h\mathbf{e}_j) - f(\mathbf{x} + h\mathbf{e}_i) - f(\mathbf{x} + h\mathbf{e}_j) + f(\mathbf{x})}{h^2}

중심 차분:

\frac{\partial^2 f}{\partial x_i \partial x_j} \approx \frac{f(\mathbf{x} + h\mathbf{e}_i + h\mathbf{e}_j) - f(\mathbf{x} + h\mathbf{e}_i - h\mathbf{e}_j) - f(\mathbf{x} - h\mathbf{e}_i + h\mathbf{e}_j) + f(\mathbf{x} - h\mathbf{e}_i - h\mathbf{e}_j)}{4h^2}

전방 차분은 O(h)의 오차를 가지며, 중심 차분은 O(h^2)의 오차를 가진다. n차원 문제에서 전체 헤시안을 계산하려면 O(n^2)회의 함수 평가가 필요하다.

자동 미분(automatic differentiation)을 사용하면 수치적 오차 없이 정확한 헤시안을 계산할 수 있으며, 특히 역방향 모드 자동 미분을 두 번 적용하면 효율적으로 헤시안-벡터 곱 \mathbf{H}\mathbf{v}를 구할 수 있다.


참고 문헌

  • Nocedal, J., & Wright, S. J. (2006). Numerical Optimization. 2nd Edition. Springer.
  • Siciliano, B., Sciavicco, L., Villani, L., & Oriolo, G. (2009). Robotics: Modelling, Planning and Control. Springer.
  • Boyd, S., & Vandenberghe, L. (2004). Convex Optimization. Cambridge University Press.
  • Ratliff, N., Zucker, M., Bagnell, J. A., & Srinivasa, S. (2009). CHOMP: Gradient Optimization Techniques for Efficient Motion Planning. IEEE International Conference on Robotics and Automation.
  • Gill, P. E., Murray, W., & Wright, M. H. (1981). Practical Optimization. Academic Press.

v 0.1