最优化方法Python计算：有约束优化应用——线性Lasso回归预测器

实际应用中，特征维度$n$通常远大于样本容量$m$（$n\ll m$），这种高维小样本场景下特征数据可能含有对标签数据$y_i$的取值不必要的成分，甚至是噪声。此时，我们希望回归模型中的优化问题
$$\{\begin{array}{l}\text{minimize}\frac{1}{2}\Vert \mathbf{F}(\mathbf{w};\mathbf{X})-\mathbf{y}{\Vert }_2^2\\ \text{s.t.}\mathbf{w}\in {\text{R}}^p\end{array}$$
其中，$p\ge n$，的最优解${\mathbf{w}}_0$是稀疏的：${\mathbf{w}}_0$中对应冗余特征项的元素为0，即过滤掉那些噪声的影响。为此，需要在上述无约束优化问题中添加一个约束条件：$\mathbf{w}$所含非零元素最少（但并非全为零）。即
$$\{\begin{array}{l}\text{minimize}\frac{1}{2}\Vert \mathbf{F}(\mathbf{w};\mathbf{X})-\mathbf{y}{\Vert }_2^2\\ \text{s.t.}\mathbf{w}\text{非零元素个数尽可能少},\mathbf{w}\in {\text{R}}^p\end{array}$$
用度量向量“稀疏性”的伪范数$L_0$，
$$\Vert \mathbf{x}{\Vert }_0=\sum _{i=1}^n\left\{\begin{array}{ll}1,& x_i\ne 0\\ 0,& x_i=0\end{array}\right\}.$$
则上述优化问题等价第表示为：
$$\{\begin{array}{l}\text{minimize}\Vert \mathbf{w}{\Vert }_0\\ \text{s.t.}\mathbf{F}(\mathbf{w};\mathbf{X})=\mathbf{y}\end{array}.$$
由于$\Vert \cdot {\Vert }_0$非连续函数，故常用的优化问题解法很难解决上述问题。可以证明，$L_1$范数$\Vert \cdot {\Vert }_1$也可用于度量向量的稀疏性。即可用
$$\{\begin{array}{l}\text{minimize}\Vert \mathbf{w}{\Vert }_1\\ \text{s.t.}\mathbf{F}(\mathbf{w};\mathbf{X})=\mathbf{y}\end{array}$$
得到回归模型的稀疏解。这是一个有等式约束的优化问题，可以通过将约束条件作为罚项，求解无约束优化问题
$$\{\begin{array}{l}\text{minimize}\Vert \mathbf{w}{\Vert }_1+\sigma \frac{1}{2}\Vert \mathbf{F}(\mathbf{w};\mathbf{X})-\mathbf{y}{\Vert }_2^2\\ \text{s.t.}\mathbf{w}\in {\text{R}}^p\end{array}$$
完成。这称为Lasso回归问题，若令$\alpha =\frac{1}{\sigma }$，上述问题的目标函数等价于
$$\text{obj}(\mathbf{w},\alpha )=\frac{1}{2}\Vert \mathbf{F}(\mathbf{w};\mathbf{X})-\mathbf{y}{\Vert }_2^2+\alpha \Vert \mathbf{w}{\Vert }_1.$$
$\sigma \to \infty$当且仅当$\alpha \to 0$。而当$\alpha =0$时，$\arg \underset{\mathbf{w}\in {\text{R}}^p}{\min }\text{obj}(\mathbf{w},0)=\arg \underset{\mathbf{w}\in {\text{R}}^p}{\min }\frac{1}{2}\Vert F(\mathbf{w};\mathbf{A})-\mathbf{y}{\Vert }_2^2$。即Lasso回归退化为原回归问题，这显然不符合预期。因此，Lasso回归需要寻求合适的$0<\alpha <+\infty$，使得${\mathbf{w}}_0=\arg \underset{\mathbf{w}\in {\text{R}}^p}{\min }\text{obj}(\mathbf{w},\alpha )$为稀疏解。以式$\frac{1}{2}\Vert \mathbf{F}(\mathbf{w};\mathbf{X})-\mathbf{y}{\Vert }_2^2+\alpha \Vert \mathbf{w}{\Vert }_1$为目标函数的优化问题也称为$L_1$范数正则化模型，其中正则化参数$\alpha$用来控制最优解${\mathbf{w}}_0$的稀疏程度。当回归函数$\mathbf{F}(\mathbf{w};\mathbf{X})$为线性函数$\mathbf{Aw}$，其中$\mathbf{A}=({\mathbf{X}}^{\top },1)$，此时$p=n+1$，称为线性Lasso回归。而当$\mathbf{F}(\mathbf{w};\mathbf{X})=\frac{1}{1+e^{-\mathbf{Aw}}}$时，称为逻辑Lasso回归。
下列代码实现线性Lasso回归模型。

import numpy as np										#导入numpy
class LineLassoModel(LineModel):						#线性lasso模型def obj(self,w):h = lambda x: np.linalg.norm(x, ord = 1)return super().obj(w) + self.alpha * h(w)
class LassoRegressor(Regression, LineLassoModel):		#线性Lasso回归预测器def __init__(self, alpha = 1.0):self.alpha = alphaself.tagVal = identity

程序中，

• 第2~5行将LineLassoModel定义为LineModel（见博文《最优化方法Python计算：无约束优化应用——线性回归模型》）的子类，其中的第3~5行重载目标函数obj。第4行定义$L_1$范数函数$h(\mathbf{x})=\Vert \mathbf{x}{\Vert }_1$。注意，此处调用Numpy的norm函数，传递ord参数为1，计算一阶范数。
• 第6~9行将线性Lasso回归预测器类LassoRegressor定义为Regression（见博文《最优化方法Python计算：无约束优化应用——线性回归预测器》）与LineLassoModel的子类。其中，第7~9行定义其构造函数__init__，第8行设置正则参数alpha，缺省值为1.0。第9行将标签值函数tagVal设置为恒等函数identity，已适用于回归预测。
  综合案例
  文件forestfires.csv（来自UC Irvine Machine Learning Repository）包含517条葡萄牙东北部地区蒙特西尼奥公园森林火灾的烧毁面积的数据记录。

X	Y	month	day	FFMC	DMC	DC	ISI	temp	RH	wind	rain	area
7	5	mar	fri	86.2	26.2	94.3	5.1	8.2	51	6.7	0	0
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
1	4	sep	tue	91	129.5	692.6	7	21.7	38	2.2	0	0.43
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
7	6	jul	tue	93.1	180.4	430.8	11	26.9	28	5.4	0	86.45
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
4	4	aug	wed	95.1	141.3	605.8	17.7	28.7	33	4	0	0
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
6	3	nov	tue	79.5	3	106.7	1.1	11.8	31	4.5	0	0

其中：

• X - 蒙特西尼奥公园地图内的x轴空间坐标：1到9
• Y - 蒙特西尼奥公园地图内的y轴空间坐标：2到9
• month - 月份-一年中的月份：“jan”到“dec”
• day - 一周中的某一天：“mon”到“sun”
• FFMC - 来自快速预警系统的FFMC指数：18.7至96.20
• DMC - 来自快速预警系统的DMC指数：1.1至291.3
• DC - 来自快速预警系统的DC指数：7.9至860.6
• ISI - 来自快速预警系统的ISI指数：0.0至56.10
• temp - 摄氏度：2.2至33.30
• RH - 相对湿度%{}：15.0至100
• wind - 风速（公里/小时）：0.40至9.40
• rain - 室外雨，以毫米/平方米为单位：0.0至6.4
• area - 森林被烧毁面积（以公顷为单位）：0.00至1090.84
  以森林烧毁面积area作为标签，其余12个数据项作为样本特征，我们用线性回归来创建、训练和测试模型。下列代码完成计算

import numpy as np
data = np.loadtxt('forestfires.csv', delimiter=',',		#读取数据文件dtype=str, skiprows=1)				#跳过表头
X = np.array(data)										#转换为数组
Y = X[:, 12].astype(float)								#读取标签数据
X=np.delete(X, [12], axis = 1)
X1 = X[:, 2]											#读取月份列
X1 = np.array([0 if x == 'jan' else						#串转换为数值1 if x == 'feb' else2 if x == 'mar' else3 if x == 'apr' else4 if x == 'may' else5 if x == 'jun' else6 if x == 'jul' else7 if x == 'aug' else8 if x == 'sept' else9 if x == 'oct' else10 if x == 'nov' else11 for x in X1])
X[:, 2] = X1											#修改月份列
X1 = X[:, 3]											#读取星期列
X1 = np.array([0 if x == 'mon' else						#串转换为数值1 if x == 'tue' else2 if x == 'wed' else3 if x == 'thu' else4 if x == 'fri' else5 if x == 'sat' else6 for x in X1])
X[:, 3] = X1											#修改季节列
X = X.astype(float)										#所有数据转换为实数型
m, n = X.shape											#读取样本个数
print('共有%d个数据样本'%m)
a = np.arange(m)										#下标集
print('随机抽取%d个样本作为训练数据。'%(m // 2))
np.random.seed(7503)									#随机种子
train = np.random.choice(a, m // 2, replace = False)	#随机取得一半作为训练集下标
test = np.setdiff1d(a, train)							#测试集下标
print('LINE REGRESSION:')
forestfires = LinearRegressor()							#构造线性回归模型
forestfires.fit(X[train], Y[train])						#训练模型
coef, inte=forestfires.coef_inte()
print('系数：', coef)
print('截距：%.4f'%inte)
score=forestfires.score(X[test], Y[test])				#测试模型
print('对其余%d个样本数据测试，评估分：%.4f'%(m - m // 2, score))

程序的2~4行读取文件，转换为数组X。第5~6行分离出标签数据Y。第7~20行将month列的字串行转换为数值。第21~29行将day列的字串转换为数值。第30行将特征数据转换为浮点型。第31~37行从数据集的下标中随机选取一半作为训练集下标train，其余的作为测试集下标test。第39行创建线性回归模型LinearRegressor（见《最优化方法Python计算：无约束优化应用——线性回归预测器》）类对象forestfires。第40行调用forestfires的fit函数用X[train]和Y[train]训练forestfires。第44行调用forestfires的score函数用X[test]和Y[test]测试训练好的forestfires。运行程序，输出

共有517个数据样本
随机抽取258个样本作为训练数据。
LINE REGRESSION:
训练中...，稍候
20次迭代后完成训练。
系数： [ 1.2764 -0.2223  4.3442  1.5441  0.2316  0.1098 -0.0589 -1.2177  1.0874-0.2069  1.4681 -9.0991]
截距：-36.2409
对其余259个样本数据测试，评估分：0.0248

意为用线性回归模型计算本问题，经20次迭代训练所得模型的测试评分为0.0248>0，属正常回归。接下来，我们对相同的训练集和测试集训练和测试一个Lasso回归模型。

……
print('LASSO REGRESSION:')
forestfires=LassoRegressor(alpha=0.035)		#构造Lasso回归模型
forestfires.fit(X[train], Y[train])			#训练模型
coef, inte=forestfires.coef_inte()
print('系数：',coef)
print('截距：%.4f'%inte)
score=forestfires.score(X[test], Y[test])	#测试模型
print('对其余%d个样本数据测试，评估分：%.4f'%(m - m // 2, score))

紧接前一程序，本程序第1行的省略号即表示前一程序的代码。第2行创建Lasso回归模型LassoRegressor（程序8.3第2行导入）类对象forestfires。第3行调用forestfires的fit函数用X[train]和Y[train]训练forestfires，注意传递给表示正则化参数$\alpha$的参数alpha的值为0.035。第7行调用forestfires的score函数用X[test]和Y[test]测试forestfires。运行程序，输出

……
LASSO REGRESSION:
训练中...，稍候
39次迭代后完成训练。
系数： [ 6.9598e-03  9.8288e-11  2.5677e-02  8.2371e-03 -9.7036e-09  1.9184e-02-2.4566e-02 -1.6390e-02  1.6269e-02 -2.0274e-02  5.5864e-03 -1.1389e-08-5.0886e-04]
截距：-0.5551
对其余259个样本数据测试，评估分：0.0260

所得模型的评分与线性回归模型的相差无几。现在我们观察Lasso回归模型训练所得的拟合函数系数，其中有若干个的绝对值小于$10^{-4}$，这意味着这些系数对应的特征数据项对标签数据并无贡献。删掉这些数据项，用筛选后的数据集再次训练、测试线性回归模型。

……
a=np.where(np.abs(coef)<1e-4)[0]			#绝对值小于1/10000的系数下标
print('LASSO回归绝对值小于1/10000的系数下标：%s'%a)
X = np.delete(X,a,axis = 1)					#删除对应特征数据列
print('LINE REGRESSION:')
forestfires = LinearRegressor()				#创建线性回归模型
forestfires.fit(X[train], Y[train])			#训练模型
coef, inte = forestfires.coef_inte()
print('系数：%s'%coef)
print('截距：%.4f'%inte)
score=forestfires.score(X[test], Y[test])	#测试模型
print('对其余%d个样本数据测试，评估分：%.4f'%(m - m // 2, score))

紧接前一程序，本程序第1行的省略号表示前两个程序的代码。第2行调用numpy的where函数计算Lasso回归模型的系数数组coef中绝对值小于$10^{-4}$的元素下标，赋予a。第4行调用numpy的delete函数删除特征数据集X中对应a的数据列。第6行创建线性回归模型类对象forestfires。第7行训练forestfires。第11行测试forestfires。运行程序，输出

……
LASSO回归绝对值小于1/10000的系数下标：[ 1  4 11]
LINE REGRESSION:
训练中...，稍候
18次迭代后完成训练。
系数：[ 1.1912  4.4249  1.5407  0.115  -0.06   -1.08    1.0462 -0.24    1.2875]
截距：-14.7009
对其余259个样本数据测试，评估分：0.0302

这意味着特征数据项Y（地图的y坐标）、预警系统的FFMC指数及下雨量rain对应的Lasso回归系数接近零。删除X的这3项数据后，用X[train]和Y[train]经18次迭代训练所得模型对X[test]和Y[test]作测试的评分为0.0302，比原数据（未删除Y、FFMC、rain数据）迭代20次所得模型测试评分0.0248还略有改善。