Processing math: 100%
Liu Shouda coder

svm derive

2017-04-16
ml

svm推导,详见pluskid博客。

基本推导

smo优化

α1,α2的取值在如图所示的直线上。假定α2的上限下限分别是H,L

y1y2时,α1α2=k

α1α2=CH 所以: H=Cα1+α2

α1α2=0L 所以: L=α2α1

同理,

y1==y2时, α1+α2=k

α1+α2=+H, 所以,H=α1+α2

α1+α2=C+L, 所以,L=α1+α2C

αnew2=α2+y2(E1E2)η

其中

η=k(x1,x2)+k(x2,x2)2k(x1,x2) Ei=predictiyi

因为,α2必须在L,H的范围内,所以如果超出范围需要进行截断clip. αnew,clip2={Hif αnew2Hαnew2if L<αnew2<HLαnew2L α1α2是线性关系,所以,

αnew1α1=y1y2(αnew,clip2α2)

所以,

αnew1=α1+y1y2(αnew,clip2α2)

求b:其实就是希望b的取值能够使得当前的预测值与标签相同。

f(x)=Nj=1yjαjK(xj,x)b

yjα1K(xj,x)=sj

所以,f(x1)=s1+s2+constantb

f(x1)new=snew1+snew2+constantbnew

求差值: f(x1)newf(x1)=snew1s1+snew2s2(bnewb)

我们希望这个这个差值能够弥补之前的预测误差,即,

E1+f(x1)newf(x1)=E1+snew1s1+snew2s2(bnewb)=0

从上式子可以求出b的一个取值b1

同理,对于样本x2同样可以获得一个取值b2.

最后采用策略选取b1,b2.

  • 如果a1在边界上,则取b1
  • 否则,如果a2在边界上,则取b2
  • 否则,取b1,b2的均值。

代码

来自MLAlgorithm

  1. # coding:utf-8
  2. from mla.base import BaseEstimator
  3. from mla.svm.kernerls import Linear
  4. import numpy as np
  5. import logging
  6. np.random.seed(9999)
  7. """
  8. References:
  9. The Simplified SMO Algorithm http://cs229.stanford.edu/materials/smo.pdf
  10. """
  11. class SVM(BaseEstimator):
  12. def __init__(self, C=1.0, kernel=None, tol=1e-3, max_iter=100):
  13. """Support vector machines implementation using simplified SMO optimization.
  14. Parameters
  15. ----------
  16. C : float, default 1.0
  17. kernel : Kernel object
  18. tol : float , default 1e-3
  19. max_iter : int, default 100
  20. """
  21. self.C = C
  22. self.tol = tol
  23. self.max_iter = max_iter
  24. if kernel is None:
  25. self.kernel = Linear()
  26. else:
  27. self.kernel = kernel
  28. self.b = 0
  29. self.alpha = None
  30. self.K = None
  31. def fit(self, X, y=None):
  32. self._setup_input(X, y)
  33. self.K = np.zeros((self.n_samples, self.n_samples))
  34. for i in range(self.n_samples):
  35. self.K[:, i] = self.kernel(self.X, self.X[i, :])
  36. self.alpha = np.zeros(self.n_samples)
  37. self.sv_idx = np.arange(0, self.n_samples)
  38. return self._train()
  39. def _train(self):
  40. iters = 0
  41. while iters < self.max_iter:
  42. iters += 1
  43. alpha_prev = np.copy(self.alpha)
  44. for j in range(self.n_samples):
  45. # Pick random i
  46. i = self.random_index(j)
  47. eta = 2.0 * self.K[i, j] - self.K[i, i] - self.K[j, j]
  48. if eta >= 0:
  49. continue
  50. L, H = self._find_bounds(i, j)
  51. # Error for current examples
  52. e_i, e_j = self._error(i), self._error(j)
  53. # Save old alphas
  54. alpha_io, alpha_jo = self.alpha[i], self.alpha[j]
  55. # Update alpha
  56. self.alpha[j] -= (self.y[j] * (e_i - e_j)) / eta
  57. self.alpha[j] = self.clip(self.alpha[j], H, L)
  58. self.alpha[i] = self.alpha[i] + self.y[i] * self.y[j] * (alpha_jo - self.alpha[j])
  59. # Find intercept
  60. b1 = self.b - e_i - self.y[i] * (self.alpha[i] - alpha_jo) * self.K[i, i] - \
  61. self.y[j] * (self.alpha[j] - alpha_jo) * self.K[i, j]
  62. b2 = self.b - e_j - self.y[j] * (self.alpha[j] - alpha_jo) * self.K[j, j] - \
  63. self.y[i] * (self.alpha[i] - alpha_io) * self.K[i, j]
  64. if 0 < self.alpha[i] < self.C:
  65. self.b = b1
  66. elif 0 < self.alpha[j] < self.C:
  67. self.b = b2
  68. else:
  69. self.b = 0.5 * (b1 + b2)
  70. # Check convergence
  71. diff = np.linalg.norm(self.alpha - alpha_prev)
  72. if diff < self.tol:
  73. break
  74. logging.info('Convergence has reached after %s.' % iters)
  75. # Save support vectors index
  76. self.sv_idx = np.where(self.alpha > 0)[0]
  77. def _predict(self, X=None):
  78. n = X.shape[0]
  79. result = np.zeros(n)
  80. for i in range(n):
  81. result[i] = np.sign(self._predict_row(X[i, :]))
  82. return result
  83. def _predict_row(self, X):
  84. k_v = self.kernel(self.X[self.sv_idx], X)
  85. return np.dot((self.alpha[self.sv_idx] * self.y[self.sv_idx]).T, k_v.T) + self.b
  86. def clip(self, alpha, H, L):
  87. if alpha > H:
  88. alpha = H
  89. if alpha < L:
  90. alpha = L
  91. return alpha
  92. def _error(self, i):
  93. """Error for single example."""
  94. return self._predict_row(self.X[i]) - self.y[i]
  95. def _find_bounds(self, i, j):
  96. """Find L and H such that L <= alpha <= H.
  97. Also, alpha must satisfy the constraint 0 <= αlpha <= C.
  98. """
  99. if self.y[i] != self.y[j]:
  100. L = max(0, self.alpha[j] - self.alpha[i])
  101. H = min(self.C, self.C - self.alpha[i] + self.alpha[j])
  102. else:
  103. L = max(0, self.alpha[i] + self.alpha[j] - self.C)
  104. H = min(self.C, self.alpha[i] + self.alpha[j])
  105. return L, H
  106. def random_index(self, z):
  107. i = z
  108. while i == z:
  109. i = np.random.randint(0, self.n_samples - 1)
  110. return i

参考:

http://blog.pluskid.org/?page_id=683

http://cs229.stanford.edu/notes/cs229-notes3.pdf

http://cs229.stanford.edu/materials/smo.pdf

https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/tr-98-14.pdf


下一篇 mxnet graph解析