이번에는 Convex Quadratic Programming(CQP)에 대해서 다뤄보려고 한다. CQP란 2차 다항식의 목적함수에 선형 제약조건을 갖는 최적화 문제를 말한다. 먼저 어떻게 생긴 문제인지 살펴보도록 하자. 모든 등식은 부등식 두 개, 예를 들어 을 와 ,로 쪼갤 수 있으므로 부등식 제약조건만을 가진  CQP에 대해서 다루도록 하겠다. 꼭 두개로 쪼개지 않더라도 제약조건을 표현하는 행렬의 Nullspace를 이용해서 부등식에 껴 넣을 수 있다.

표준형은 다음과 같다.

여기서 는 의 행렬이고 , , 의 벡터들이다.

이렇게 표현해 놓으니까 뭔지 잘 와닿지 않을 것 같다. 좀더 쉬운 예를 들어보도록 하겠다. 다음과 같은 문제를 한번 보자.

이 문제를 한번 초딩처럼 모눈종이에 그려봐? ^^

퍼런 부분이 제약 조건 네개로 이루어진 영역이고 동그라미 들은 목적함수의 등고선(contour)이다. 감좀 잡았나?  이거 고등학교 3학년(?) 때 하던 선형 계획의 확장 판이다. ^^ 수능에 자주 등장하는 그 문제.

그럼 이 문제를 어떻게 풀 것인가? 어떻게 푸는지 파해쳐 주려고 쓰는 것 아닌가? 우선 원래 문제의 Lagrangian을 구한다. 다음과 같다.

해가 있는 곳은 다음과 같다.

1. Lagrangian의 Gradient가 이 되는 곳, (이면서)

2. 제약조건과 의 Complementarity가 성립되는 곳, 즉 for , , , and .

즉, 다음이 성립되는 곳에 해가 있다.

이를 KKT (Karush-Kuhn-Tucker) 조건이라 한다. 여기서 는 사실상 를 등호 조건으로 변경하기 위해 도입한 변수라 slack 변수라 부른다.

KKT condition의 Geometry 해석은 다음 그림과 같다.

(그림)

이 KKT 조건 문제를 해결하면 최적화 문제이 답을 찾게 되는데 원래 문제의 답만 찾는게 아니라 다음과 같은 Dual 문제의 답도 찾게 된다.

초보들이 보통 개념없이 프로그램 혹은 엑셀 쉬트를 작성하면서 행렬 연산을 하다보면 극악 무도한 행동을 좀 한다 (내가 자주 목격했다!!!). 그게 뭔고 하니 Nonsingular인  Linear System 를 푸는데

를 한다는 것이다!!!

“도대체 이게 뭐가 문제지?”하고 생각했다면 당신 x잡고 반성해라. 혹시 당신 수학자거나 과학자 혹은 퀀트라면 이거 정말 심각하다. 이게 얼마나 극악 무도한 짓인지 알려주겠다.

일단 정상적으로 리니어 시스템을 풀려면 어떻게 하는지 알아둘 필요가 있다. 물론 아까도 얘기했듯이 Nonsingular인 경우다. Nonsingular가 무슨 말이냐면 매트릭스 가 Full rank이다는 말이고 이 말은 다시 의 역행렬이 존재한다는 얘기다. 그럼 도대체 어떻게 푸는 것이 좋은 건지 빨리 알고 싶다고? 행렬 분해를 해서 풀어야 한다. A가 Symmetric Postiive Definite이 아니라면 LU 분해를 사용하게 된다. LU 분해가 뭔고 하니
가 되는 Lower triangular matrix 과 Upper triangular matrix 를 찾는 것이다. 왜 이런 짓을 하냐고? 이런 행렬을 찾아내면 를 찾아내기가 쉬워서 그렇다. 자세한 과정은 수치해석 교재나 Wikipedia를 참조하자. 아무튼 이렇게 찾아냈으면 forward-substitution 즉

를 풀어 를 찾고 또 backward-substitution 즉

를 풀어 를 찾아낸다. 아 참 귀찮다고?
혹시 MATLAB을 사용한다면

를 하면 바로 위 과정을 알아서 해준다.

이제 역행렬을 구해서 푸는 게 뭐가 문젠지 알아보자. 첫째는 계산 시간 문제다. 역행렬을 구하기 위해 어떤 작업을 하는지 알아둘 필요가 있다.
역행렬을 찾아내는 것은 위에 매틀랩식 식으로 나타내면

인데 여기서 는 Identity 행렬로 행렬의 대각은 1로 채워져 있고 나머지 원소들은 0인 행렬이다. 이 문제는

를 에 대해 푸는 것과 같은데 와 는 각각 행렬 와 의 번째 열이다.

자 왜 시간이 오래 걸리는지 알겠지?

시간만 문제가 되는게 아니다. 계산 오류도 커진다. 특히 Matrix 가 singular에 가까우면 더 심각해진다. 사실 이 때는 LU분해도 문제긴 하다. 그래도 LU분해가 좀 낫긴 하다. 이 경우 를 구한 다음에 를 해보면 와 상당히 동떨어진 결과를 내 놓는다.

이 때는 Rank revealing QR 분해나 SVD 분해를 사용해서 최대한 근사치를 구하는 방식을 쓴다. 가 full rank가 아닐 때와 같은 방식으로 푼다는 얘기다. 의 rank가 거의 이라고 하면 의 레인지(Range)가 차원 공간에서 거의 차원의 subspace에 걸친다는 얘긴다. 이 때는 A의 Rank revealing QR알고리즘을 써서 QR 팩터를 구하고 Q(:,1:r)와 R(1:r,:)를 곱해 A의 근사 매트릭스 B를 만든다음 를 최소로 하는 를 찾으면 된다.

그런데 이런 모든 것을 MATLAB의 \ 연산자가 알아서 처리해준다. 그러니까 MATLAB 사용자라면 x = inv(A)*b대신에 x = A\b를 써라.

Linear algebra package  별로 지원하는 포맷과 기능을 정리해 놨다. 공짜로 구할 수 있는 코드도 상당히 많으니 직접 만들지 말고 갖다 쓰는게 더 좋을 듯. 항상 말하지만 행렬 계산 라이브러리 전문가들이 만들어 놓을 것을 사용하는게 직접 아무 생각 없이 만드는 것보다 좋다.

http://www.netlib.org/utk/people/JackDongarra/la-sw.html

그리고 다양한 행렬 계산 라이브러리를 윈도에서 사용하는 방법이 여기 있다.

http://www.matrixprogramming.com