最优化方法Python计算：无约束优化应用——线性回归模型

一、监督学习模型

回归算法是典型的监督学习模型之一。 需要指出的是，在机器学习模型中，无论是训练部分还是预测部分，都需要对数据作一些规范化处理。首先，将样本特征向量${\mathbf{x}}_i=\left(\begin{array}{c}{x}_{i1}\\ x_{i2}\\ ⋮\\ x_{in}\end{array}\right)\in {\text{R}}^n,i=1,2,\cdots ,m$将其组织成一个$m\times n$矩阵
$$\mathbf{X}=\left(\begin{array}{c}{\mathbf{x}}_1^{\top }\\ {\mathbf{x}}_2^{\top }\\ ⋮\\ {\mathbf{x}}_m^{\top }\end{array}\right)=\left(\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1n}\\ x_{21}& x_{22}& \cdots & x_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ x_{m1}& x_{m2}& \cdots & x_{mn}\end{array}\right)$$
使得表达式展开时更简洁。其次，要对数据进行“标准化”处理。以消除不同量纲单位带来的数据偏差，并使数据指标处于同一数量级，更适合综合对比评价。此处以归一化作为规范化处理方式：
(1)训练阶段。此时，样本特征数据${\mathbf{x}}_i,i=1,2,\cdots ,m$且$m>1$。对每个$1\le j\le n$，第$j$列数据$\left(\begin{array}{c}{x}_{1j}\\ x_{2j}\\ ⋮\\ x_{mj}\end{array}\right)$表示$m$个样本第$j$个特征数据。计算每一列的最小值与最大值 min ⁡ x j = min ⁡ i { x i j } \min x_j=\min\limits_{i}\{x_{ij}\} minxj​=imin​{xij​}及 max ⁡ x j = max ⁡ i { x i j } \max x_j=\max\limits_{i}\{x_{ij}\} maxxj​=imax​{xij​}。并记
$$\{\begin{array}{l}\min \mathbf{x}=(\min x_1,\min x_2,\cdots ,\min x_n)\\ \max \mathbf{x}=(\max x_1,\max x_2,\cdots ,\max x_n)\end{array}$$
令$\Delta \mathbf{x}=\max \mathbf{x}-\min \mathbf{x}=(\max x_1-\min x_1,\max x_2-\min x_2,\cdots ,\max x_n-\min x_n)$，以
$$\frac{{\mathbf{x}}_i^{\top }-\min \mathbf{x}}{\Delta \mathbf{x}}=\left(\frac{x_1-\min x_1}{\max x_1-\min x_1},\frac{x_2-\min x_2}{\max x_2-\min x_2},\cdots ,\frac{x_n-\min x_n}{\max x_n-\min x_n}\right)$$
作为${\mathbf{x}}_i^{\top }$归一化后的向量，$i=1,2,\cdots ,m$。
相仿地，对接收到的标签数据$y_i$，$i=1,2,\cdots ,m$，记 min ⁡ y = min ⁡ i { y i } , max ⁡ y = max ⁡ i { y i } \min y=\min\limits_{i}\{y_i\},\max y=\max\limits_{i}\{y_i\} miny=imin​{yi​},maxy=imax​{yi​}，以
$$\frac{y_i-\min y}{\max y-\min y}$$
作为$y_i$归一化后的值。数据归一化后，其值均介于0,1之间。
(2)预测阶段。此时，只有新样本的特征向量${\mathbf{x}}_i$，$i=1,2,\cdots ,m$。进行归一化时，需要用训练时算得的${\min }_j$和${\max }_j$进行计算。这是因为，我们认为训练数据是总体的简单样本，其统计特征值（最小值、最大值）表示总体的近似分布。而检测数据也来自同一总体，出于一致性考虑，使用训练时取得的最大值、最小值对检测数据作归一化计算。下列代码将监督学习模型表示为一个Python类

import numpy as np										#导入numpy
from abc import ABC, abstractmethod						#导入abc模块
from scipy.optimize import minimize						#导入minimize
class SupervisedLearningModel(ABC):						#监督学习模型抽象类
    def xnormalize(self, x, trained):					#样本特征数据归一化方法
        if not trained:									#训练前
            xmin = np.min(x,axis = 0)					#按列计算最小值
            xmax = np.max(x,axis = 0)					#按列计算最大值
            self.xmin = xmin							#记录计算结果
            self.xmax = xmax
        else:											#测试或预测前
            xmin = self.xmin							#使用训练时记录的数据
            xmax = self.xmax
        return (x - xmin) / (xmax - xmin) 
    def ynormalize(self, y, trained):					#样本标签归一化方法
        if not trained:									#训练前
            self.ymin = np.min(y)						#记录最小值
            self.ymax = np.max(y)						#记录最大值
        return (y - self.ymin)/(self.ymax - self.ymin)
    def pretreat(self, x, y = None, trained = False):	#数据预处理函数
        if isinstance(y, np.ndarray):					#需处理样本标签
            y = self.ynormalize(y, trained)				#归一化标签
        if not isinstance(x, np.ndarray):				#一个一元样本
            x = np.array([x]).reshape(1, 1)
        else:
            if len(x.shape) == 1:						#一维数组
                if self.scalar:							#多个一元样本
                    x = x.reshape(x.size,1)
                else:									#一个多元样本
                    x = x.reshape(1,x.size)
        x = self.xnormalize(x, trained)					#归一化样本特征
        return x, y
    def predict(self, X):								#预测方法
        X, _ = self.pretreat(X, trained = True)			#归一化样本特征数据
        yp = self.F(X)									#计算拟合函数值
        if yp.size == 1:								#单样本预测值
            yp = yp[0]
        return self.tagVal(yp * (self.ymax - self.ymin) + self.ymin)
    @abstractmethod
    def F(self, w = None, x = None):					#决策函数
        pass
    @abstractmethod
    def obj(self, w):									#目标函数
        pass
    @abstractmethod
    def fit(self, X, Y, w = None):						#训练函数
        pass
    @abstractmethod
    def score(self, x, y):								#模型测试函数
        pass

程序的第4~50行定义的监督学习模型类SupervisedLearningModel包含了七个函数：xnormalize、ynormalize、pretreat、predict、F、obj和fit均为对象函数（含有参数self）。其中，决策函数F、目标函数obj和训练函数fit由指令@abstractmethod（第2行导入）指引，为抽象函数。所以，SupervisedLearningModel是一个抽象类ABC（第2行导入）。程序中

• 第5~14行的xnormalize函数用于对样本特征数据进行归一化处理。输入参数x是一个组织成$m\times n$矩阵表示样本数据。布尔型参数trained指示归一化操作是在训练前还是训练后进行。若trained=False，则表示训练前的归一化操作；若trained=True，则表示测试或预测前的归一化操作。函数体内第7~10行针对训练前的归一化操作：第7、8行对x分别调用Numpy的min及max函数，计算各列（由参数axis=0决定）的最小值和最大值，存于数组xmin和xmax中。第9、10行将算得的xmin和xmax记录为对象的同名属性（self.xmin、self.xmax）。第12、13行则是针对测试或预测前样本特征数据的归一化操作：简记训练时记录下的最小值、最大值为xmin、xmax。第14行计算x的归一化，并返回。
• 第15~19行的ynormalize函数用于归一化标签数据y。若这一操作仅在训练前进行，第17、18行直接算得其最小值、最大值记录为对象属性self.ymin和self.ymax。第19行计算标签数据的归一化，并返回。
• 第20~32行的pretreat函数负责对数据进行预处理，即归一化操作。参数x表示样本特征数据；y表示样本标签数据，（缺省值为None，即预测阶段无需传递标签数据）；布尔型参数trained指示预处理操作是在训练前还是训练后进行，缺省值为False，表示训练前的预处理。第21~22行在训练或测试时调用ynormalize函数归一化标签数据y。第23~30行将表示样本特征数据的参数x归一化为$m\times n$矩阵。具体而言，第23行检测到样本仅含一个特征，且x含一个样本，即x为一数值时，第24行将其转换成1行1列矩阵[[x]]。第26行测得x为一维数组，这又分两种情形：第27行检测表示是否为一元（样本仅含一个特征）模型的属性scalar。若是，第28行则将x转换为$m\times 1$矩阵。否则，即x表示一个$n$元样本，第30行将x转换为$1\times n$矩阵。然后，第31行调用xnormalize函数，对x作归一化操作。最终，第32行将归一化后的x和y作为返回值返回。
• 第33~38行的predict函数是监督学习模型的一个核心操作：训练后，用所得模型参数${\mathbf{w}}_0$代入回归函数，得到拟合函数$F({\mathbf{w}}_0;\mathbf{x})$。对新的样本特征$\mathbf{x}$，预测对应值$y=F({\mathbf{w}}_0;\mathbf{x})$。函数只有一个外部参数X表示新的样本特征数据。第34行调用pretreat函数，对X作规范化处理。同时注意传递给参数trained的值为True，表示此处的数据预处理是在训练后进行的。第35行将X传递给拟合函数F（利用训练时算得的模型参数w0），计算预测值赋予yp。由于计算得到的是一个数组，如果算得结果仅含一个数值，第37行将其转换成数值。为了与原问题保持一致性，第38行将预测值yp通过逆归一化还原到原始数据的量级，并作为标签值函数tagVal的参数返回。标签值函数tagVal的意义在后面的内容展开中解释。
• 第39~41行的抽象函数F是监督学习模型的决策函数。该函数有两个参数：表示样本特征数据的x和表示模型参数向量$\mathbf{w}$的初始点w，它们的缺省值均为None。该函数的具体实现由子类实现。
• \item第42~44行的抽象函数obj是监督学习模型中计算最优模型参数的目标函数。该函数仅含一个参数w，表示模型参数向量。该函数的具体实现由子类实现。
• 第45~47行的抽象函数fit是监督学习模型的训练函数。该函数有三个外部参数：表示样本数据的X，标签数据的Y，和模型参数向量$\mathbf{w}$的初始点w，其缺省值为None。该函数的具体实现由子类实现。
• 第48~50行的抽象函数score是监督学习模型的测试函数。该函数有两个参数：表示样本特征数据的x和表示样本标签数据的y。该函数的具体实现由子类实现。

二、最小二乘法模型

回归是一种统计学方法，用于根据样本数据$({\mathbf{x}}_i,y_i)$，$i=1,2,\cdots ,m$，探究变量$\mathbf{x}$与$y$之间的关系。具体而言，回归模型的任务是找出拟合函数$F(\mathbf{x})$，使得
$$y_i\approx F({\mathbf{x}}_i),i=1,2,\cdots ,m$$
并用$F(\mathbf{x})$来对新的输入$\mathbf{x}\in {\text{R}}^n$预测对应的输出$y\in \text{R}$。寻求拟合函数的过程，通常是选择一个具有待定参数$\mathbf{w}$的函数$F(\mathbf{w};\mathbf{x})$，其中$\mathbf{w}\in {\text{R}}^p,p\in \text{N}$。然后计算使得$y_i\approx F({\mathbf{w}}_0;{\mathbf{x}}_i),i=1,2,\cdots ,m$最“合适”的参数${\mathbf{w}}_0$——称为该学习模型的模式，并将$F({\mathbf{w}}_0;\mathbf{x})$作为拟合函数。
此处所谓最“合适”的模式，可以从不同的角度去认知。如果从欧氏空间的集合观点，可得出回归过程最常用的最小二乘法。给定序列$({\mathbf{x}}_i,y_i)$，$i=1,2,\cdots ,m$，最小二乘法对所选含有待定参数$\mathbf{w}\in {\text{R}}^p$的函数$F(\mathbf{w};\mathbf{x})$，记$\mathbf{F}(\mathbf{w})=\left(\begin{array}{c}F(\mathbf{w};{\mathbf{x}}_1)\\ F(\mathbf{w};{\mathbf{x}}_2)\\ ⋮\\ F(\mathbf{w};{\mathbf{x}}_m)\end{array}\right)$，$\mathbf{y}=\left(\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right)$，解无约束最优化问题
$$\{\begin{array}{l}\min \Vert \mathbf{F}(\mathbf{w})-\mathbf{y}{\Vert }^2\\ \text{s.t}\mathbf{w}\in {\text{R}}^p\end{array},$$
设${\mathbf{w}}_0=\arg \underset{\mathbf{w}\in {\text{R}}^p}{\min }\Vert \mathbf{F}(\mathbf{w})-\mathbf{y}{\Vert }^2$，则$F({\mathbf{w}}_0,\mathbf{x})$即为所求的拟合函数。
下列代码将最小二乘法模型实现为SupervisedLearningModel的子类：

import numpy as np
from scipy.optimize import minimize
class LeastSquar(SupervisedLearningModel):				#最小二乘法模型
    def obj(self,w):									#目标函数
        return np.linalg.norm(self.F(w) - (self.Y)) ** 2
    def fit(self, X, Y, w = None):						#训练方法
        print("训练中...，稍候")
        self.scalar = (len(X.shape) == 1)				#记录一元样本特征
        self.X, self.Y = self.pretreat(X, Y)			#预处理训练数据
        p = self.w0len()								#计算模型参数长度
        if not isinstance(w1, np.ndarray):
            if w == None:								#未传递初始参数向量w
                w = np.random.random(p)					#随机产生初始向量w
            else:										#w是常数
                w = np.array([w] * p)
        res = minimize(self.obj, w)						#解最优化问题
        self.w0 = res.x									#记录模型参数
        print("%d次迭代后完成训练。"%res.nit)

第3~18行定义的LeastSquar类作为SupervisedLearningModel的子类。其中

• 第4~5行的obj函数实现了最小二乘法的优化目标函数。该函数有一个参数w，表示模型参数向量$\mathbf{w}$。第5行调用numpy.linalg模块中计算向量范数的norm函数，计算$\Vert F(\mathbf{w};\mathbf{A})-\mathbf{y}{\Vert }^2$。其中，$\mathbf{A}$为训练用的样本特征增广矩阵$\left(\begin{array}{c}{\mathbf{x}}_1^{\top },1\\ {\mathbf{x}}_2^{\top },1\\ ⋮\\ {\mathbf{x}}_m^{\top },1\end{array}\right)$，$\mathbf{y}$为训练用的标签数据向量$\left(\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right)$。
• 第6~18行的fit函数是回归模型中的核心操作之一：用最优化方法确定回归函数$F(\mathbf{w};\mathbf{x})$中的待定参数${\mathbf{w}}_0$。该函数有三个外部参数：表示样本数据的X，标签数据的Y，和模型参数向量$\mathbf{w}$的初始点w，其缺省值为None。第8行根据X的形状确定模型是否是一元的（样本仅含一个特征），记录于scalar。第9行调用pretreat函数规范化特征数据X和标签数据Y，记录于对象的X和Y属性中。注意此处对trained参数，使用起缺省值False，表示数据预处理是在训练前进行的。第10行调用w0len函数计算模型参数长度p。第12~13行对未确定模型参数向量初始点（w == None）情形，第13行调用Numpy.random的random函数生成具有p个介于0，1之间的随机数的数组赋予w。第14~15针对传递进来的w为一常数的情形：将w生成为含有p个指定常数值的数组。第16行调用minimize函数（第2行导入），用默认的系统方法计算目标函数obj，从初始点w开始的最优解。第17行记录算得的最优模型参数为对象属性w0。

三、线性回归模型

回归模型中，若拟合函数选择为线性函数，即
$$y=F(\mathbf{w};\mathbf{x})=\sum _{i=1}^n{w}_i{x}_i+w_{n+1}=({\mathbf{x}}^{\top },1)\mathbf{w}$$
则称为线性回归模型。其中，$w_1,w_2,\cdots ,w_n$为$x_1,x_2,\cdots ,x_n$的加权和$\sum _{i=1}^n{w}_i{x}_i$的系数，$w_{n+1}$为偏移量。线性回归模型可图示化为

下面，将线性回归模型定义成LeastSquar的子类LineModel。

import numpy as np										#导入numpy
class LineModel(LeastSquar):							#线性回归模型
    def pretreat(self, x, y = None, trained = False):	#重载预处理函数
        X, Y = super().pretreat(x, y, trained)			#执行父类预处理函数
        X = np.hstack((X, np.ones((X.shape[0], 1))))	#构造特征数据增广矩阵
        return X, Y
    def w0len(self):									#模型参数维度函数
        return self.X.shape[1]
    def F(self, w = None, x = None):					#线性拟合函数
        if w is None:
            w = self.w0
        if x is None:
            x = self.X
        return x @ w
    def coef_inte(self):								#系数与截距
        xmin, xmax = self.xmin, self.xmax				#特征数据最大/最小值
        ymin, ymax = self.ymin, self.ymax				#标签数据最大最小值
        dx = xmax - xmin								#x增量
        dy = ymax - ymin								#y增量
        w = self.w0										#模型参数
        n = self.w0len()								#模型参数长度
        coef = dy / dx * w[0 : n - 1]					#系数
        inte = dy * w[n - 1] + ymin						#截距
        - dy * np.dot(xmin / dx, w[0 : n - 1])
        return coef, inte
    def fit(self, X, Y, w = None):						#重载训练函数
        super().fit(X, Y, w)							#执行父类训练函数
        self.coef_, self.intercept_ = self.coef_inte()	#系数与截距

程序的第2~28行定义的线性模型LineModel类是程序4.2定义的最小二乘法模型类LeastSquar的子类。其中，

• 第3~6行重载了LeastSquar类的预处理函数pretreat。该函数的参数x表示特征数据，y表示标签数据，trained表示是否已经训练。第5行调用父类的预处理函数pretreat，计算增广矩阵X和标签数据Y。第5行将X的最后一列添加1，构造增广矩阵$\mathbf{X}=\left(\begin{array}{c}{\mathbf{x}}_1^{\top },1\\ {\mathbf{x}}_2^{\top },1\\ ⋮\\ {\mathbf{x}}_m^{\top },1\end{array}\right)$。第6行返回增广矩阵X和标签数据Y。
  -第7~8行的w0len函数定义中，返回增广矩阵$\mathbf{X}$的列数，即模型参数长度。
• 第9~14行的拟合函数F定义中，参数x表示增广矩阵$\mathbf{X}$，w表示拟合函数的待定参数$\mathbf{w}$。第14行计算拟合函数$F(\mathbf{w};\mathbf{x})=({\mathbf{x}}^{\top },1)\mathbf{w}$并返回。
• 第15~25行的coef_inte函数定义中，第16、17行分别读取对象的关于特征数据最小值/最大值xmin/xmax及标签数据最小值/最大值ymin/ymax。第18、19行计算$\Delta x$和$\Delta y$，分别赋予dx，dy。第20行读取对象的模型参数赋予w。第21行计算模型参数长度赋予n。第22、23~24行计算线性回归模型的系数coef与截距inte。
• 第26~28行重载训练函数fit。其中，第27行执行父类的训练函数fit，第28行调用coef_inte函数计算拟合函数的系数和截距分别赋予本对象属性coef_和intercept_。
  写博不易，敬请支持：
  如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！