如何使用Matlab的quadprog实现软间隔支持向量机模型？

Question

假设我们有一个训练数据集 {yᵢ, xᵢ}，其中 i = 1, ..., n，yᵢ 可以是 -1 或 1，而 xᵢ 可以是二维或三维点等。

一般情况下，当输入点是线性可分的时，SVM 模型可以定义如下：

min 1/2*||w||²
w,b

受限于约束条件（对于 i = 1, ..., n）

yᵢ*(w*xᵢ - b) >= 1

这通常被称为硬间隔SVM模型，因此是一个约束最小化问题，其中未知数是w和b。我们也可以在要最小化的函数中省略1/2，因为它只是一个常数。

现在，关于Matlab的quadprog的文档说明如下：

x = quadprog(H, f, A, b) 最小化 1/2*x'*H*x + f'*x，同时满足限制条件 A*x ≤ b。其中A是一个双精度矩阵，b是一个双精度向量。

我们可以使用quadprog函数来实现硬间隔SVM模型，以获得权重向量w，如下所示：

• H 变成了一个单位矩阵。
• f' 变成了一个零矩阵。
• A 是约束条件的左侧。
• b 等于 -1，因为原始约束条件为 >= 1，当我们在两侧乘以 -1 时，它变成了 <= -1。

---

现在，我正在尝试实现一个软间隔支持向量机模型。这里的最小化方程式是：

min (1/2)*||w||² + C*(∑ ζᵢ)
w,b

受限于约束条件（对于 i = 1, ..., n）

yᵢ*(w*xᵢ - b) >= 1 - ζᵢ

这个优化问题需要满足ζᵢ >= 0，其中∑是求和符号，ζᵢ = max(0, 1 - yᵢ*(w*xᵢ - b))，C是一个超参数。

如何使用Matlab的quadprog函数解决这个优化问题？我不清楚该方程应如何映射到quadprog函数的参数。

软间隔支持向量机模型的"原始形式"（即上述定义）可以转换为"对偶形式"。我已经完成了转换，并能够获取拉格朗日变量值（在对偶形式中）。但是，我想知道是否可以直接使用quadprog来解决原始形式，而无需转换为对偶形式。

- user1067334

Answer 1

我认为这不会成为问题。让我们定义向量 z 包含 (2n + 1) 个变量：

z = (w, eps, b)

然后，H 变成一个对角矩阵，其对角线上的前 n 个值为 1，最后一个值为 n+1 的值为零：

H = diag([ones(1, n), zeros(1, n + 1)])

向量f可以表示为：

f = [zeros(1, n), C * ones(1, n), 0]'

第一组限制条件如下：

Aineq = [A1, eye(n), zeros(n, 1)]
bineq = ones(n, 1)

其中A1与原始形式中的矩阵相同。

第二组约束成为下限：

lb = (inf(n, 1), zeros(n, 1), inf(n, 1))

然后你可以调用MATLAB:

z = quadprog(H, f, Aineq, bineq, [], [], lb);

备注：在某些小细节上我可能会犯错，但总体思路是正确的。

- vharavy

---

"n"变量的使用方式让我更加困惑，但我大致明白了。在很多地方，n变量并不匹配。 - user1067334

@user1067334 n代表的是特征数还是训练样本数？ - tusharfloyd

Answer 2

我想澄清@vharavy的答案，因为在尝试推断他代码中'n'的含义时可能会迷失方向。根据他的回答和SVM维基百科文章，这是我的版本。我假设我们有一个名为"test.dat"的文件，其中包含测试点的坐标以及它们在最后一列的类成员资格。

"test.dat"的示例内容如下，其中包含3D点：

-3,-3,-2,-1
-1,3,2,1
5,4,1,1
1,1,1,1
-2,5,4,1
6,0,1,1
-5,-5,-3,-1
0,-6,1,-1
-7,-2,-2,-1

这里是代码：

data = readtable("test.dat");
tableSize = size(data);
numOfPoints = tableSize(1);
dimension = tableSize(2) - 1;
PointsCoords = data(:, 1:dimension);
PointsSide = data.(dimension+1);

C = 0.5; %can be changed
n = dimension;
m = numOfPoints; %can be also interpretet as number of constraints
%z = [w, eps, b]; number of variables in 'z' is equal to n + m + 1
H = diag([ones(1, n), zeros(1, m + 1)]);
f = [zeros(1, n), C * ones(1, m), 0];
Aineq = [-diag(PointsSide)*table2array(PointsCoords), -eye(m), PointsSide];
bineq = -ones(m, 1);
lb = [-inf(1, n), zeros(1, m), -inf];
z = quadprog(H, f, Aineq, bineq, [], [], lb);

- JAJA

Answer 3

完整的代码。思路与上面相同。

n = size(X,1);
m = size(X,2);
H = diag([ones(1, m), zeros(1, n + 1)]);
f = [zeros(1,m+1) c*ones(1,n)]';
p = diag(Y) * X;
A = -[p Y eye(n)];
B = -ones(n,1);
lb = [-inf * ones(m+1,1) ;zeros(n,1)];
z = quadprog(H,f,A,B,[],[],lb);
w = z(1:m,:);
b = z(m+1:m+1,:);
eps = z(m+2:m+n+1,:);

- ashwin kumar

---

这份代码无法运行。您能提供一个可运行的示例吗？例如，您应该提供 X 等。 - nbro

Answer 4

如果让z = (w; w0; eps)T成为一个具有n+1+m个元素的长向量。(其中m是点的数量) 那么，

H= diag([ones(1,n),zeros(1,m+1)]).
f = [zeros(1; n + 1); ones(1;m)]

不等式约束可以指定为：

A = -diag(y)[X; ones(m; 1); zeroes(m;m)] -[zeros(m,n+1),eye(m)],

其中X是原始形式下的n x m输入矩阵。在A的两个部分中，第一部分是为了w0，第二部分是为了eps。

b = ones(m,1) 

等式约束：

Aeq = zeros(1,n+1 +m)
beq = 0

边界：

lb = [-inf*ones(n+1,1); zeros(m,1)]
ub = [inf*ones(n+1+m,1)]

现在，z=quadprog(H,f,A,b,Aeq,beq,lb,ub)。

- akuriako