一般步骤
- 将一般形式的待求解目标函数化成标准形式。
标准形式如下:
{ min x = 1 2 x T P x + q T x s . t . G x ≤ h A x = b \begin{cases} \min \limits_{x} = \frac{1}{2}x^{T}Px + q^{T}x \\ s.t. \ \ \ \ Gx \leq h \\ Ax = b \end{cases} ⎩⎪⎨⎪⎧xmin=21xTPx+qTxs.t. Gx≤hAx=b
- 带入cvxopt包中的solvers方法求解
引例
【例】求如下的二次规划问题
min
x
f
(
x
)
=
1
2
x
1
2
+
x
2
2
−
x
1
x
2
−
2
x
1
−
6
x
2
,
s
.
t
.
{
x
1
+
x
2
≤
2
−
x
1
+
2
x
2
≤
2
2
x
1
+
x
2
≤
3
x
1
,
x
2
≥
0
\min \limits_{x}f(x) = \frac{1}{2}x_{1}^{2} + x_{2}^{2} -x_{1}x_{2} -2x_{1} - 6x_{2},s.t.\ \begin{cases} x_{1} + x_{2} \leq 2\\ -x_{1} + 2x_{2} \leq 2\\ 2x_{1} + x_{2}\leq3\\ x_{1}, x_{2} \geq 0 \end{cases}
xminf(x)=21x12+x22−x1x2−2x1−6x2,s.t. ⎩⎪⎪⎪⎨⎪⎪⎪⎧x1+x2≤2−x1+2x2≤22x1+x2≤3x1,x2≥0
首先,我们将上式化成标准形式。
向量
x
\boldsymbol{x}
x 很容易写出来,因为
f
(
x
)
f(x)
f(x) 包含两个变量
x
1
x_{1}
x1,
x
2
x_{2}
x2,因此
x
=
[
x
1
x
2
]
x = \begin{bmatrix} x_{1} \\ x_{2} \end{bmatrix}
x=[x1x2]
向量
q
q
q只与两个变量
x
1
x_{1}
x1,
x
2
x_{2}
x2的一次项有关,所以
q
T
x
=
−
2
x
1
−
6
x
2
q^T \boldsymbol{x} = -2x_{1} - 6x_{2}
qTx=−2x1−6x2,因此
q
=
[
−
2
−
6
]
q = \begin{bmatrix} -2 \\ -6 \end{bmatrix}
q=[−2−6]
最后,矩阵
P
P
P 只与两个变量
x
1
x_{1}
x1,
x
2
x_{2}
x2的二次项有关,所以
1
2
x
T
P
x
=
1
2
x
1
2
+
x
2
2
−
x
1
x
2
\frac{1}{2}x^TPx = \frac{1}{2}x_{1}^2 + x_{2}^2 - x_{1}x_{2}
21xTPx=21x12+x22−x1x2,这里要注意的是不同于二次型,这里有个系数
1
2
\frac{1}{2}
21,所以矩阵
P
P
P的元素是二次型中的矩阵元素大小的两倍。给出一个规律:设矩阵
P
P
P第
i
i
i行第
j
j
j列的元素大小为
P
(
i
,
j
)
P(i, j)
P(i,j),二次项
x
i
,
y
j
x_{i}, y_{j}
xi,yj的系数为
a
(
i
,
j
)
a(i, j)
a(i,j),
P
(
i
,
j
)
=
{
2
a
(
i
,
j
)
,
i
=
j
a
(
i
,
j
)
,
i
≠
j
P(i, j) = \begin{cases} 2a(i, j), i = j \\ a(i, j), i \neq j \end{cases}
P(i,j)={2a(i,j),i=ja(i,j),i=j
本例中,
P
(
i
,
j
)
=
[
1
−
2
−
2
2
]
P(i, j ) = \begin{bmatrix} 1 &-2\\ -2&2 \end{bmatrix}
P(i,j)=[1−2−22],这是由于
x
1
x_{1}
x1的平方项(即
x
1
2
x_{1}^2
x12),所以第1行的第1列的元素为
1
=
2
∗
(
1
2
)
1 = 2*(\frac{1}{2})
1=2∗(21),
x
2
x_{2}
x2的平方项(即
x
2
2
x_{2}^2
x22)系数为1,所以第2行第2列的元素为
2
=
2
∗
1
2 = 2*1
2=2∗1,
x
1
x
2
x_{1}x_{2}
x1x2即(
x
2
x
1
x_{2}x_{1}
x2x1)的系数为-1,所以第1行第2列和第2行第1列的元素均为-2。
再看约束条件部分,约束条件应该写成以下标准形式:
{
G
x
≤
h
A
x
=
b
\begin{cases} Gx \leq h \\ Ax = b \end{cases}
{Gx≤hAx=b
本例中约束条件只有不等式约束,因此
A
=
∅
A = \emptyset
A=∅,对于
G
G
G 和
h
h
h很容易就看得出来:
G
=
[
1
1
−
1
2
2
1
−
1
0
0
−
1
]
,
h
=
[
2
2
3
0
0
]
G = \begin{bmatrix} 1 & 1\\ -1 & 2\\ 2 & 1\\ -1 & 0\\ 0 & -1 \end{bmatrix}, h = \begin{bmatrix} 2\\ 2\\ 3\\ 0\\ 0 \end{bmatrix}
G=⎣⎢⎢⎢⎢⎡1−12−101210−1⎦⎥⎥⎥⎥⎤,h=⎣⎢⎢⎢⎢⎡22300⎦⎥⎥⎥⎥⎤
注意,当
x
≤
0
x \leq 0
x≤0时,乘以一个
−
1
-1
−1,变成
−
x
≥
0
-x \geq 0
−x≥0。
python 代码求解:1)要是约束条件中有的没有,则不输入给qb函数即可,比如此例子中的A,b。2)输入必须是float类型,所以用tc=‘d’,转一下数据类型。
def sqp(paraP, paraq, paraG, parah):
P = matrix(np.array(paraP), tc='d')
q = matrix(np.array(paraq), tc='d')
G = matrix(np.array(paraG), tc='d')
h = matrix(np.array(parah), tc='d')
# A = matrix(np.array(paraA), tc='d')
# b = matrix(np.array(parab), tc='d')
result = solvers.qp(P, q, G, h)
print('x\n', result['x'])
P = [[1, -1], [-1, 2]]
q = [-2, 6]
G = [[1, 1], [-1, 2], [2, 1], [-1, 0], [0, -1]]
h = [2, 2, 3, 0, 0]
sqp(P, q, G, h)
求解结果:
SVDD求解
SVDD原理是在特征空间寻找一个体积最小的超球体,为了构造这样一个最小体积的超球体,SVDD需要解决以下优化问题:
{
min
x
,
R
,
ξ
R
2
+
C
∑
i
=
1
n
ξ
i
s
.
t
.
∣
∣
ϕ
(
x
i
)
−
a
∣
∣
2
≤
R
+
ξ
i
,
ξ
i
≥
0
,
∀
i
=
1
,
2
,
⋯
,
n
\begin{cases} \min \limits_{x, R, \xi} R^2 + C\sum\limits_{i = 1}^n\xi_{i} \\ s.t. \ \ \ \ ||\phi(x_{i}) - a ||^2 \leq R + \xi_{i}, \xi_{i} \geq 0, \forall i = 1, 2, \cdots , n \end{cases}
⎩⎨⎧x,R,ξminR2+Ci=1∑nξis.t. ∣∣ϕ(xi)−a∣∣2≤R+ξi,ξi≥0,∀i=1,2,⋯,n
式中,
R
R
R 是超球体半径,
a
a
a 是超球体的球心,
ξ
\xi
ξ 是松弛因子,
C
C
C 是一个权衡超球体体积和误分率的惩罚参数,
C
C
C 大,则表示惩罚越大,为了使目标函数最小,则
ξ
\xi
ξ 变小,同时
R
R
R 变大。相反,
C
C
C 越小,
R
R
R 越小。
单分类
对于单分类来说,结合拉格朗日乘子法(具体过程详见参考文献),原问题的对偶问题为:
{
min
α
i
∑
i
=
1
n
∑
j
=
1
n
α
i
α
j
K
(
x
i
,
x
j
)
−
∑
i
=
1
n
α
i
K
(
x
i
,
x
j
)
s
.
t
.
0
≤
α
i
≤
C
,
∑
i
=
1
n
α
i
=
1
\begin{cases} \min \limits_{\alpha_{i}}\sum\limits_{i=1}^n\sum\limits_{j=1}^n\alpha_{i}\alpha_{j} K(x_{i}, x_{j}) - \sum\limits_{i = 1}^n\alpha_{i}K(x_{i}, x_{j})\\ s.t. \ \ \ \ 0 \leq \alpha_{i} \leq C, \ \ \sum\limits_{i =1}^n\alpha_{i} = 1 \end{cases}
⎩⎪⎪⎨⎪⎪⎧αimini=1∑nj=1∑nαiαjK(xi,xj)−i=1∑nαiK(xi,xj)s.t. 0≤αi≤C, i=1∑nαi=1
把所有的
α
i
,
α
j
\alpha_{i}, \alpha_{j}
αi,αj 看成是一个向量
α
\boldsymbol{\alpha}
α,根据二次规划的标准形式可以得到:
P
=
2
∑
i
=
1
n
∑
j
=
1
n
K
(
x
i
,
x
j
)
,
q
=
(
−
∑
i
=
1
n
K
(
x
i
,
x
i
)
)
T
G
=
[
1
0
0
⋯
0
0
1
0
⋯
0
⋮
⋱
⋱
⋱
⋮
0
0
0
⋯
1
−
1
0
0
⋯
0
0
−
1
0
⋯
0
⋮
⋱
⋱
⋱
⋮
0
0
0
⋯
−
1
]
,
h
=
[
C
⋮
C
0
⋮
0
]
A
=
[
1
1
⋯
1
]
b
=
1
P = 2\sum\limits_{i = 1}^n\sum\limits_{j = 1}^nK(x_{i}, x_{j}),\ q =(-\sum\limits_{i =1}^nK(x_{i}, x_{i}))^T\\ G = \begin{bmatrix} 1 & 0 &0 &\cdots &0\\ 0 & 1& 0&\cdots &0\\ \vdots & \ddots& \ddots&\ddots & \vdots\\ 0 & 0& 0&\cdots &1\\ -1 & 0 &0 &\cdots &0\\ 0 & -1& 0&\cdots &0\\ \vdots & \ddots& \ddots&\ddots & \vdots\\ 0 & 0& 0&\cdots &-1 \end{bmatrix}, \ h = \begin{bmatrix} C \\ \vdots\\ C\\ 0\\ \vdots\\ 0\end{bmatrix}\\ A = \begin{bmatrix} 1 & 1& \cdots &1 \end{bmatrix} \ b =1
P=2i=1∑nj=1∑nK(xi,xj), q=(−i=1∑nK(xi,xi))TG=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡10⋮0−10⋮001⋱00−1⋱000⋱000⋱0⋯⋯⋱⋯⋯⋯⋱⋯00⋮100⋮−1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤, h=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡C⋮C0⋮0⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤A=[11⋯1] b=1
设样本总数是
n
n
n,其中
G
G
G 大小是
2
n
∗
n
2n * n
2n∗n,
h
h
h 的大小是
2
n
∗
1
2n * 1
2n∗1,
A
A
A 的大小是
1
∗
n
1 * n
1∗n,
b
b
b 的大小时1。
二分类
对于二分类,在正类训练集中加入了少数的负类样本来防止过拟合情况,假设正样例和负样例的标签分别是
y
i
=
+
1
y_{i} = +1
yi=+1 和
y
j
=
−
1
y_{j} = -1
yj=−1,则原优化问题的对偶问题变为:
{
min
α
i
∑
i
=
1
n
∑
j
=
1
n
y
i
y
j
α
i
α
j
K
(
x
i
,
x
j
)
−
∑
i
=
1
n
y
i
α
i
K
(
x
i
,
x
j
)
s
.
t
.
0
≤
α
i
≤
C
1
,
0
≤
α
i
≤
C
2
,
∑
i
=
1
n
y
i
α
i
=
1
\begin{cases} \min \limits_{\alpha_{i}}\sum\limits_{i=1}^n\sum\limits_{j=1}^n y_{i}y_{j} \alpha_{i}\alpha_{j} K(x_{i}, x_{j}) - \sum\limits_{i = 1}^ny_{i}\alpha_{i}K(x_{i}, x_{j})\\ s.t. \ \ \ \ 0 \leq \alpha_{i} \leq C_{1}, \ 0 \leq \alpha_{i} \leq C_{2}, \ \ \sum\limits_{i =1}^ny_{i}\alpha_{i} = 1 \end{cases}
⎩⎪⎪⎨⎪⎪⎧αimini=1∑nj=1∑nyiyjαiαjK(xi,xj)−i=1∑nyiαiK(xi,xj)s.t. 0≤αi≤C1, 0≤αi≤C2, i=1∑nyiαi=1
两类比一类会多一个
y
y
y,此时可以的到:
P
=
2
∑
i
=
1
n
∑
j
=
1
n
y
i
y
j
K
(
x
i
,
x
j
)
,
q
=
(
−
∑
i
=
1
n
y
i
K
(
x
i
,
x
i
)
)
T
G
=
[
1
0
0
⋯
0
0
1
0
⋯
0
⋮
⋱
⋱
⋱
⋮
0
0
0
⋯
1
−
1
0
0
⋯
0
0
−
1
0
⋯
0
⋮
⋱
⋱
⋱
⋮
0
0
0
⋯
−
1
]
,
h
=
[
C
1
⋮
C
2
⋮
0
]
A
=
[
y
1
,
y
2
,
⋯
,
y
n
]
,
b
=
1
P = 2\sum\limits_{i = 1}^n\sum\limits_{j = 1}^n y_{i} y_{j}K(x_{i}, x_{j}),\ q =(-\sum\limits_{i =1}^ny_{i}K(x_{i}, x_{i}))^T\\ G = \begin{bmatrix} 1 & 0 &0 &\cdots &0\\ 0 & 1& 0&\cdots &0\\ \vdots & \ddots& \ddots&\ddots & \vdots\\ 0 & 0& 0&\cdots &1\\ -1 & 0 &0 &\cdots &0\\ 0 & -1& 0&\cdots &0\\ \vdots & \ddots& \ddots&\ddots & \vdots\\ 0 & 0& 0&\cdots &-1 \end{bmatrix}, \ h = \begin{bmatrix}C_{1} \\ \vdots\\ C_{2} \\ \vdots\\ 0\end{bmatrix}\\ A = \begin{bmatrix} y_{1}, y_{2}, \cdots, y_{n} \end{bmatrix}, \ b =1
P=2i=1∑nj=1∑nyiyjK(xi,xj), q=(−i=1∑nyiK(xi,xi))TG=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡10⋮0−10⋮001⋱00−1⋱000⋱000⋱0⋯⋯⋱⋯⋯⋯⋱⋯00⋮100⋮−1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤, h=⎣⎢⎢⎢⎢⎢⎢⎡C1⋮C2⋮0⎦⎥⎥⎥⎥⎥⎥⎤A=[y1,y2,⋯,yn], b=1
与单分类不同的是
h
h
h 的前
n
n
n 部分的
C
C
C 是对应于正负样本的
C
1
,
C
2
C_{1}, C_{2}
C1,C2,总大小和单分类一样为
2
n
∗
1
2n *1
2n∗1。
A
A
A 等于标签的值。
代码分析
label = np.mat(label)
K = np.multiply(label*label.T, K)
# P
n = K.shape[0]
P = K+K.T
# q
q = -np.multiply(label, np.mat(np.diagonal(K)).T)
# G
G1 = -np.eye(n)
G2 = np.eye(n)
G = np.append(G1, G2, axis=0)
# h
h1 = np.mat(np.zeros(n)).T # lb
h2 = np.mat(np.ones(n)).T
if self.labeltype == 'single':
h2[label == 1] = self.parameters["positive penalty"]
if self.labeltype == 'hybrid':
h2[label == 1] = self.parameters["positive penalty"]
h2[label == -1] = self.parameters["negative penalty"]
h = np.append(h1, h2, axis=0)
# A, b
A = np.mat(np.ones(n) * np.array(label).reshape(1, n))
#A = np.mat(np.ones(n))
b = 1.
#
P = matrix(P)
q = matrix(q)
G = matrix(G)
h = matrix(h)
A = matrix(A)
b = matrix(b)
#
sol =solvers.qp(P, q, G, h, A, b)
alf = np.array(sol['x'])
print(alf)