FPGA视觉从入门到放弃——懒人的支持向量机

支持向量机曾是机器学习领域中的主流方法。针对小样本，现在用起来依然很方便。同时，该方面的工具和教程多得数不清。所以这里引用老师木的话就很合适：“经常，有些事还没做就已经知道它无意义，于是就没做；有些事做不做都知道无意义，还是蠢蠢欲动”。

本篇围绕“用现有的库把自己的代码工作量减到最小”为话题，以比大多教程尽量简单的支持向量机理论和实践为例，讲述软件与硬件(DSP或FPGA等)结合时的偷懒方法。

1. 基本原理 5

(1) 超平面

假设数据集中有$l$个样本，每个样本包括属性向量和标签。所以，有样本$\lbrace (X^1,y^1)...(X^i,y^i),...(X^l,y^l) \rbrace$，$l$为样本的个数，$y^i\in \lbrace -1,1\rbrace$为标签，第$i$个样本$X^{i}=[x^i_1,...,x^i_n]$，$n$为每个样本的属性数目。

超平面的方程为：

$$H=w^TX+b=w_0x_0+w_1x_1+...+w_nx_n+b=0$$
其中，$w$为$n$维权重，b为偏置。

如果数据可分，则存在很多可分的超平面。贴着正负样本时的边缘超平面为$H=\pm1$，中间超平面为$H=0$。最优的超平面会使边缘超平面和中间超平面之间的间隔最大。

(2) 优化

第i个样本的属性向量$X^i$距离超平面$H$的距离为：

$$d(X^i,H)=\frac{w^TX^i}{||w||}$$

最大化间隔等价于优化问题：

$$\begin{align} &\min_{w,b} \frac{1}{2}||w||^2 \\&s.t.y_i(w^TX^i+b)\ge1, i=1,...,l \end{align}$$

• 超平面$H=\pm1$之间的距离为$\frac{2}{||w||}$，所以最大化间隔等价于最小化$\frac{||w||^2}{2}$。
• 当$y_i=1$时，$w^TX^i+b\ge1$，$x^i$位于$H=1$超平面的上侧；
• 当$y_i=-1$时，$w^TX^i+b\le1$，$x^i$位于$H=-1$超平面的下侧。

(3) 二次规划函数

X = quadprog(H,f,A,c)

二次规划函数求解关于向量X的二次规划问题：

$$\begin{align} &\min_{X}\frac{1}{2}X^THX+fX \\&s.t.AX\le c \end{align}$$

X仅表示二次规划要求解的项，与SVM中的样本数据没有关系。

(4) 编码优化问题

a. 目标函数

$$\begin{align} &\frac{1}{2}w^Tw \\&=\frac{1}{2}\left(\begin{matrix}w&b\end{matrix}\right) \left(\begin{matrix}I&0\\0&0\end{matrix}\right) \left(\begin{matrix}w\\b\end{matrix}\right) \\&=\frac{1}{2}\left(\begin{matrix}w_1&w_2...&w_n&b\end{matrix}\right) \left(\begin{matrix}1&0&0&0&0\\0&1&0&0&0\\0&0&...&0&0\\0&0&0&1&0\\0&0&0&0&0\end{matrix}\right) \left(\begin{matrix}w_1\\w_2\\...\\w_n\\b\end{matrix}\right) \end{align}$$

b. 约束条件

$$\begin{align} &\left(\begin{matrix}y_1(w^Tx^1+b)\\y_2(w^Tx^2+b)\\...\\y_l(w^Tx^l+b)\end{matrix}\right)=\left(\begin{matrix}y_1&0&0&0\\0&y_2&0&0\\0&0&...&0\\0&0&0&y_l\end{matrix}\right)\left(\begin{matrix}x_1^1&...&x_n^1&1\\x_1^2&...&x_n^2&1\\...&...&...&1\\x_1^l&...&x_n^l&1\end{matrix}\right)\left(\begin{matrix}w_1\\...\\w_n\\b\end{matrix}\right)\ge\left(\begin{matrix}1\\...\\1\\1\end{matrix}\right) \end{align}$$

2. 实现

(1) Matlab 2维版

这里把原优化问题作为二次规划问题求解。同理，也可以对它的对偶问题按二次规划问题求解。然而，该问题的规模正比于训练样本数，为避开高维数据计算的巨大开销，通常采用SMO方法[见 周志华, “机器学习”, pp. 124]。线性可分时，二次规划的解与工具箱的解 1相似。

clc;
clear;
close all;

%%

load fisheriris

X = meas(1:100,[2,3]);
group = species(1:100);
[groupIdx, groupStr] = grp2idx(group);

y = (groupIdx==1) * 2 - 1;

[l,n] = size(X);

H = eye(n + 1);
H(n + 1,n + 1) = 0;
f = zeros(n + 1,1);

Z = [X ones(l,1)];
A = -diag(y) * Z;
c = -1 * ones(l,1);

% 二次规划求解
w = quadprog(H,f,A,c);

% 调整横轴范围，便于与工具箱效果比较
X1 = [2:5];
w1 = w(1,1);
w2 = w(2,1);
b = w(3,1);

% 可分超平面：w1x1+w2x2+b=0 -> x2=-(w1*x1+b)/w2
Y1 = -(w1 * X1 + b) / w2;

XPos = X(find(y==1),:);
XNeg = X(find(y==-1),:);

% 超平面上界：w1x1+w2x2+b=1
YUP = (1 - w1 * X1 - b) / w2;
% 超平面上界：w1x1+w2x2+b=-1
YLOW = (-1 - w1 * X1 - b) / w2;

figure(1);
set(gcf,'Color',[1,1,1]);
subplot(1,2,1);
plot(XPos(:,1),XPos(:,2),'r+'); hold on;
plot(XNeg(:,1),XNeg(:,2),'g*'); hold on;
plot(X1,Y1,'k-'); hold on;
plot(X1,YUP,'m:'); hold on;
plot(X1,YLOW,'m:'); 
legend('setosa','versicolor','hyperplane','upper margin','lower margin');

% 调整纵轴范围，便于与工具箱效果比较
ylim([1,5.5]);

title('svm trained with quadratic programming');
xlabel('sepal length');
ylabel('sepal width');
grid on;

subplot(1,2,2);
svmStruct = svmtrain(X,group,'ShowPlot',true);
title('svm trained with matlab toolbox');
xlabel('sepal length');
ylabel('sepal width');
grid on;

注： 这里二次规划求解要求训练样本线性可分，否则找不到可行解。

(2) Matlab n维核函数版

工具箱自带预测器用到的参数 2有：

Alpha($\alpha$)——拉格朗日乘子
Beta($\beta$)——线性预测器的系数
Bias($b$)——偏置
KernelParameters.Scale($s$)——缩放因子

如果KernelParameters.Function为”Linear“，则输出为：

$$f(x)=\left(\frac{x}{s}\right)'\beta+b$$

更一般地，带核函数的预测输出可表示为：

$$f(x)=\sum^{m}_{i=1}\alpha_iy_iK(x,x_i)+b$$
其中， $y_i$为第$i$个支持向量$x_i$的标签，$m$为训练得到的支持向量的总数，$K(x,x_i)$为测试样本$x$与$x_i$间的核函数输出。

clc;
clear;
close all;
%%

load fisheriris

X = meas(1:100,:);
group = species(1:100);
[groupIdx, groupStr] = grp2idx(group);

Y = (groupIdx==2) * 2 - 1;

XPos = X(find(Y==1),:);
XNeg = X(find(Y==-1),:);
%%

cv = cvpartition(Y,'k',10);
err = zeros(cv.NumTestSets,1);

for i = 1:cv.NumTestSets

    trainIdx = cv.training(i);
    testIdx = cv.test(i);
    Xtrain = X(trainIdx,:);
    Ytrain = Y(trainIdx);
    Xtest = X(testIdx,:);
    Ytest = Y(testIdx);

    %% SVM 训练
    clf = fitcsvm(Xtrain,Ytrain,'KernelFunction','polynomial','ClassNames',[-1,1]);
    [m,n] = size(clf.SupportVectors);

    %% SVM 验证

    % (0) 线性核 'linear'
    % Ypred = (Xtest / clf.KernelParameters.Scale) * clf.Beta + clf.Bias > 0;

    Ypred = zeros(length(Xtest),1);
    for j=1:m
        % (1) 线性核 'linear'
        % Kernel = Xtest * clf.SupportVectors(j,:)'; 
        % (2) 高斯核 'rbf'
        % Kernel = exp(-sum((Xtest - repmat(clf.SupportVectors(j,:),length(Xtest),1)).^2,2));
        % (3) 多项式核 'polynomial'
        Kernel = (1 + Xtest * clf.SupportVectors(j,:)').^clf.ModelParameters.KernelPolynomialOrder;

        Ypred = Ypred + clf.Alpha(j) * clf.SupportVectorLabels(j) * Kernel;
    end
    Ypred = (Ypred + clf.Bias > 0) * 2 - 1;    

    err(i) = sum(~strcmp(Ypred,Ytest));
end
cvError = sum(err)/sum(cv.TestSize)
%%

figure(1);
set(gcf,'Color',[1,1,1]);
%subplot(1,1,1);
plot(XPos(:,1),XPos(:,2),'ro'); hold on;
plot(XNeg(:,1),XNeg(:,2),'go'); hold on;
plot(Xtest(:,1),Xtest(:,2),'bo'); hold on;
plot(clf.SupportVectors(:,1),clf.SupportVectors(:,2),'k^'); hold on;

[l,n] = size(Xtest);

for i=1:l
    if Ypred(i)==1
        plot(Xtest(i,1),Xtest(i,2),'r*'); hold on;
    else
        plot(Xtest(i,1),Xtest(i,2),'g*'); hold on;        
    end
end

title('prediction with svm classifier');
xlabel('sepal length');
ylabel('sepal width');
grid on;

10折交叉验证 3后最后1折多项式核 4的预测结果如下图。其中三角形为支持向量，蓝色圆圈为验证样本，星型颜色和周围样本的颜色一致则分类正确。

(3) 转换

Matlab离线训练后得到分类器clf。线性核与高斯核预测时，保存支持向量即可；多项式核预测时，另外保存多项式的阶数。根据需要对预测部分转换成项目约束的语言即可。

---