最优化方法实验四--最小二乘法实验

方不可方可

于 2024-02-18 13:47:01 发布

阅读量1.2k

点赞数 34

分类专栏：最优化方法文章标签：最小二乘法算法线性代数机器学习

本文链接：https://blog.csdn.net/m0_68553640/article/details/136144754

版权

最优化方法专栏收录该内容

5 篇文章 0 订阅

订阅专栏

本文介绍了最小二乘法及其正规方程求解，探讨了梯度下降法在无约束最小二乘问题中的应用，包括QR分解和Gram矩阵的使用，以及如何通过线性搜索确定最优步长。通过实例展示了梯度下降法求解的迭代过程和结果与正规方程的对比。

摘要由CSDN通过智能技术生成

一、实验目的与要求

1.熟练最小二乘法优化模型的意义和求解手段；

2.掌握最小二乘法的正规方程，能实现代码对其求解；

3.掌握最速梯度下降法求解无约束最小二乘法问题。

二、问题

三、模型建立及求解

解决问题思路，模型建立、性能分析，存在问题等方面进行阐述；梯度下降法迭代求解，可以设置迭代次数或相邻迭代解之间“相对接近程度”，如 $\left \| x^{k}-x^{k+1} \right \|_{2}/\left \| x^{k} \right \|_{2}$ ，作为迭代停止条件；代码不要放在报告里面，可以作为附件提交！

1、定理、定义引入

Gram矩阵的定义：

若矩阵 $B=A^{T}A$ ，则B为Gram矩阵。每个Gram矩阵都是半正定的，即

$\forall x:x^{T}Bx=x^{T}A^{T}Ax=\left \| Ax \right \|_{2}^{2}\geq 0$

若要使Gram矩阵为正定的，则要满足：

$\forall x\neq 0:x^{T}Bx=x^{T}A^{T}Ax=\left \| Ax \right \|_{2}^{2}> 0$

即A是列向量无关的。正定矩阵都是非奇异的。

矩阵QR分解的定义（具体参考文章矩阵QR分解）：

QR分解是将一个列向量无关的矩阵 $A\in R^{m\times n}$ 分解成具有标准正交列向量的矩阵Q和上三角矩阵R（对角线元素不为0）的矩阵分解方法，即A=QR：

$A=\begin{bmatrix} a_{1} & a_{2} & \cdots & a_{n} \end{bmatrix} =\begin{bmatrix} q_{1} & q_{2} & \cdots & q_{n} \end{bmatrix}\begin{bmatrix} R_{11}& R_{12} & \cdots & R_{1n}\\ 0& R_{22} & \cdots & R_{2n}\\ \vdots & \vdots & \ddots & \vdots \\ 0& 0 & \cdots & R_{nn} \end{bmatrix}$

$a_{1},a_{2} , \cdots , a_{n}$ 为A的列且线性独立， $q_{1},q_{2} , \cdots , q_{n}$ 为Q的列且两两正交，所以有：

$Q^{T}Q=\begin{bmatrix} q_{1}^{T}q_{1} & q_{1}^{T}q_{2} & \cdots &q_{1}^{T}q_{n} \\ q_{2}^{T}q_{1} & q_{2}^{T}q_{2} & \cdots &q_{2}^{T}q_{n} \\ \vdots & \vdots & \ddots & \vdots \\ q_{n}^{T}q_{1} & q_{n}^{T}q_{2} & \cdots &q_{n}^{T}q_{n} \end{bmatrix}=\begin{bmatrix} 1 & 0 &\cdots & 0\\ 0 & 1 & \cdots & 0\\ \vdots & \vdots &\ddots &\vdots \\ 0 & 0&\cdots & 1 \end{bmatrix}=I$

Q矩阵可逆且 $Q^{T}$ 与Q互为逆矩阵， $R=Q^{-1}A=Q^{T}A$

定理1：设函数f(x)在 $\hat{x}$ 处可微，则 $\triangledown f(\hat{x})=0$ 是 $\hat{x}=\underset{x\in R^{n}}{argmin}f(x)$ 的必要条件

证明：将函数f(x)在 $\hat{x}$ 处一阶泰勒展开，有：

$f(x)=f(\hat{x})+\left \langle \triangledown f(\hat{x}),x-\hat{x} \right \rangle+o(\left \| x-\hat{x} \right \|_{2})$

令 $\tilde{x}=\hat{x}-t\triangledown f(\hat{x})$ ， $t> 0$ ，可得：

$f(\tilde{x})=f(\hat{x})-t\left \| \triangledown f(\hat{x}) \right \|_{2}^{2}+o(t\left \| \triangledown f(\hat{x}) \right \|_{2})$

当 $t\rightarrow 0$ ，则 $t\left \| \triangledown f(\hat{x}) \right \|_{2}^{2}\rightarrow 0$ ， $o(t\left \| \triangledown f(\hat{x}) \right \|_{2})\rightarrow 0$ ，当t足够小，存在 $t\left \| \triangledown f(\hat{x}) \right \|_{2}> o(t\left \| \triangledown f(\hat{x}) \right \|_{2})$ ，即：

$-t\left \| \triangledown f(\hat{x}) \right \|_{2}>+o(t\left \| \triangledown f(\hat{x}) \right \|_{2})<0$

若 $\triangledown f(\hat{x})\neq 0$ ，则： $f(\tilde{x})=f(\hat{x})-t\left \| \triangledown f(\hat{x}) \right \|_{2}^{2}+o(t\left \| \triangledown f(\hat{x}) \right \|_{2})<f(\hat{x})$ ，与 $\hat{x}=\underset{x\in R^{n}}{argmin}f(x)$ 矛盾，所以 $\triangledown f(\hat{x})=0$ 得证。

$\triangledown f(\hat{x})=0$ 是最优解问题 $\hat{x}=\underset{x\in R^{n}}{argmin}f(x)$ 的必要条件但不是充分条件，因为 $\triangledown f(\hat{x})=0$ 时， $f(\hat{x})$ 有可能是最大值，故引入定理2：

定理2：若可微函数f(x)是凸函数，则 $\triangledown f(\hat{x})=0$ 是 $\hat{x}=\underset{x\in R^{n}}{argmin}f(x)$ 的充要条件

定理1已证必要条件，现证充分条件，即 $\triangledown f(\hat{x})=0\Rightarrow \hat{x}=\underset{x\in R^{n}}{argmin}f(x)$ ：由于f(x)是可微且凸的，则 $\forall x\in R^{n}$ ，有： $f(x)\geq f(\hat{x})+\left \langle \triangledown f(\hat{x}),x-\hat{x} \right \rangle\geq f(\hat{x})$ 。 $\hat{x}=\underset{x\in R^{n}}{argmin}f(x)$ 得证。

2、最小二乘法

2.1最小二乘问题

给定 $A\in R^{m\times n}$ ， $b\in R^{m}$ ，若线性方程组 $Ax=b$ 无解，寻找该方程组近似解，并尽可能逼近方程组的目标。设残差向量 $r=Ax-b$ ，求解方程组的近似解，即使残差向量 $r$ 某种度量下最小。当用二范数度量残差大小： $\left \| r \right \|_{2}^{2}=\underset{x}{min}\left \| Ax-b \right \|_{2}^{2}$ ，称其为最小二乘法问题，二乘即指残差二范的平方。

2.2模型建立和求解

设函数 $f(x)=\left \| Ax-b \right \|_{2}^{2}=\sum_{i=1}^{m}(\sum_{j=1}^{n}A_{ij}x_{j}-b_{i})^{2}$ 。要使残差向量 $r$ 的二范数平方最小，则最小二乘问题的最优解 $\hat{x}=\underset{x\in R^{n}}{argmin}f(x)$ ，根据定理2，f(x)是可微且凸的，则有： $\triangledown f(\hat{x})=0$ ，即：

$\triangledown f(\hat{x})=\begin{bmatrix} \frac{\partial f(\hat{x})}{\partial x_{1}}\\ \vdots \\ \frac{\partial f(\hat{x})}{\partial x_{n}} \end{bmatrix}=\begin{bmatrix} 0\\ \vdots \\ 0 \end{bmatrix}$

令 $g_{i}(x)=\sum_{j=1}^{n}A_{ij}x_{j}-b_{i}$ ，则有：

$f(x)=\left \| Ax-b \right \|_{2}^{2}=\sum_{i=1}^{m}(\sum_{j=1}^{n}A_{ij}x_{j}-b_{i})^{2}=\sum_{i=1}^{m}(g_{i}(x))^{2}$

函数f(x)对 $x_{k}$ 的偏导为：

$\frac{\partial f(x)}{\partial x_{k}}=\sum_{i=1}^{m}\left [ (2g_{i}(x))\frac{\partial g_{i}(x)}{\partial x_{k}} \right ]$

而 $\frac{\partial g_{i}(x)}{\partial x_{k}}=A_{ik}$ ，所以：

$\frac{\partial f(x)}{\partial x_{k}}=\sum_{i=1}^{m}2g_{i}(x)A_{ik}=2\sum_{i=1}^{m}\left [ (\sum_{j=1}^{n}A_{ij}x_{j}-b_{i})A_{ik} \right ]$

$=2\sum_{i=1}^{m}\left [ (\sum_{j=1}^{n}A_{ij}x_{j})A_{ik} \right ]-2\sum_{i=1}^{m}b_{i}A_{ik}$

根据矩阵乘法 $\sum_{j=1}^{n}A_{ij}x_{j}$ 即矩阵A第i行乘以向量x，而 $\sum_{i=1}^{m}\left [ (\sum_{j=1}^{n}A_{ij}x_{j})A_{ik} \right ]$ 即为向量Ax的第i个元素乘以 $A_{ik}$ ，再求向量元素之和，如图1所示。

设 $a_{1}$ 为矩阵A的第i列，则有： $\sum_{i=1}^{m}\left [ (\sum_{j=1}^{n}A_{ij}x_{j})A_{ik} \right ]=a_{k}^{T}Ax$

而 $\sum_{i=1}^{m}b_{i}A_{ik}$ 则为向量b第i个元素乘以向量 $a_{k}$ 的第i个元素之和，即向量内积，等于 $a_{k}^{T}b$ 。所以原式等于：

$2\sum_{i=1}^{m}\left [ (\sum_{j=1}^{n}A_{ij}x_{j})A_{ik} \right ]-2\sum_{i=1}^{m}b_{i}A_{ik}=2a_{k}^{T}Ax-2a_{k}^{T}b$ $=2a_{k}^{T}(Ax-b)$

则 $\triangledown f(x)$ 等于：

$\triangledown f(x)=\begin{bmatrix} \frac{\partial f(x)}{\partial x_{1}}\\ \vdots \\ \frac{\partial f(x)}{\partial x_{n}} \end{bmatrix}=2\begin{bmatrix} a_{1}^{T}(Ax-b)\\ \vdots \\ a_{n}^{T}(Ax-b) \end{bmatrix}=2A^{T}(Ax-b)$

令 $f(x)=2A^{T}(Ax-b)=0$ ，则 $A^{T}Ax=A^{T}b$ ，称该等式为正规方程。

根据Gram矩阵的定义，当A列向量无关时， $A^{T}A$ 是正定的，即非奇异的，所以：

$\hat{x}=(A^{T}A)^{-1}A^{T}b$

该模型的几何解释为，如图2所示，range(A)为A的向量张成空间， $A\hat{x}$ 为该空间的一个向量。b为空间中任一向量。若要使 $A\hat{x}$ 和b两个向量的差的二范数平方最小，可知 $A\hat{x}$ 即为b在range(A)上的投影： $A(A^{T}A)^{-1}A^{T}b$ 。

2.3根据求解的模型解决最小二乘问题

由上述分析可得，最小二乘问题 $\underset{x}{argminf(x)}$ 的最优解为： $\hat{x}=(A^{T}A)^{-1}A^{T}b$ ，根据该等式，可以根据构造Gram矩阵或QR分解A矩阵来求解 $\hat{x}$ 。

①构造Gram矩阵方法

根据正规方程 $A^{T}Ax=A^{T}b$ ，首先构造Gram矩阵 $A^{T}A$ ，然后计算 $A^{T}b$ ，再通过高斯消元求解线性方程组的方法求解 $\hat{x}$ 。

计算Gram矩阵 $A^{T}A$ 的时间复杂度约为 $2mn^{2}$ flops，计算 $A^{T}b$ 的复杂度为2mnflops，高斯消元求解的复杂度为约为 $n^{3}$ 。总的时间复杂度约为 $(2mn^{2}+n^{3})$ flops。

读取附件“MatrixA_b.mat”文件中的矩阵A和向量b，按上述方法求解，循环20次，计算每次的残差向量大小： $\left \| r \right \|_{2}^{2}=\left \| Ax-b \right \|_{2}^{2}$ ，画出折线图如图3所示，可以看到平均残差约1.1738。

②QR分解方法

当A列向量无关时，首先对A进行QR分解，根据QR分解的定义， $A=QR,Q^{T}Q=I$ ，所以根据等式 $\hat{x}=(A^{T}A)^{-1}A^{T}b$ 有：

$\hat{x}=(A^{T}A)^{-1}A^{T}b=((QR)^{T}(QR))^{-1}(QR)^{T}b$

$=(R^{T}Q^{T}QR)^{-1}R^{T}Q^{T}b=(R^{T}R)^{-1}R^{T}Q^{T}b=R^{-1}Q^{T}b$

得到线性方程组： $Rx=Q^{T}b$ 。然后计算 $Q^{T}b$ ，由于R是一个上三角矩阵，所以可以直接回代求解线性方程组。

矩阵A的QR分解时间复杂度约为 $2mn^{2}$ flops，计算 $Q^{T}b$ 的复杂度约为2mnflops，回代求解 $Rx=Q^{T}b$ 的复杂度约为 $n^{2}$ flops。总的时间复杂度约为 $2mn^{2}$ flops。

由图3和图4可以看到，两种方法求解的 $\hat{x}$ ，计算出的残差数值相等，但是通过Gram矩阵求解的方法构造的Gram矩阵在经过四舍五入后，由于精度问题有可能是奇异矩阵，所以QR分解求解更稳定，且时间复杂度更低。

3、梯度下降法

3.1沿梯度方向下降

当 $A\in R^{m\times n}$ 列向量相关或n非常大时，通过构造Gram矩阵或QR分解的方法无法求解或需要消耗较多时间。故引入梯度下降法求解。

梯度下降法的思想为：通过迭代求解目标的最优解过程： $x^{1},x^{2},...,x^{k},....\rightarrow \hat{x}$ ，如图5所示。 $x^{k}$ 是第k步迭代，期望更新 $x^{k+1}$ ，满足 $f(x^{k+1})<f(x^{k})$ 。

迭代退出条件，可以设置最大迭代次数，或设置相邻迭代之间的相对接近程度，如当 $\left \| x^{k}-x^{k+1} \right \|_{2}/\left \| x^{k} \right \|_{2}$ 小于设置的阈值时，可以认为目标函数值收敛到一个稳定值，算法已经达到最优解。

在图3所示函数图像中，导数代表切线斜率，还可以代表梯度方向。函数从局部最小值处沿着梯度逐渐增大。所以要最小化函数f(x)，可以从 $x_{1}$ 处沿着梯度的负方向减小函数值，且沿着梯度方向函数值下降得最快，则需要求出 $x_{1}$ 处的导数取其负值，然后再乘以一个系数 $\alpha _{1}$ ， $\alpha _{1}>0$ ，，通常叫做步长或学习率。 $\alpha _{k}$ 为第k次迭代的步长。令 $x^{2}=x^{1}-\alpha _{1}\triangledown f(x)$ 。

可以证明当 $x^{2}=x^{1}-\alpha _{1}\triangledown f(x)$ 时满足 $f(x^{2})<f(x^{1})$ ：函数f(x)可微，根据泰勒公式，在 $x^{1}$ 的一阶公式为：

$f(x^{2})=f(x^{1})+\left \langle \triangledown f(x^{1}),x^{2}-x^{1} \right \rangle+o(\left \| x^{2}-x^{1} \right \|)$

当 $\left \| x^{2}-x^{1} \right \|$ 足够小，则有：

$f(x^{2})-f(x^{1})\approx \left \langle \triangledown f(x^{1}),x^{2}-x^{1} \right \rangle=\left \langle \triangledown f(x^{1}),-\alpha _{1}\triangledown f(x) \right \rangle$ $=-\alpha _{1}\left \| \triangledown f(x) \right \|_{2}^{2}<0$

同理，在第k步迭代时，令 $x^{k+1}=x^{k}-\alpha _{k}\triangledown f(x)$ ，可证 $f(x^{k+1})<f(x^{k})$ 。f(x)在 $x^{k}$ 的一阶泰勒展开公式为：

$f(x^{k+1})=f(x^{k})+\left \langle \triangledown f(x^{k}),x^{k+1}-x^{k} \right \rangle+o(\left \| x^{k+1}-x^{k} \right \|)$

$f(x^{k+1})-f(x^{k})\approx \left \langle \triangledown f(x^{k}),x^{k+1}-x^{k} \right \rangle=\left \langle \triangledown f(x^{k}),-\alpha _{k}\triangledown f(x) \right \rangle$

$=-\alpha _{k}\left \| \triangledown f(x) \right \|_{2}^{2}<0$

3.2最优步长

步长 $\alpha _{k}$ 的取值影响求最优解的速度。如果步长过小，可能导致算法收敛缓慢，还可能使得算法无法跳出局部最优点，从而无法找到全局最优解；如果步长过大，可能导致算法无法收敛，还可能使得算法在搜索空间中跳过最优解，无法稳定地接近最优解。

可以通过线性搜索估计步长 $\alpha _{k}$ ： $\alpha _{k}=\underset{\alpha \in R}{argmin}f(x^{k}-\alpha _{k}\triangledown f(x))$ ，即 $\alpha _{k}$ 为最优步长。令 $g(\alpha _{k})=f(x^{k}-\alpha _{k}\triangledown f(x))$ ， $g(\alpha )$ 是关于 $\alpha$ 的可微凸函数，根据定理2有：当 $g^{'}(\hat{\alpha _{k}})=0$ 时， $\hat{\alpha _{k}}=\underset{\alpha \in R}{argmin}f(x^{k}-\alpha _{k}\triangledown f(x))$ 。由于 $\underset{x}{min}\left \| Ax-b \right \|_{2}^{2}$ 与 $\underset{x}{min}\frac{1}{2}\left \| Ax-b \right \|_{2}^{2}$ 有相同解，不妨设：

$f(x)=\frac{1}{2}\left \| Ax-b \right \|_{2}^{2}, \triangledown f(x)=A^{T}(Ax-b)$

则 $g(\alpha _{k})$ 和 $g^{'}({\alpha _{k}})$ 为：

$g(\alpha _{k})=f(x^{k}-\alpha _{k}\triangledown f(x))=\left \| A(x^{k}-\alpha _{k}A^{T}(Ax^{k}-b))-b \right \|_{2}^{2}$

$=\left \| Ax^{k}-\alpha _{k}AA^{T}(Ax^{k}-b)-b \right \|_{2}^{2}=\left \| (I-\alpha _{k}AA^{T})(Ax^{k}-b) \right \|_{2}^{2}$

$=(Ax^{k}-b)^{T}(I-\alpha _{k}AA^{T})^{T}(I-\alpha _{k}AA)(Ax^{k}-b)$

$=(Ax^{k}-b)^{T}(I-2\alpha _{k}AA^{T}+\alpha _{k}^{2}AA^{T}AA^{T})(Ax^{k}-b)$

$=(Ax^{k}-b)^{T}I(Ax^{k}-b)-(Ax^{k}-b)^{T}2\alpha _{k}AA^{T}(Ax^{k}-b)$

$+(Ax^{k}-b)^{T}\alpha _{k}^{2}AA^{T}AA^{T}(Ax^{k}-b)$

$=\left \| Ax^{k}-b \right \|_{2}^{2}-2\alpha _{k}\left \|A^{T} (Ax^{k}-b) \right \|_{2}^{2}+\alpha _{k}^{2}\left \|AA^{T} (Ax^{k}-b) \right \|_{2}^{2}$

$g^{'}(\alpha _{k})=-2\left \|A^{T} (Ax^{k}-b) \right \|_{2}^{2}+2\alpha _{k}\left \|AA^{T} (Ax^{k}-b) \right \|_{2}^{2}$

令 $g^{'}(\alpha _{k})=0$ ，求得： $\alpha _{k}=\frac{\left \|A^{T} (Ax^{k}-b) \right \|_{2}^{2}}{\left \|AA^{T} (Ax^{k}-b) \right \|_{2}^{2}}=\frac{\left \|\triangledown f(x^{k}) \right \|_{2}^{2}}{\left \|A\triangledown f(x^{k}) \right \|_{2}^{2}}$ 。

3.3算法步骤

总结用梯度下降法求解的步骤如下：

3.4测试结果

读取附件“MatrixA_b.mat”文件中的矩阵A和向量b，使用梯度下降法求解 $\hat{x}$ ，设置最大迭代次数为500，计算每一次迭代中的残差向量的大小： $\left \| r \right \|_{2}^{2}=\left \| Ax-b \right \|_{2}^{2}$ ，画出折线图如图6所示。

可以看到随着迭代残差向量逐渐变小，当迭代次数达到100次的时候残差向量的大小趋于平缓，退出迭代时其大小与使用正规方程求解的残差1.1738稍微有一点差异，增加迭代次数至800次，残差大小与1.1738几乎相等。

计算梯度下降法求解出的近似解与使用正规方程求解得到的 $\hat{x}$ 每个元素之差的绝对值，结果如图7所示，可以看到差异较小，证明使用梯度下降法计算的近似解是可行的。

更改迭代退出条件为 $\left \| x^{k}-x^{k+1} \right \|_{2}/\left \| x^{k} \right \|_{2}\leq p$ ，p为阈值。当 $p=0.00001$ 时求得的近似解与准确解十分接近。计算每一次迭代中的残差向量的大小，画出折线图如图8所示，可以看到共迭代了761次。计算近似解与准确解的 $\hat{x}$ 每个元素之差的绝对值，结果如图9所示，可以看到差异较小。