一般来说，一个完整的机器学习项目分为以下步骤：
# 引入相关科学计算包
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use("ggplot")
import seaborn as sns
(1) 收集数据集并选择合适的特征： （第二步）
数据集依旧是Boston房价数据集，简单，处理快。
from sklearn import datasets
boston = datasets.load_boston() # 返回一个类似于字典的类
X = boston.data #特征X的矩阵
y = boston.target # 因变量的向量
features = boston.feature_names # 特征名称
boston_data = pd.DataFrame(X,columns=features)
boston_data["Price"] = y
boston_data.head()
| CRIM | ZN | INDUS | CHAS | NOX | RM | AGE | DIS | RAD | TAX | PTRATIO | B | LSTAT | Price | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.00632 | 18.0 | 2.31 | 0.0 | 0.538 | 6.575 | 65.2 | 4.0900 | 1.0 | 296.0 | 15.3 | 396.90 | 4.98 | 24.0 |
| 1 | 0.02731 | 0.0 | 7.07 | 0.0 | 0.469 | 6.421 | 78.9 | 4.9671 | 2.0 | 242.0 | 17.8 | 396.90 | 9.14 | 21.6 |
| 2 | 0.02729 | 0.0 | 7.07 | 0.0 | 0.469 | 7.185 | 61.1 | 4.9671 | 2.0 | 242.0 | 17.8 | 392.83 | 4.03 | 34.7 |
| 3 | 0.03237 | 0.0 | 2.18 | 0.0 | 0.458 | 6.998 | 45.8 | 6.0622 | 3.0 | 222.0 | 18.7 | 394.63 | 2.94 | 33.4 |
| 4 | 0.06905 | 0.0 | 2.18 | 0.0 | 0.458 | 7.147 | 54.2 | 6.0622 | 3.0 | 222.0 | 18.7 | 396.90 | 5.33 | 36.2 |
各个特征的相关解释：
(2) 选择度量模型性能的指标： （第三步）
均方误差反映估计量与被估计量之间差异程度的一种度量,对无偏估计使用方差，对有偏估计使用均方误差。均方误差是评价点估计的最一般的标准。
平均绝对误差:所有单个观测值与算术平均值的偏差的绝对值的平均,与平均误差相比，平均绝对误差由于离差被绝对值化，不会出现正负相抵消的情况，因而，平均绝对误差能更好地反映预测值误差的实际情况。
决定系数:反应因变量的全部变异能通过回归关系被自变量解释的比例。如R平方为0.8，则表示回归关系可以解释因变量80%的变异。反之，如果我们能控制自变量不变，则因变量的变异程度会减少80%。R方介于0~1之间，越接近1，回归拟合效果越好，一般认为超过0.8的模型拟合优度比较高。
如果R2=解释方差得分，则意味着：平均误差=0，平均误差反映了我们估计的趋势，即：有偏vs无偏估计。
https://scikit-learn.org/stable/modules/model_evaluation.html#regression-metrics
在这个案例中，我们使用MSE均方误差为模型的性能度量指标。
(3) 选择具体的模型并进行训练（第四步）
回归分析是一种预测性的建模技术，它研究的是因变量（目标）和自变量（特征）（多个）之间的因果/显著关系（如价格变动与促销活动数量之间联系）。通常使用曲线/线来拟合数据点，目标是使曲线到数据点的距离差异最小。
而线性回归就是回归问题中的一种，线性回归假设目标值与特征之间线性相关，即满足一个多元一次方程。通过构建损失函数，来求解损失函数最小时的参数w ：
假设：数据集
D
=
{
(
x
1
,
y
1
)
,
.
.
.
,
(
x
N
,
y
N
)
}
D = \{(x_1,y_1),...,(x_N,y_N) \}
D={(x1,y1),...,(xN,yN)}，
x
i
∈
R
p
,
y
i
∈
R
,
i
=
1
,
2
,
.
.
.
,
N
x_i \in R^p,y_i \in R,i = 1,2,...,N
xi∈Rp,yi∈R,i=1,2,...,N，
X
=
(
x
1
,
x
2
,
.
.
.
,
x
N
)
T
,
Y
=
(
y
1
,
y
2
,
.
.
.
,
y
N
)
T
X = (x_1,x_2,...,x_N)^T,Y=(y_1,y_2,...,y_N)^T
X=(x1,x2,...,xN)T,Y=(y1,y2,...,yN)T
假设X和Y之间存在线性关系，模型的具体形式为
y
^
=
f
(
w
)
=
w
T
x
\hat{y}=f(w) =w^Tx
y^=f(w)=wTx
(a) 最小二乘估计：
我们需要衡量真实值
y
i
y_i
yi与线性回归模型的预测值
w
T
x
i
w^Tx_i
wTxi之间的差距，在这里我们和使用二范数的平方和L(w)来描述这种差距，即：
L
(
w
)
=
∑
i
=
1
N
∣
∣
w
T
x
i
−
y
i
∣
∣
2
2
=
∑
i
=
1
N
(
w
T
x
i
−
y
i
)
2
=
(
w
T
X
T
−
Y
T
)
(
w
T
X
T
−
Y
T
)
T
=
w
T
X
T
X
w
−
2
w
T
X
T
Y
+
Y
Y
T
因
此
，
我
们
需
要
找
到
使
得
L
(
w
)
最
小
时
对
应
的
参
数
w
，
即
：
w
^
=
a
r
g
m
i
n
L
(
w
)
为
了
达
到
求
解
最
小
化
L
(
w
)
问
题
，
我
们
应
用
高
等
数
学
的
知
识
，
使
用
求
导
来
解
决
这
个
问
题
：
∂
L
(
w
)
∂
w
=
2
X
T
X
w
−
2
X
T
Y
=
0
,
因
此
：
w
^
=
(
X
T
X
)
−
1
X
T
Y
L(w) = \sum\limits_{i=1}^{N}||w^Tx_i-y_i||_2^2=\sum\limits_{i=1}^{N}(w^Tx_i-y_i)^2 = (w^TX^T-Y^T)(w^TX^T-Y^T)^T = w^TX^TXw - 2w^TX^TY+YY^T\\ 因此，我们需要找到使得L(w)最小时对应的参数w，即：\\ \hat{w} = argmin\;L(w)\\ 为了达到求解最小化L(w)问题，我们应用高等数学的知识，使用求导来解决这个问题： \\ \frac{\partial L(w)}{\partial w} = 2X^TXw-2X^TY = 0,因此： \\ \hat{w} = (X^TX)^{-1}X^TY
L(w)=i=1∑N∣∣wTxi−yi∣∣22=i=1∑N(wTxi−yi)2=(wTXT−YT)(wTXT−YT)T=wTXTXw−2wTXTY+YYT因此，我们需要找到使得L(w)最小时对应的参数w，即：w^=argminL(w)为了达到求解最小化L(w)问题，我们应用高等数学的知识，使用求导来解决这个问题：∂w∂L(w)=2XTXw−2XTY=0,因此：w^=(XTX)−1XTY
(b) 几何解释：
在线性代数中，我们知道两个向量a和b相互垂直可以得出：
<
a
,
b
>
=
a
.
b
=
a
T
b
=
0
<a,b> = a.b = a^Tb = 0
<a,b>=a.b=aTb=0,而平面X的法向量为Y-Xw，与平面X互相垂直，因此：
X
T
(
Y
−
X
w
)
=
0
X^T(Y-Xw) = 0
XT(Y−Xw)=0，即：
w
=
(
X
T
X
)
−
1
X
T
Y
w = (X^TX)^{-1}X^TY
w=(XTX)−1XTY
© 概率视角：
假设噪声
ϵ
∽
N
(
0
,
σ
2
)
,
y
=
f
(
w
)
+
ϵ
=
w
T
x
+
ϵ
\epsilon \backsim N(0,\sigma^2),y=f(w)+\epsilon=w^Tx+\epsilon
ϵ∽N(0,σ2),y=f(w)+ϵ=wTx+ϵ，因此：
y
∣
x
i
,
w
N
(
w
T
x
,
σ
2
)
y|x_i,w ~ N(w^Tx,\sigma^2)
y∣xi,w N(wTx,σ2)
我们使用极大似然估计MLE对参数w进行估计：
L
(
w
)
=
l
o
g
P
(
Y
∣
X
;
w
)
=
l
o
g
∏
i
=
1
N
P
(
y
i
∣
x
i
;
w
)
=
∑
i
=
1
N
l
o
g
P
(
y
i
∣
x
i
;
w
)
=
∑
i
=
1
N
l
o
g
(
1
2
π
σ
e
x
p
(
−
(
y
i
−
w
T
x
i
)
2
2
σ
2
)
)
=
∑
i
=
1
N
[
l
o
g
(
1
2
π
σ
)
−
1
2
σ
2
(
y
i
−
w
T
x
i
)
2
]
a
r
g
m
a
x
w
L
(
w
)
=
a
r
g
m
i
n
w
[
l
(
w
)
=
∑
i
=
1
N
(
y
i
−
w
T
x
i
)
2
]
因
此
：
线
性
回
归
的
最
小
二
乘
估
计
<
=
=
>
噪
声
ϵ
∽
N
(
0
,
σ
2
)
的
极
大
似
然
估
计
L(w) = log\;P(Y|X;w) = log\;\prod_{i=1}^N P(y_i|x_i;w) = \sum\limits_{i=1}^{N} log\; P(y_i|x_i;w)\\ = \sum\limits_{i=1}^{N}log(\frac{1}{\sqrt{2\pi \sigma}}exp(-\frac{(y_i-w^Tx_i)^2}{2\sigma^2})) = \sum\limits_{i=1}^{N}[log(\frac{1}{\sqrt{2\pi}\sigma})-\frac{1}{2\sigma^2}(y_i-w^Tx_i)^2] \\ argmax_w L(w) = argmin_w[l(w) = \sum\limits_{i = 1}^{N}(y_i-w^Tx_i)^2]\\ 因此：线性回归的最小二乘估计<==>噪声\epsilon\backsim N(0,\sigma^2)的极大似然估计
L(w)=logP(Y∣X;w)=logi=1∏NP(yi∣xi;w)=i=1∑NlogP(yi∣xi;w)=i=1∑Nlog(2πσ
1exp(−2σ2(yi−wTxi)2))=i=1∑N[log(2π
σ1)−2σ21(yi−wTxi)2]argmaxwL(w)=argminw[l(w)=i=1∑N(yi−wTxi)2]因此：线性回归的最小二乘估计<==>噪声ϵ∽N(0,σ2)的极大似然估计
下面，我们使用sklearn的线性回归实例来演示：
https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html#sklearn.linear_model.LinearRegression
from sklearn import linear_model # 引入线性回归方法
lin_reg = linear_model.LinearRegression() # 创建线性回归的类
lin_reg.fit(X,y) # 输入特征X和因变量y进行训练
print("模型系数：",lin_reg.coef_) # 输出模型的系数
print("模型得分：",lin_reg.score(X,y)) # 输出模型的决定系数R^2
# 此处spss ，Excel 都可以进行相关的操作
模型系数： [-1.07170557e-01 4.63952195e-02 2.08602395e-02 2.68856140e+00
-1.77957587e+01 3.80475246e+00 7.51061703e-04 -1.47575880e+00
3.05655038e-01 -1.23293463e-02 -9.53463555e-01 9.39251272e-03
-5.25466633e-01]
模型得分： 0.7406077428649427
(b) 广义可加模型(GAM)：
广义可加模型GAM实际上是线性模型推广至非线性模型的一个框架，在这个框架中，每一个变量都用一个非线性函数来代替，但是模型本身保持整体可加性。GAM模型不仅仅可以用在线性回归的推广，还可以将线性分类模型进行推广。具体的推广形式是：
标准的线性回归模型：
y
i
=
w
0
+
w
1
x
i
1
+
.
.
.
+
w
p
x
i
p
+
ϵ
i
y_i = w_0 + w_1x_{i1} +...+w_px_{ip} + \epsilon_i
yi=w0+w1xi1+...+wpxip+ϵi
GAM模型框架：
y
i
=
w
0
+
∑
j
=
1
p
f
j
(
x
i
j
)
+
ϵ
i
y_i = w_0 + \sum\limits_{j=1}^{p}f_{j}(x_{ij}) + \epsilon_i
yi=w0+j=1∑pfj(xij)+ϵi
GAM模型的优点与不足：
- 优点：简单容易操作，能够很自然地推广线性回归模型至非线性模型，使得模型的预测精度有所上升；由于模型本身是可加的，因此GAM还是能像线性回归模型一样把其他因素控制不变的情况下单独对某个变量进行推断，极大地保留了线性回归的易于推断的性质。
- 缺点：GAM模型会经常忽略一些有意义的交互作用，比如某两个特征共同影响因变量，不过GAM还是能像线性回归一样加入交互项
x
(
i
)
×
x
(
j
)
x^{(i)} \times x^{(j)}
x(i)×x(j)的形式进行建模；但是GAM模型本质上还是一个可加模型，如果我们能摆脱可加性模型形式，可能还会提升模型预测精度，详情请看后面的算法。
(a) 多项式回归实例介绍：
https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PolynomialFeatures.html?highlight=poly#sklearn.preprocessing.PolynomialFeatures
sklearn.preprocessing.PolynomialFeatures(degree=2, *, interaction_only=False, include_bias=True, order=‘C’):
from sklearn.preprocessing import PolynomialFeatures
# 这个类可以进行特征的构造，构造的方式就是特征与特征相乘（自己与自己，自己与其他人），这种方式叫做使用多项式的方式。
# 例如：有 a、b 两个特征，那么它的 2 次多项式的次数为 [1,a,b,a2,ab,b2]
X_arr = np.arange(6).reshape(3, 2)
print("原始X为：\n",X_arr)
# 第一种
poly = PolynomialFeatures(2)
print("2次转化X：\n",poly.fit_transform(X_arr))
# 第二种
poly = PolynomialFeatures(interaction_only=True) #指定为 True，不会有特征自己和自己结合的项，组合的特征中没有 a2 和 b2
print("2次转化X：\n",poly.fit_transform(X_arr))
# 两者比较 ，第一种，interaction_only=False,（自己与自己/其他人）
# 第二种，interaction_only=True (自己与其他人,不会自己和自己结合)
原始X为：
[[0 1]
[2 3]
[4 5]]
2次转化X：
[[ 1. 0. 1. 0. 0. 1.]
[ 1. 2. 3. 4. 6. 9.]
[ 1. 4. 5. 16. 20. 25.]]
2次转化X：
[[ 1. 0. 1. 0.]
[ 1. 2. 3. 6.]
[ 1. 4. 5. 20.]]
(b) GAM模型实例介绍：
安装pygam：pip install pygam
https://github.com/dswah/pyGAM/blob/master/doc/source/notebooks/quick_start.ipynb
# 建立模型
from pygam import LinearGAM
gam = LinearGAM().fit(boston_data[boston.feature_names], y)
# 查看结果。summary()提供了很多统计量，比如AIC，UBRE，修正R^2。
gam.summary()
LinearGAM
=============================================== ==========================================================
Distribution: NormalDist Effective DoF: 103.2768
Link Function: IdentityLink Log Likelihood: -1589.695
Number of Samples: 506 AIC: 3387.9436
AICc: 3442.7341
GCV: 13.7688
Scale: 8.8256
Pseudo R-Squared: 0.9168
==========================================================================================================
Feature Function Lambda Rank EDoF P > x Sig. Code
================================= ==================== ============ ============ ============ ============
s(0) [0.6] 20 11.4 2.02e-11 ***
s(1) [0.6] 20 12.6 8.10e-02 .
s(2) [0.6] 20 13.4 2.53e-03 **
s(3) [0.6] 20 3.7 2.65e-01
s(4) [0.6] 20 11.4 1.11e-16 ***
s(5) [0.6] 20 10.2 1.11e-16 ***
s(6) [0.6] 20 10.4 8.29e-01
s(7) [0.6] 20 8.5 6.66e-16 ***
s(8) [0.6] 20 3.6 6.03e-03 **
s(9) [0.6] 20 3.6 1.25e-09 ***
s(10) [0.6] 20 1.9 3.33e-03 **
s(11) [0.6] 20 6.3 5.39e-02 .
s(12) [0.6] 20 6.5 1.11e-16 ***
intercept 1 0.0 9.89e-14 ***
==========================================================================================================
Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
建立回归树的过程大致可以分为以下两步：
- 将自变量的特征空间(即
x
(
1
)
,
x
(
2
)
,
x
(
3
)
,
.
.
.
,
x
(
p
)
x^{(1)},x^{(2)},x^{(3)},...,x^{(p)}
x(1),x(2),x(3),...,x(p))的可能取值构成的集合分割成J个互不重叠的区域
R
1
,
R
2
,
.
.
.
,
R
j
R_1,R_2,...,R_j
R1,R2,...,Rj。
- 对落入区域
R
j
R_j
Rj的每个观测值作相同的预测，预测值等于
R
j
R_j
Rj上训练集的因变量的简单算术平均。
具体来说，就是：
a. 选择最优切分特征j以及该特征上的最优点s：
遍历特征j以及固定j后遍历切分点s，选择使得下式最小的(j,s)
m
i
n
j
,
s
[
m
i
n
c
1
∑
x
i
∈
R
1
(
j
,
s
)
(
y
i
−
c
1
)
2
+
m
i
n
c
2
∑
x
i
∈
R
2
(
j
,
s
)
(
y
i
−
c
2
)
2
]
min_{j,s}[min_{c_1}\sum\limits_{x_i\in R_1(j,s)}(y_i-c_1)^2 + min_{c_2}\sum\limits_{x_i\in R_2(j,s)}(y_i-c_2)^2 ]
minj,s[minc1xi∈R1(j,s)∑(yi−c1)2+minc2xi∈R2(j,s)∑(yi−c2)2]
b. 按照(j,s)分裂特征空间：
R
1
(
j
,
s
)
=
{
x
∣
x
j
≤
s
}
和
R
2
(
j
,
s
)
=
{
x
∣
x
j
>
s
}
,
c
^
m
=
1
N
m
∑
x
∈
R
m
(
j
,
s
)
y
i
,
m
=
1
,
2
R_1(j,s) = \{x|x^{j} \le s \}和R_2(j,s) = \{x|x^{j} > s \},\hat{c}_m = \frac{1}{N_m}\sum\limits_{x \in R_m(j,s)}y_i,\;m=1,2
R1(j,s)={x∣xj≤s}和R2(j,s)={x∣xj>s},c^m=Nm1x∈Rm(j,s)∑yi,m=1,2
c. 继续调用步骤1，2直到满足停止条件，就是每个区域的样本数小于等于5。
d. 将特征空间划分为J个不同的区域，生成回归树：
f
(
x
)
=
∑
m
=
1
J
c
^
m
I
(
x
∈
R
m
)
f(x) = \sum\limits_{m=1}^{J}\hat{c}_mI(x \in R_m)
f(x)=m=1∑Jc^mI(x∈Rm)
如以下生成的关于运动员在棒球大联盟数据的回归树：
回归树与线性模型的比较：
线性模型的模型形式与树模型的模型形式有着本质的区别，具体而言，线性回归对模型形式做了如下假定：
f
(
x
)
=
w
0
+
∑
j
=
1
p
w
j
x
(
j
)
f(x) = w_0 + \sum\limits_{j=1}^{p}w_jx^{(j)}
f(x)=w0+j=1∑pwjx(j)，而回归树则是
f
(
x
)
=
∑
m
=
1
J
c
^
m
I
(
x
∈
R
m
)
f(x) = \sum\limits_{m=1}^{J}\hat{c}_mI(x \in R_m)
f(x)=m=1∑Jc^mI(x∈Rm)。那问题来了，哪种模型更优呢？这个要视具体情况而言，如果特征变量与因变量的关系能很好的用线性关系来表达，那么线性回归通常有着不错的预测效果，拟合效果则优于不能揭示线性结构的回归树。反之，如果特征变量与因变量的关系呈现高度复杂的非线性，那么树方法比传统方法更优。
树模型的优缺点：
- 树模型的解释性强，在解释性方面可能比线性回归还要方便。
- 树模型更接近人的决策方式。
- 树模型可以用图来表示，非专业人士也可以轻松解读。
- 树模型可以直接做定性的特征而不需要像线性回归一样哑元化。
- 树模型能很好处理缺失值和异常值，对异常值不敏感，但是这个对线性模型来说却是致命的。
- 树模型的预测准确性一般无法达到其他回归模型的水平，但是改进的方法很多。
sklearn使用回归树的实例：
https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeRegressor.html?highlight=tree#sklearn.tree.DecisionTreeRegressor
sklearn.tree.DecisionTreeRegressor(*, criterion=‘mse’, splitter=‘best’, max_depth=None, min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0, max_features=None, random_state=None, max_leaf_nodes=None, min_impurity_decrease=0.0, min_impurity_split=None, presort=‘deprecated’, ccp_alpha=0.0）
from sklearn.tree import DecisionTreeRegressor
reg_tree = DecisionTreeRegressor(criterion = "mse",min_samples_leaf = 5) # 使用均方误差
reg_tree.fit(X,y)
reg_tree.score(X,y) #衡量模型好坏程度
0.9369572966959968
支持向量机回归(SVR)
在介绍支持向量回归SVR之前，我们先来了解下约束优化的相关知识：
KKT条件(最优解的一阶必要条件)
因为KKT条件是最优化的相关内容，在本次开源学习中并不是重点，因此在这里我用一个更加简单的例子说明KKT条件，严格的证明请参见凸优化相关书籍。
在这个例子中，我们考虑：(
x
∗
x^*
x∗为我们的最优解)
m
i
n
f
(
x
)
s
.
t
.
g
1
(
x
)
≤
0
,
x
∈
R
n
g
2
(
x
)
≤
0
g
3
(
x
)
≤
0
minf(x)\\ s.t.\;g_1(x) \le 0,\;x \in R^n\\ \;\;\;g_2(x) \le 0\\ \;\;\;g_3(x) \le 0
minf(x)s.t.g1(x)≤0,x∈Rng2(x)≤0g3(x)≤0
我们可以看到：
−
∇
f
(
x
∗
)
-\nabla f(x^*)
−∇f(x∗)可以由
∇
g
1
(
x
∗
)
\nabla g_1(x^*)
∇g1(x∗)与
∇
g
2
(
x
∗
)
\nabla g_2(x^*)
∇g2(x∗)线性表出，因此有：
−
∇
f
(
x
∗
)
=
λ
1
∇
g
1
(
x
∗
)
+
λ
2
∇
g
2
(
x
∗
)
-\nabla f(x^*) = \lambda_1 \nabla g_1(x^*) + \lambda_2 \nabla g_2(x^*)
−∇f(x∗)=λ1∇g1(x∗)+λ2∇g2(x∗)，其中
λ
1
,
λ
2
≥
0
\lambda_1,\lambda_2 \ge 0
λ1,λ2≥0，即：
∇
f
(
x
∗
)
+
λ
1
∇
g
1
(
x
∗
)
+
λ
2
∇
g
2
(
x
∗
)
=
0
,
其
中
λ
1
,
λ
2
≥
0
\nabla f(x^*) + \lambda_1 \nabla g_1(x^*) + \lambda_2 \nabla g_2(x^*) = 0,\;\;\;其中\lambda_1,\lambda_2 \ge 0
∇f(x∗)+λ1∇g1(x∗)+λ2∇g2(x∗)=0,其中λ1,λ2≥0
我们把没有起作用的约束
g
3
(
x
)
g_3(x)
g3(x)也放到式子里面去，目的也就是为了书写方便，即要求：
∇
f
(
x
∗
)
+
λ
1
∇
g
1
(
x
∗
)
+
λ
2
∇
g
2
(
x
∗
)
+
λ
3
∇
g
3
(
x
∗
)
=
0
,
其
中
λ
1
,
λ
2
≥
0
,
λ
3
=
0
\nabla f(x^*) + \lambda_1 \nabla g_1(x^*) + \lambda_2 \nabla g_2(x^*) + \lambda_3 \nabla g_3(x^*)= 0,\;\;\;其中\lambda_1,\lambda_2 \ge 0,\lambda_3 = 0
∇f(x∗)+λ1∇g1(x∗)+λ2∇g2(x∗)+λ3∇g3(x∗)=0,其中λ1,λ2≥0,λ3=0
由于点
x
∗
x^*
x∗位于方程
g
1
(
x
)
=
0
g_1(x)=0
g1(x)=0与
g
2
(
x
)
=
0
g_2(x)=0
g2(x)=0上，因此：
λ
1
g
1
(
x
∗
)
=
0
,
λ
2
g
2
(
x
∗
)
=
0
,
λ
3
g
3
(
x
∗
)
=
0
\lambda_1 g_1(x^*) = 0,\lambda_2 g_2(x^*) = 0 , \lambda_3 g_3(x^*)= 0
λ1g1(x∗)=0,λ2g2(x∗)=0,λ3g3(x∗)=0
因此，KKT条件就是：假设
x
∗
x^*
x∗为最优化问题§的局部最优解，且
x
∗
x^*
x∗ 在某个适当的条件下 ,有：
∇
f
(
x
∗
)
+
∑
i
=
1
m
λ
i
∇
g
(
x
∗
)
+
∑
j
=
1
l
μ
j
∇
h
j
(
x
∗
)
=
0
(
对
偶
条
件
)
λ
i
≥
0
,
i
=
1
,
2
,
.
.
.
,
m
(
对
偶
条
件
)
g
i
(
x
∗
)
≤
0
(
原
问
题
条
件
)
h
j
(
x
∗
)
=
0
(
原
问
题
条
件
)
λ
i
g
(
x
∗
)
=
0
(
互
补
松
弛
定
理
)
\nabla f(x^*) + \sum\limits_{i=1}^{m}\lambda_i \nabla g(x^*) + \sum\limits_{j=1}^{l}\mu_j \nabla h_j(x^*) = 0(对偶条件)\\ \lambda_i \ge 0,\;i = 1,2,...,m(对偶条件)\\ g_i(x^*) \le 0(原问题条件)\\ h_j(x^*) = 0(原问题条件)\\ \lambda_i g(x^*) = 0(互补松弛定理)
∇f(x∗)+i=1∑mλi∇g(x∗)+j=1∑lμj∇hj(x∗)=0(对偶条件)λi≥0,i=1,2,...,m(对偶条件)gi(x∗)≤0(原问题条件)hj(x∗)=0(原问题条件)λig(x∗)=0(互补松弛定理)
对偶理论：
为什么要引入对偶问题呢？是因为原问题与对偶问题就像是一个问题两个角度去看，如利润最大与成本最低等。有时侯原问题上难以解决，但是在对偶问题上就会变得很简单。再者，任何一个原问题在变成对偶问题后都会变成一个凸优化的问题，这点我们后面会有介绍。下面我们来引入对偶问题：
首先，我们的原问题§是：
m
i
n
f
(
x
)
s
.
t
.
g
i
(
x
)
≤
0
,
i
=
1
,
2
,
.
.
.
,
m
h
j
(
x
)
=
0
,
j
=
1
,
2
,
.
.
.
,
l
min f(x) \\ s.t.\;\;\;g_i(x) \le 0,\; i=1,2,...,m\\ \;\;\;\;\; h_j(x) = 0,\; j=1,2,...,l
minf(x)s.t.gi(x)≤0,i=1,2,...,mhj(x)=0,j=1,2,...,l
引入拉格朗日函数：
L
(
x
,
λ
,
μ
)
=
f
(
x
)
+
∑
i
=
1
m
λ
i
g
i
(
x
)
+
∑
j
=
1
l
μ
j
h
j
(
x
)
L(x,\lambda,\mu) = f(x) + \sum\limits_{i=1}^{m}\lambda_i g_i(x) + \sum\limits_{j=1}^{l}\mu_j h_j(x)
L(x,λ,μ)=f(x)+i=1∑mλigi(x)+j=1∑lμjhj(x)
拉格朗日对偶函数：
d
(
λ
,
μ
)
=
m
i
n
x
∈
X
{
f
(
x
)
+
∑
i
=
1
m
λ
i
g
i
(
x
)
+
∑
j
=
1
l
μ
j
h
j
(
x
)
}
,
其
中
X
为
满
足
条
件
的
x
变
量
≤
m
i
n
x
∈
S
{
f
(
x
)
+
∑
i
=
1
m
λ
i
g
i
(
x
)
+
∑
j
=
1
l
μ
j
h
j
(
x
)
}
,
由
于
g
i
(
x
)
≤
0
,
h
j
(
x
)
=
0
,
λ
i
≥
0
,
其
中
S
为
可
行
域
≤
m
i
n
x
∈
S
{
f
(
x
)
}
d(\lambda,\mu) = min_{x\in X}\{ f(x) + \sum\limits_{i=1}^{m}\lambda_i g_i(x) + \sum\limits_{j=1}^{l}\mu_j h_j(x)\} ,其中X为满足条件的x变量\\ \le min_{x\in S}\{ f(x) + \sum\limits_{i=1}^{m}\lambda_i g_i(x) + \sum\limits_{j=1}^{l}\mu_j h_j(x) \},由于g_i(x) \le 0,h_j(x) = 0,\lambda_i \ge 0 ,其中S为可行域\\ \le min_{x\in S}\{f(x) \}
d(λ,μ)=minx∈X{f(x)+i=1∑mλigi(x)+j=1∑lμjhj(x)},其中X为满足条件的x变量≤minx∈S{f(x)+i=1∑mλigi(x)+j=1∑lμjhj(x)},由于gi(x)≤0,hj(x)=0,λi≥0,其中S为可行域≤minx∈S{f(x)}
因此：拉格朗日对偶函数
d
(
λ
,
μ
)
d(\lambda,\mu)
d(λ,μ)是原问题最优解的函数值
p
∗
p^*
p∗的下界，即每个不同的
λ
\lambda
λ与
μ
\mu
μ确定的
d
(
λ
,
μ
)
d(\lambda,\mu)
d(λ,μ)都是
p
∗
p^*
p∗的下界，但是我们希望下界越大越好，因为越大就更能接近真实的
p
∗
p^*
p∗。因此：
拉格朗日对偶问题(D)转化为：
m
a
x
λ
,
μ
d
(
λ
,
μ
)
s
.
t
.
λ
i
≥
0
,
i
=
1
,
2
,
.
.
.
,
m
也
就
是
：
m
a
x
λ
≥
0
,
μ
m
i
n
x
∈
S
L
(
x
,
λ
,
μ
)
max_{\lambda,\mu}d(\lambda,\mu)\\ s.t. \lambda_i \ge 0,i = 1,2,...,m\\ 也就是：\\ max_{\lambda \ge 0,\mu}\;min_{x \in S} L(x,\lambda,\mu)
maxλ,μd(λ,μ)s.t.λi≥0,i=1,2,...,m也就是：maxλ≥0,μminx∈SL(x,λ,μ)
我们可以观察到，对偶问题是关于
λ
\lambda
λ和
μ
\mu
μ的线性函数，因此对偶问题是一个凸优化问题，凸优化问题在最优化理论较为简单。
弱对偶定理：对偶问题(D)的最优解
D
∗
D^*
D∗一定小于原问题最优解
P
∗
P^*
P∗，这点在刚刚的讨论得到了充分的证明，一定成立。
强对偶定理：对偶问题(D)的最优解
D
∗
D^*
D∗在一定的条件下等于原问题最优解
P
∗
P^*
P∗，条件非常多样化且不是唯一的，也就是说这是个开放性的问题，在这里我给出一个最简单的条件，即：
f
(
x
)
f(x)
f(x)与
g
i
(
x
)
g_i(x)
gi(x)为凸函数，
h
j
(
x
)
h_j(x)
hj(x)为线性函数，X是凸集，
x
∗
x^*
x∗满足KKT条件，那么
D
∗
=
P
∗
D^* = P^*
D∗=P∗。
在线性回归的理论中，每个样本点都要计算平方损失，但是SVR却是不一样的。SVR认为：落在
f
(
x
)
f(x)
f(x)的
ϵ
\epsilon
ϵ邻域空间中的样本点不需要计算损失，这些都是预测正确的，其余的落在
ϵ
\epsilon
ϵ邻域空间以外的样本才需要计算损失，因此：
m
i
n
w
,
b
,
ξ
i
,
ξ
^
i
1
2
∣
∣
w
∣
∣
2
+
C
∑
i
=
1
N
(
ξ
i
,
ξ
^
i
)
s
.
t
.
f
(
x
i
)
−
y
i
≤
ϵ
+
ξ
i
y
i
−
f
(
x
i
)
≤
ϵ
+
ξ
^
i
ξ
i
,
ξ
^
i
≤
0
,
i
=
1
,
2
,
.
.
.
,
N
min_{w,b,\xi_i,\hat{\xi}_i} \frac{1}{2}||w||^2 +C \sum\limits_{i=1}^{N}(\xi_i,\hat{\xi}_i)\\ s.t.\;\;\; f(x_i) - y_i \le \epsilon + \xi_i\\ \;\;\;\;\;y_i - f(x_i) \le \epsilon +\hat{\xi}_i\\ \;\;\;\;\; \xi_i,\hat{\xi}_i \le 0,i = 1,2,...,N
minw,b,ξi,ξ^i21∣∣w∣∣2+Ci=1∑N(ξi,ξ^i)s.t.f(xi)−yi≤ϵ+ξiyi−f(xi)≤ϵ+ξ^iξi,ξ^i≤0,i=1,2,...,N
引入拉格朗日函数：
L
(
w
,
b
,
α
,
α
^
,
ξ
,
ξ
,
μ
,
μ
^
)
=
1
2
∥
w
∥
2
+
C
∑
i
=
1
N
(
ξ
i
+
ξ
^
i
)
−
∑
i
=
1
N
ξ
i
μ
i
−
∑
i
=
1
N
ξ
^
i
μ
^
i
+
∑
i
=
1
N
α
i
(
f
(
x
i
)
−
y
i
−
ϵ
−
ξ
i
)
+
∑
i
=
1
N
α
^
i
(
y
i
−
f
(
x
i
)
−
ϵ
−
ξ
^
i
)
\begin{array}{l} L(w, b, \alpha, \hat{\alpha}, \xi, \xi, \mu, \hat{\mu}) \\ \quad=\frac{1}{2}\|w\|^{2}+C \sum_{i=1}^{N}\left(\xi_{i}+\widehat{\xi}_{i}\right)-\sum_{i=1}^{N} \xi_{i} \mu_{i}-\sum_{i=1}^{N} \widehat{\xi}_{i} \widehat{\mu}_{i} \\ \quad+\sum_{i=1}^{N} \alpha_{i}\left(f\left(x_{i}\right)-y_{i}-\epsilon-\xi_{i}\right)+\sum_{i=1}^{N} \widehat{\alpha}_{i}\left(y_{i}-f\left(x_{i}\right)-\epsilon-\widehat{\xi}_{i}\right) \end{array}
L(w,b,α,α^,ξ,ξ,μ,μ^)=21∥w∥2+C∑i=1N(ξi+ξ
i)−∑i=1Nξiμi−∑i=1Nξ
iμ
i+∑i=1Nαi(f(xi)−yi−ϵ−ξi)+∑i=1Nα
i(yi−f(xi)−ϵ−ξ
i)
再令
L
(
w
,
b
,
α
,
α
^
,
ξ
,
ξ
,
μ
,
μ
^
)
L(w, b, \alpha, \hat{\alpha}, \xi, \xi, \mu, \hat{\mu})
L(w,b,α,α^,ξ,ξ,μ,μ^)对
w
,
b
,
ξ
,
ξ
^
w,b,\xi,\hat{\xi}
w,b,ξ,ξ^求偏导等于0，得：
w
=
∑
i
=
1
N
(
α
^
i
−
α
i
)
x
i
w=\sum_{i=1}^{N}\left(\widehat{\alpha}_{i}-\alpha_{i}\right) x_{i}
w=∑i=1N(α
i−αi)xi。
上述过程中需满足KKT条件，即要求：
{
α
i
(
f
(
x
i
)
−
y
i
−
ϵ
−
ξ
i
)
=
0
α
i
^
(
y
i
−
f
(
x
i
)
−
ϵ
−
ξ
^
i
)
=
0
α
i
α
^
i
=
0
,
ξ
i
ξ
^
i
=
0
(
C
−
α
i
)
ξ
i
=
0
,
(
C
−
α
^
i
)
ξ
^
i
=
0
\left\{\begin{array}{c} \alpha_{i}\left(f\left(x_{i}\right)-y_{i}-\epsilon-\xi_{i}\right)=0 \\ \hat{\alpha_{i}}\left(y_{i}-f\left(x_{i}\right)-\epsilon-\hat{\xi}_{i}\right)=0 \\ \alpha_{i} \widehat{\alpha}_{i}=0, \xi_{i} \hat{\xi}_{i}=0 \\ \left(C-\alpha_{i}\right) \xi_{i}=0,\left(C-\widehat{\alpha}_{i}\right) \hat{\xi}_{i}=0 \end{array}\right.
⎩⎪⎪⎪⎨⎪⎪⎪⎧αi(f(xi)−yi−ϵ−ξi)=0αi^(yi−f(xi)−ϵ−ξ^i)=0αiα
i=0,ξiξ^i=0(C−αi)ξi=0,(C−α
i)ξ^i=0
SVR的解形如:
f
(
x
)
=
∑
i
=
1
N
(
α
^
i
−
α
i
)
x
i
T
x
+
b
f(x)=\sum_{i=1}^{N}\left(\widehat{\alpha}_{i}-\alpha_{i}\right) x_{i}^{T} x+b
f(x)=∑i=1N(α
i−αi)xiTx+b
sklearn中使用SVR实例：
sklearn.svm.SVR(*, kernel=‘rbf’, degree=3, gamma=‘scale’, coef0=0.0, tol=0.001, C=1.0, epsilon=0.1, shrinking=True, cache_size=200, verbose=False, max_iter=-1)
https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVR.html?highlight=svr#sklearn.svm.SVR
# 支持向量机回归
from sklearn.svm import SVR
from sklearn.preprocessing import StandardScaler # 标准化数据
from sklearn.pipeline import make_pipeline # 使用管道，把预处理和模型形成一个流程
reg_svr = make_pipeline(StandardScaler(), SVR(C=1.0, epsilon=0.2))
reg_svr.fit(X, y)
reg_svr.score(X,y)
0.7024053100500254