이 글은 고려대학교 통계학과 신승준 교수님의 2018년 1학기 통계 계산 방법론(대학원) 강의를 바탕으로 작성되었습니다.
참고할 만한 링크들
첫번째 블로그, 핸즈온머신러닝, 패터인식(오일석) 참고하여 쓰인 글
Optimization with inequality constraints
Lagrangian method를 통해 equality constraint가 있는 minimization/maximization 문제를 풀 수 있다고 했다. 여기서는 inequality constraint를 포함하는 최적화 문제를 다룬다.
이를 일반화 시켜서 표현하면 다음과 같다.
\[\underset{x} {min} f(x) \quad s.t \quad h(x) \le 0, \;\; l(x)=0\]
- constraint가 여러 개 있어도 상관 없고, Lagrangian method와 같은 방식으로 풀 수 있다.
Primal problem
여기서 Lagrangian primal function이라는 것을 정의할 수 있는데 이는 다음과 같다.
\[L(x, u, v) = f(x) + u\;h(x) + v\;l(x)\]
- inequality constraint에 대해 \(\quad u\ge0\)라는 조건이 필요하다.
- x에 대한 함수로 볼 수 있고, 이 문제를 풀 수 있는 x를 주는 적절한 u와 v가 존재한다.
- 기본적으로 x에 대한 optimization 문제로 볼 수 있고, 이것을
Primal problem이라고 한다. - 이것을 u와 v에 대한 optimization 문제로 바꿔줄 수 있는데, 이것을
Dual problem이라고 한다. - Primal problem과 Dual problem은 적절한 조건을 만족하면 같은 solution을 준다.
- Primal problem이 쉬울 때는 굳이 Dual problem으로 넘어갈 필요가 없지만, Primal problem을 푸는 것이 어렵다면 Dual problem으로 바꾼 후 푸는 것이 편한 경우가 있다.
Dual problem을 정의하기 전에 먼저 Primal solution \(x^{*}\)를 정의할 필요가 있다.
\[x^{*}=\underset{x \in C} {argmin} \; f(x)\]
C는 feasible region(해의 영역)이라고 하고, 다음과 같이 정의된다.
\[C=\{ x; h(x) \le 0 \;\; \& \;\; l(x)=0 \}\]
이 조건 하에서는 아래의 두 가지 부등식이 성립한다.
(1)
\[f(x^{*}) \ge \underset{x \in C} {min} L(x, u, v)\]
\(L(x, u, v) = f(x) + u\;h(x) + v\;l(x)\)이므로 이를 minimize 하면 \(f(x^{*}) + u\;h(x) + v\;l(x)\)가 되는데, C안에서 \(h(x) \le 0\)이고 \(l(x)=0\)이므로 \(u\;h(x) + v\;l(x)\)는 non-positive가 된다. 따라서 위의 부등식이 성립한다.
(2)
\[ \underset{x \in C} {min} f(x) \ge \underset{x}{min} L(x,u,v)\]
- global region에서 minimize한 solution이 C에서 minimize한 solution보다 작거나 같다는 의미이다.
- 이런 경우에 \(\underset{x}{min} L(x,u,v)\)를 \(g(u,v)\)로 정의하고, 이를 Dual function이라고 한다.
Dual problem
이제 Dual problem을 아래와 같이 정의할 수 있다.
\[(u^{*}, v^{*})=\underset{u, v}{argmax} \; g(u,v) \quad s.t \;\; u \ge 0\]
- \(u^{*}, v^{*}\): Dual solution
\[f(x^{*}) \ge \underset{x \in C} {min} L(x, u, v)\]
\(f(x^{*}) \ge \underset{x \in C} {min} L(x, u, v)\)과
\(\underset{x \in C} {min} f(x) \ge \underset{x}{min} L(x,u,v)\)에서
\(f(x^{*}) \ge \underset{x}{min} L(x,u,v)\)를 얻을 수 있고, \(\underset{x}{min} L(x,u,v)\)를 \(g(u,v)\)로 정의했으므로, \(f(x^{*}) \ge g(u,v)\)가 성립한다.
따라서 u,v를 아무리 maximize해도 \(f(x^{*})\)를 넘을 수 없고, 이를 아래와 같이 표현할 수 있다.
\[f(x^{*}) \ge g(u^{*}, v^{*})\]
Primal problem은 \(L(x, u, v)\)를 x에 대해 minimize하는 문제였는데, Dual problem은 \(g(u, v)\)를 u, v에 대해 maximize하는 문제로 바뀌었다.
어떤 조건을 만족하면 \(f(x^{*}) = g(u^{*}, v^{*})\)이 성립하는데, 이 조건을 Slater's condition이라고 한다.
Duality Gap
\(f(x^{*}) - g(u^{*}, v^{*})\)를 Duality Gap으로 정의하는데,
\(DG \ge 0\)이면 weak duality라 하고,
\(DG = 0\)이면 strong duality라 한다.
- 통계학에서 다루는 대부분의 dual problem은 strong duality가 성립한다. (ex. SVM)
In mathematical optimization, the Karush<u+2013>Kuhn<u+2013>Tucker (KKT) conditions</u+2013></u+2013>, also known as the Kuhn<u+2013>Tucker conditions</u+2013>, are first-order necessary conditions for a solution in nonlinear programming to be optimal, provided that some regularity conditions are satisfied.
Allowing inequality constraints, the KKT approach to nonlinear programming generalizes the method of Lagrange multipliers, which allows only equality constraints.
The system of equations and inequalities corresponding to the KKT conditions is usually not solved directly, except in the few special cases where a closed-form solution can be derived analytically.
In general, many optimization algorithms can be interpreted as methods for numerically solving the KKT system of equations and inequalities.
The KKT conditions were originally named after Harold W. Kuhn, and Albert W. Tucker, who first published the conditions in 1951. Later scholars discovered that the necessary conditions for this problem had been stated by William Karush in his master’s thesis in 1939.
KKT condition은 x, u, v를 계산해내고 싶을 때 만족해야 하는 필요충분 조건이다. 단, strong duality를 만족하는 상황에서만 필요충분 조건이 성립한다.
\[KKT\; condition \Leftrightarrow \;\; solution\]
(i) \(\frac{\partial L}{\partial x} = 0\) - Lagrange method와 같음, Stationary condition(참고: stationary point)
(ii) \(uh(x) = 0\), \(vl(x) = 0\)
(iii) \(h(x) \le 0\)
(iv) \(l(x) = 0\)
- 이 4가지 조건을 만족하는 x, u, v는 이 문제의 primal/dual solution이고, 반대로 primal/dual solution이라고 이 조건을 만족해야 한다.
- (iii)과 (iv)는 feasible region이면 성립하므로 feasible region에서 (i)과 (ii)를 만족하는 x를 찾으면 그게 \(x^{(*)}\)
- 경제학에서 많이 쓰이고, 머신러닝에서는 SVM에 쓰인다.
Support Vector Machine(SVM)
Optimal Separating Hyperplane
\((x_i, y_i)\), \(i=1,\cdots, n\),
\(x_i \in R^p\) => p-dimensional space
\(y_i= \{-1, 1\}\) => binary outcome
이라 하자.
p-1 차원 hyperplane \(f(x)=\beta_o + \beta^Tx\)를 찾는 것이 목적이고, \(\hat{y_i}=sign(f(x))\)로 y를 예측한다.
이거 맞나? -->
- SVM은 두 집단을 가장 잘 분리시키는 hyperplane(2차원에서는 선)을 찾는 알고리즘이다. 이것을
optimal Separating Hyperplane이라고 한다. - 직선과 집단(data cloud) 사이의 거리를 재서 그 거리를 가장 크게 만들어주는 직선을 찾는 것.
- 직선에서 각 집단에 속한 가장 가까운 점에 내린 수선의 길이를 거리라고 정의한다.
- 직선 양쪽으로의 거리 합을 최대화하는데, 두 거리의 합이 일정하므로 기울기가 변하지 않으면 직선을 평행이동 시켜도 거리는 일정하다. 두 개의 거리가 같은 값(M)을 갖도록 정의하면 하나의 직선을 특정할 수 있다.
이제 이 문제를 optimization problem으로 표현해보자. 그러려면 M을 \(\beta\)와 \(\beta_0\)의 함수로 표현해야 한다.
- M을 maximize하는 \(\beta\)와 \(\beta_0\)를 찾아야 한다.
- 원점에서 직선에 내린 수선을 \(\beta\)라고 할 수 있고, \(\beta\)와 \(x_1-x_0\)는 직교하므로 내적은 0이다.
M을 \(\beta\)와 \(\beta_0\)의 함수로 표현하기 위해서 두 벡터 사이의 각도 \(\theta\)의 정의와 직각 삼각형에서 \(cos\theta\)의 정의를 이용한다.
두 벡터 사이의 각도 \(\theta\)의 정의에 의해 \(cos\theta = \frac{<\beta,\;\; x-x_0>}{||\beta|| \;||x-x_0||} \cdots (a)\)이고,
직각삼각형에서 \(cos\theta\)의 정의에 의해 \(cos\theta= \frac{M}{||x-x_0||} \cdots (b)\)이다.
(b)를 \(M=cos\theta ||x-x_0||\)로 표현할 수 있고, 여기에 (a)을 대입하면
\[M=cos\theta \; ||x-x_0|| = \frac{1}{||\beta||}(\beta^Tx-\beta^Tx_0)\]
인데 \(\beta^Tx_0=-\beta_0\)이므로
\[M=\frac{1}{||\beta||}(\beta_0 + \beta^Tx)\]
이 된다. 그런데 여기서 M은 거리이기 때문에 양수여야 한다. M은 직선의 위쪽에 있는 점들에 대해서는 양수값을 갖고, 직선의 아래쪽에 있는 점들에 대해서는 음수 값을 갖는다. 이를 방향에 관계없이 표현하기 위해
\[M=\frac{y}{||\beta||}(\beta_0 + \beta^Tx)\]
로 표현할 수 있는데, 여기서 y는 response로 점들이 직선 위쪽에 있을 때는 1, 직선 아래쪽에 있을 때는 -1의 값을 갖는다. 그래서 y를 0, 1이 아니라 -1, 1로 정의하는 것이다. 로지스틱 회귀의 경우 Bernoulli distribution에서 온 것이기 때문에 0, 1이 자연스럽지만 SVM은 distribution에서 나온 것이 아니라 위에서 설명한 geometry에서 나온 것이기 때문에 -1, 1로 코딩하는 것이 더 좋다.
Maximal Margin Classifier - linearly separable case
- M이 정해지면 직선과 나머지 점들 사이의 거리는 M보다 커야 한다. 이를 식으로 표현하면 아래와 같다.
\[\underset{\beta_0, \beta_1} {max} M \quad s.t \quad \frac{y_i}{||\beta||}(\beta_0+\beta^Tx_i)\ge M\]
여기서 \(\beta\)를 unique하게 표현하려면 norm을 어떤 값으로 정해줘야 하는데 그 값을 \(\frac{1}{M}\)로 하면 풀어야 하는 식을 \(\beta\)에 대한 식으로 바꿀 수 있다.
\[||\beta||=\frac{1}{M} \Leftrightarrow M=\frac{1}{||\beta||}\]
M을 maximize하는 것은 \(||\beta||\)를 minimize하는 것과 같으므로
\[\underset{\beta_0, \beta}{min}\frac{1}{2}||\beta||^2 \quad s.t \quad y_i(\beta_0 + \beta^Tx_i)\ge1\]
\[\Leftrightarrow\]
\[\underset{\beta_0, \beta}{min}\frac{1}{2}\beta^T\beta \quad s.t \quad y_i(\beta_0 + \beta^Tx_i)\ge1\]
와 같이 표현할 수 있다. 여기서 \(||\beta||\)는 norm이므로 \(||\beta||^2\)로 표현할 수 있고, \(1 \over 2\)를 곱한 것은 미분을 할 때 계산을 편하게 하기 위함이다.
- \(\quad y_i(\beta_0 + \beta^Tx_i)\ge1\) 이라는 constraint 때문에 y를 -1과 1로 사용하는 것이다.
Soft Margin Classifier(linear SVM) - linearly nonseparable case
Figure 12.7 in Element of Statistical Learning
Maximal margin classifier는 점들이 완전히 선형적으로 분리되는 경우에만 적용할 수 있다. 위의 그림과 같이 점들이 완전히 선형적으로 분리되지 않는 경우에는 Maximal margin classifier의 constraint를 조금 느슨하게 해줘야 한다.
\[\underset{\beta_0, \beta}{min}\frac{1}{2}\beta^T\beta + C \overset{n}{\underset{i=1}{\Sigma}} \xi_i \quad s.t \quad y_i(\beta_0 + \beta^Tx_i) \ge 1 -\xi_i, \quad \xi_i \ge 0\] Soft margin classifier는 몇몇 점들이 hyperplane 반대쪽으로 넘어가는 것을 허용하는데 너무 많은 점들이 넘어가는 것을 막기 위하여 penalty로 \(C \overset{n}{\underset{i=1}{\Sigma}} \xi_i\)를 준다.
\(\xi_i\)를 slack variable이라고 하는데
- i번째 관측치가 margin 안쪽 있으면 correct side of margin에 있다고 하고, 이때 \(\xi_i=0\)이 된다.
- i번째 관측치가 margin 바깥쪽에 있지만 hyperplane 안쪽에 있으면 wrong side of margin에 있다고 하고, 이때 \(\xi_i>0\)이 된다.
- i번째 관측치가 hyperplane 바깥쪽에 있으면 wrong side of margin에 있다고 하고, 이때 \(\xi_i>1\)이 된다.
C(Cost)는 nonnegative tuning parameter로 관측치들이 hyperplane 반대쪽으로 넘어가는 것에 얼만큼 penalty를 줄 것인가를 결정한다. 일반적으로 cross-validation을 통해 결정한다.
- When C is small, we seek narrow margins that are rarely violated; this amounts to a classifier that is highly fit to the data, which may have low bias but high variance.
- On the other hand, when C is larger, the margin is wider and we allow more violations to it; this amounts to fitting the data less hard and obtaining a classifier that is potentially more biased but may have lower variance.
이제 \(\underset{\beta_0, \beta}{min}\frac{1}{2}\beta^T\beta + C \overset{n}{\underset{i=1}{\Sigma}} \xi_i \quad s.t \quad y_i(\beta_0 + \beta^Tx_i) \ge 1 -\xi_i, \quad \xi_i \ge 0\)를 풀어야 하는데 여기서 KKT condition을 포함한 앞에서 언급한 개념들이 사용된다. Primal variable은 \(\beta_0\)(1개), \(\beta\)(p개), \(\xi\)(n개)이고, 2n개의 constriant가 존재한다.
Lagrangian primal function은 다음과 같다.
\[L_P(\beta_0, \beta, \xi_i)= \frac{1}{2}\beta^T\beta + C \overset{n}{\underset{i=1}{\Sigma}} \xi_i +\overset{n}{\underset{i=1}{\Sigma}}\alpha_i (1-y_i(\beta_0+\beta^Tx_i)-\xi_i)-\overset{n}{\underset{i=1}{\Sigma}}\mu_i\xi_i \quad \cdots(1)\]
\(\alpha_i \ge 0, \quad \mu_i \ge 0\)
\(L_P(\beta_0, \beta, \xi_i)\)를 minimize하는 \(\beta_0, \;\beta, \;\xi_i\)는
- \(\frac{\partial L}{\partial \beta} = \beta - \overset{n}{\underset{i=1}{\Sigma}}\alpha_iy_ix_i=0\) \(\Leftrightarrow\) \(\beta = \overset{n}{\underset{i=1}{\Sigma}}\alpha_iy_ix_i\quad \cdots(2)\)
- \(\frac{\partial L}{\partial \beta_0} = \overset{n}{\underset{i=1}{\Sigma}}\alpha_iy_i=0\quad \cdots(3)\)
- \(\frac{\partial L}{\partial \xi_i}\) \(\Leftrightarrow\) \(\alpha_i = C-\mu_i, \quad i=1,\cdots,n\quad \cdots(4)\)
를 통해 얻을 수 있고, constraint는 아래와 같다.
\(\alpha_i (1-y_i(\beta_0+\beta^Tx_i)-\xi_i)=0, \quad i=1,\cdots,n\)
\(\mu_i\xi_i=0 \quad i=1,\cdots,n\)
(1)에 (2), (3), (4)를 대입해서 풀면 Lagrangian dual function을 얻을 수 있다.
\(L_D=\frac{1}{2}\overset{n}{\underset{i=1}{\Sigma}}\overset{n}{\underset{j=1}{\Sigma}} \alpha_i y_i x_i^Tx_j y_j \alpha_j+C \overset{n}{\underset{i=1}{\Sigma}}\xi_i+\overset{n}{\underset{i=1}{\Sigma}}\alpha_i-\overset{n}{\underset{i=1}{\Sigma}}\alpha_iy_i\beta_0-\overset{n}{\underset{i=1}{\Sigma}}\alpha_iy_i\beta^Tx_i-\overset{n}{\underset{i=1}{\Sigma}}\alpha_i\xi_i-\overset{n}{\underset{i=1}{\Sigma}}\mu_i\xi_i\)
\(L_D=\frac{1}{2}\overset{n}{\underset{i=1}{\Sigma}}\overset{n}{\underset{j=1}{\Sigma}} \alpha_i y_i x_i^Tx_j y_j \alpha_j+C \overset{n}{\underset{i=1}{\Sigma}}\xi_i+\overset{n}{\underset{i=1}{\Sigma}}\alpha_i-\beta_0\overset{n}{\underset{i=1}{\Sigma}}\alpha_iy_i-\overset{n}{\underset{i=1}{\Sigma}}\overset{n}{\underset{j=1}{\Sigma}} \alpha_i y_i x_i^Tx_j y_j \alpha_j-\overset{n}{\underset{i=1}{\Sigma}}\alpha_i\xi_i-\overset{n}{\underset{i=1}{\Sigma}}\mu_i\xi_i\)
(3)에서 \(\overset{n}{\underset{i=1}{\Sigma}}\alpha_iy_i=0\)이고, (4)에서 \(\alpha_i = C-\mu_i\)이므로
\(L_D=\overset{n}{\underset{i=1}{\Sigma}}\alpha_i - \frac{1}{2}\overset{n}{\underset{i=1}{\Sigma}}\overset{n}{\underset{j=1}{\Sigma}} \alpha_i y_i x_i^Tx_j y_j \alpha_j\)
를 얻을 수 있고, primal variable이 모두 사라진 것을 확인할 수 있다.
이제 \(\alpha\)에 대해 다음 식을 풀면 된다.
\[\hat{\alpha}=\underset{\alpha}{argmax}[\overset{n}{\underset{i=1}{\Sigma}}\alpha_i-\frac{1}{2}\overset{n}{\underset{i=1}{\Sigma}}\overset{n}{\underset{j=1}{\Sigma}} \alpha_i y_i x_i^Tx_j y_j \alpha_j]\]
- 이 식을 다른 형태로 바꾸면 \(1^T\alpha-\frac{1}{2}\alpha^TK\alpha\)의 꼴이 되고, \(\alpha\)에 대해 linear constraint를 준다. 이런 형태를
Quadratic objective function이라고 하고, 이것을 푸는 것을Quadratic Programming문제라고 한다.
linear constraint는 아래와 같다.
\(C \ge \alpha_i \ge 0, \;\; i=1, \cdots, n\) \(\quad \quad\) <= \(\quad \quad\) \(\mu_i\)가 남긴 유산
\(\overset{n}{\underset{i=1}{\Sigma}}y_i\alpha_i=0\)
\(\alpha\)를 구한 후에는 \(\beta\)와 \(\beta_0\)를 구해야하는데
\(\beta\)는 \(\beta = \overset{n}{\underset{i=1}{\Sigma}}\alpha_iy_ix_i\quad \cdots(2)\)를 통해 구하면 된다.
\(\beta_0\)를 구하는 것이 조금 까다로운데
i) \(y_i(\beta_0 + \beta^Tx_i)>1\)인 경우(관측값이 support vector 반대쪽으로 넘어가지 않은 경우; 제대로 분류된 경우)
\(y_i(\beta_0 + \beta^Tx_i) \ge 1-\xi_i\)에서 \(y_i(\beta_0 + \beta^Tx_i)>1\)이 되려면 \(\xi_i=0\)이어야 한다.
그러면 KKT condition \(\alpha_i (1-y_i(\beta_0+\beta^Tx_i)-\xi_i)=0\)에서 \(\xi_i=0\)이면 \(\alpha_i\)가 0이 되어야 한다.
ii) \(y_i(\beta_0 + \beta^Tx_i)<1\)인 경우(관측값이 support vector 반대쪽으로 넘어간 경우)
\(y_i(\beta_0 + \beta^Tx_i) \ge 1-\xi_i\)에서 \(y_i(\beta_0 + \beta^Tx_i)<1\)이 되려면 \(\xi_i>0\)이어야 한다.
그러면 KKT condtion \(\mu_i\xi_i=0\)에서 \(\mu_i=0\)이어야 한다.
이때, \(\alpha_i = C-\mu_i, \quad i=1,\cdots,n\quad \cdots(4)\)에서 \(\alpha=C\)를 얻을 수 있다.
iii) \(y_i(\beta_0 + \beta^Tx_i)=1\)인 경우(관측값이 support vector)
그러면 KKT condition \(\alpha_i (1-y_i(\beta_0+\beta^Tx_i)-\xi_i)=0\)에서 \(\alpha_i (1-1-\xi_i)=0\)이므로 \(\alpha_i\xi_i=0\)을 얻을 수 있다.
또한 \(y_i(\beta_0 + \beta^Tx_i) \ge 1-\xi_i\)에서 \(1 \ge 1-\xi_i\) \(\Leftrightarrow\) \(\xi_i \ge 0\)을 얻을 수 있다.
따라서 \(\alpha_i\)나 \(\xi_i\)는 0이상의 어떤 값이 될 수 있다.
그런데, \(C \ge \alpha_i \ge 0\)라는 linear constraint가 있으므로 \(\alpha_i\)의 범위는 \(0 \le \alpha_i \le C\)이다.
\(\alpha_i\)가 이 범위에 들어오는 것이 있으면 해당되는 인덱스의 \(x_i\)와 \(y_i\)를 찾아서 \(\beta\)를 구하고,
\(y_i(\beta_0 + \beta^Tx_i)=1\)에서 양변에 \(y_i\)를 곱하면 \(y_iy_i=1\)이므로
\(y_iy_i(\beta_0 + \beta^Tx_i)=y_i\) \(\Leftrightarrow\) \((\beta_0 + \beta^Tx_i)=y_i\) \(\Leftrightarrow\) \(\beta_0=y_i-\beta^Tx_i\)
을 통해 \(\beta_0\)를 구할 수 있다.
Point
- \(\alpha\)를 구하면 \(\beta\)와 \(\beta_0\)를 구할 수 있다.
- SVM을 푼다 = Quadratic Programming 문제를 푼다.
R code
Data setting
library(quadprog)
set.seed(1)
n <- 100 # sample size
p <- 2 # predictor dimension
C <- 1 # cost for relaxation
# predictor
x <- matrix(rnorm(n*p), n, p)
# true parameters
true.beta0 <- 0
true.beta0## [1] 0true.beta <- rep(2, p)
true.beta## [1] 2 2true desision function을 만든다.
\(f(x)=\beta_0 + \beta^T x\)
fx <- true.beta0 + x %*% true.betaseparable하지 않게 하기 위해 error를 추가한다.
e <- rnorm(n)output label을 만든다.
y <- c(sign(fx + e))
y## [1] -1 1 -1 1 -1 1 1 1 1 1 1 1 1 -1 1 1 -1 1 1 1 -1 1 -1
## [24] -1 -1 1 -1 -1 -1 1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 1 -1 -1 -1 -1
## [47] 1 1 -1 -1 1 -1 -1 -1 1 1 -1 -1 -1 1 1 -1 1 1 -1 1 -1 -1 -1
## [70] 1 1 -1 1 1 -1 1 -1 1 1 1 -1 1 1 -1 1 -1 1 -1 1 -1 -1 1
## [93] 1 1 1 -1 1 -1 -1 -1plot(x[,1], x[,2], type = "n", xlab = "x1", ylab = "x2",
main = "Linear Support Vector Machine")
points(x[y == 1, 1], x[y == 1, 2], pch = 1, col = 2)
points(x[y != 1, 1], x[y != 1, 2], pch = 3, col = 4)
abline(a = -true.beta0/true.beta[2], b = -true.beta[1]/true.beta[2], lty = 2)
- plot을 그려보면 hyperplane을 넘어가는 관측값들이 존재해서 nonseparable case인 것을 알 수 있다.
- \(\beta_0 + \beta_1x_1 + \beta_2x_2=0 \Leftrightarrow x_2=-\frac{\beta_0}{\beta_2}-\frac{\beta_1}{\beta_2}x_1\)
Quadratic Programming for SVM
\[\hat{\alpha}=\underset{\alpha}{argmax}[\overset{n}{\underset{i=1}{\Sigma}}\alpha_i-\frac{1}{2}\overset{n}{\underset{i=1}{\Sigma}}\overset{n}{\underset{j=1}{\Sigma}} \alpha_i y_i x_i^Tx_j y_j \alpha_j]\]
\[1^T\alpha-\frac{1}{2}\alpha^TK\alpha\]
K <- (x %*% t(x)) * (outer(y, y)) + 1.0e-8 * diag(rep(1, n))
dim(K)## [1] 100 100- 1.0e-8 * diag(rep(1, n))는 numerically singular가 되는 것을 방지한기 위한 term이다.
d <- rep(1, n)
A <- cbind(y, diag(n), -diag(n)) # 100 x 201
b0 <- c(0, rep(0, n), rep(-C, n))
obj <- solve.QP(K, d, A, b0, meq = 1)
alpha <- obj$solution # dual solutionsolve.QP {quadprog}
Description
This routine implements the dual method of Goldfarb and Idnani (1982, 1983) for solving quadratic programming problems of the form \(min(-d^T b + \frac{1}{2} b^T D b)\) with the constraints \(A^T b \ge b_0\).
solve.QP(Dmat, dvec, Amat, bvec, meq=0, factorized=FALSE)
Arguments
Dmat: matrix appearing in the quadratic function to be minimized.
dvec: vector appearing in the quadratic function to be minimized.
Amat: matrix defining the constraints under which we want to minimize the quadratic function.
bvec: vector holding the values of \(b_0\) (defaults to zero).
meq: the first meq constraints are treated as equality constraints, all further as inequality constraints (defaults to 0).
\(min(-d^T b + \frac{1}{2} b^T D b)\) \(\Leftrightarrow\) \(1^T\alpha-\frac{1}{2}\alpha^TK\alpha\)에서
- \(d(dvec) = 1\), \(b = \alpha\), \(D(Dmat) = K\), \(b_0(bvec): b0\)
- \(A\) : \(\underset{\alpha}{argmax}[\overset{n}{\underset{i=1}{\Sigma}}\alpha_i-\frac{1}{2}\overset{n}{\underset{i=1}{\Sigma}}\overset{n}{\underset{j=1}{\Sigma}} \alpha_i y_i x_i^Tx_j y_j \alpha_j]\) 를 풀기위한 linear constraint
- \(\overset{n}{\underset{i=1}{\Sigma}}y_i\alpha_i=0\) \(\Leftrightarrow\) \(y^T\alpha=0\)
- \(0 \le \alpha \le C\) \(\quad \Leftrightarrow \quad\) \(0 \le \alpha \Rightarrow I\alpha \quad -C \le -\alpha \Rightarrow -I\alpha\)
- meq=1 => \(y^T\alpha=0\)을 만들기 위해 필요
이제 \(\beta_0\)를 recover하기 위해 support vector의 인덱스를 찾아야 한다.
round(alpha, 3)## [1] 0.000 1.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000
## [12] 1.000 0.000 0.000 0.000 1.000 1.000 1.000 0.000 1.000 1.000 0.000
## [23] 1.000 0.000 1.000 1.000 1.000 0.000 0.000 1.000 0.000 0.000 0.000
## [34] 0.000 0.000 0.000 0.000 1.000 1.000 0.154 0.000 0.000 0.000 1.000
## [45] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 1.000
## [56] 0.000 1.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000
## [67] 0.000 1.000 1.000 0.000 0.000 0.154 0.000 1.000 0.000 1.000 1.000
## [78] 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000
## [89] 1.000 0.000 0.000 0.000 1.000 0.000 1.000 1.000 1.000 0.000 0.000
## [100] 0.000alpha를 보면 거의 0아니면 1인데 아닌 것들이 몇 개 보인다. 이 인덱스를 찾겠다는 것.
# indices for sv
sv.id <- which(1.0e-3 < alpha & alpha < (C - 1.0e-3))
sv.id## [1] 40 72- 확실하게 안쪽으로 들어와 있는 관측값들을 찾기 위해 작은 수를 더해줬다.
\(\beta = \overset{n}{\underset{i=1}{\Sigma}}\alpha_iy_ix_i\)
beta <- apply(alpha * y * x, 2, sum)
beta## [1] 1.526371 1.527540- 비가 같으면 기울기가 같으므로 두 개의 비만 같으면 된다.
\(\beta_0=y_i-\beta^Tx_i\)
temp <- y[sv.id] - x[sv.id,,drop = F] %*% beta
temp## [,1]
## [1,] -0.07797688
## [2,] -0.07797605beta0 <- mean(temp) # support vector들은 모두 동일한 b0
beta0## [1] -0.07797646plot(x[,1], x[,2], type = "n", xlab = "x1", ylab = "x2",
main = "Linear Support Vector Machine")
points(x[y == 1, 1], x[y == 1, 2], pch = 1, col = 2)
points(x[y != 1, 1], x[y != 1, 2], pch = 3, col = 4)
abline(a = -true.beta0/true.beta[1], b = -true.beta[1]/true.beta[2], lty = 2)
abline(a = -beta0/beta[1], -beta[1]/beta[2], col = 3, lwd = 2)
legend("topright", c("true", "estimated"), col = c(1,3), lty = c(2,1))
Reference
Trevor Hastie · Robert Tibshirani · Jerome Friedman(2009). Elment of Statistical Learning. Springer.
Gareth James , Daniela Witten, Trevor Hastie, Robert Tibshirani(2013). Introduction to Statistiacal Learning. Springer.