基于广义逆矩阵理论求解线性最小二乘问题

在介绍广义逆矩阵之前，我们首先讨论如下所示的一个线性方程组的求解问题：

其中A是m×n的矩阵。根据线性代数的知识，我们知道方程组有解的充要条件为：

此时我们称该方程组为相容方程组。如果方程组Ax=b中的A是一个可逆矩阵，那么方程组的解可表示为：

并且解是唯一的。这样的可逆首先要求矩阵A是方阵，并且A的行列式不为0。对于矛盾方程组，矩阵A可能并不为方阵或者其行列式等于0，显然在这种情况下矩阵A并不存在逆矩阵。因此我们希望推广可逆的概念，并且从某个角度看它与传统定义本质上式吻合的。也就是说，我们希望找到矩阵A的广义逆矩阵G，使得线性方程组的解仍然能够表示为：

在曲线拟合和多元线性回归分析等问题中，常常要计算矛盾方程组Ax=b的最小二乘解。广义逆矩阵理论可以使得求矛盾方程组Ax=b的最小二乘解的方法标准化、规范化，整个求解过程归结为求A的广义逆，而不用求误差平方和的极值等一套繁琐的过程。本文将主要给出广义逆矩阵的一些基本概念及定理，并重点介绍广义逆矩阵的计算方法及其在线性最小二乘问题中的应用。

的全部或者一部分，则称G是A的广义逆矩阵，也称为伪逆。

上式中条件(3)(4)中的上标H表示转置共轭，它是针对复矩阵（矩阵中含有复数）而言的；对于实矩阵，其与转置（上标T表示）是等价的。

按照这个定义，广义逆矩阵可以是满足上述条件的一部分，或者是满足全部条件的矩阵，因此共有15类广义矩阵。应用最多的是满足全部条件的情况，记为A⁺，称为加号逆或极小最小二乘广义逆；满足条件(1)和(4)的情况，记为A_l^-，称为最小二乘广义逆；满足条件(1)的情况，记为A^-，称为减号逆。

广义逆矩阵A^-起源于线性方程组Ax=b（其中A是m×n矩阵）的求解问题。

定理1

对任意的b∈R(A)，存在矩阵G使得x=Gb都是方程Ax=b的解的充分必要条件是G满足

称G是A的减号逆，记为A^-。

在实际计算中，可以采用如下定理计算矩阵A的减号逆。

定理2

设A是m×n矩阵，rank(A)=r，非奇异矩阵P∈C^m×m，Q∈C^n×n使得

则G∈C^n×m使得AGA=A的充分必要条件是

其中G₁₂，G₂₁，G₂₂分别是r×(m-r),(n-r)×r,(n-r)×(m-r)的任意矩阵，E_r为秩为r的单位矩阵。

由上述定理可知，对于任意的A∈C^m×n，它的减号逆总存在，并且并不唯一。

计算示例

下面给出一个简单的示例，试求A的减号逆A^-，其中A为：

首先通过对矩阵A做初等行变换将其化为阶梯形矩阵：

由于变换后矩阵A的非零行为2，因此矩阵A的秩为rank(A) =2。

因此我们需要寻找非奇异矩阵P和Q使得

为此我们构造如下矩阵并首先做初等行变换：

为了将左上角的原矩阵A继续化为最简形，我们继续做初等列变换：

因此有

因此矩阵A的减号逆A^-可表示为：

式中：g₁₃, g₂₃, g₃₁, g₃₂, g₃₃为任意常数。

在计算得到矩阵A的任意一个减号逆A^-后，则相容方程组Ax=b的通解可表示为：

从上一节可以看到，相容方程组Ax=b有无穷多个解。在所有的解中2-范数最小的解x₀，即满足

的解x₀，我们称为相容方程组Ax=b的极小范数解。

因此我们希望找到矩阵G，使得对于任意的b∈R(A)，x=Gb都是相容方程组Ax=b的极小范数解，这样的矩阵G我们将其称为矩阵A的极小范数广义逆，记为A_m^-。可以证明，这样的矩阵G必须满足：

可以通过如下所示的定理来计算矩阵A的极小范数广义逆：

设Ax=b是矛盾（不相容）方程，它在通常意义下无解。现要求这样的解x₀，使得：

这样的解x₀，称为矛盾方程Ax=b的最小二乘解。

现在，我们考虑的问题是：是否存在矩阵G，使得对任意的b∈Cⁿ，满足x=Gb都是方程Ax=b的最小二乘解？

定理1

设A∈C^m×n，G∈C^n×m。对任意的b∈Cⁿ，x=Gb都是方程Ax=b的最小二乘解的充分必要条件是G满足：

此时称G是A的最小二乘广义逆，记作A_l^-。

可以证明G满足

定理2

矛盾方程组Ax=b的通解是

可以看到，该通解一般情况下并不唯一。但如果矩阵A满足列满秩，则Ax=b的最小二乘解唯一存在，且满足：

在曲线拟合和多元线性回归分析中，常常要计算Ax=b的最小二乘解。基于广义逆矩阵理论，可以将该问题转化为求解矩阵A的最小二乘广义逆矩阵，从而使得求解过程变得规范化。

由前面给出的定理可知，方程Ax=b的最小二乘解一般并不唯一，其所有最小二乘解中范数最小的解，称为方程Ax=b的极小最小二乘解。

因此对于任意的b∈C^m，我们希望找到G∈C^n×m，使得x=Gb是方程Ax=b的极小最小二乘解。

定理1

对于任意的b∈C^m，方程Ax=b的最小二乘解是：

可以证明对于A∈C^m×n，G=A_m^-AA_l^-满足全部的Penrose-Moore方程，我们将其记为A⁺，称为加号逆或Penrose-Moore逆矩阵。

可以证明加号逆存在如下性质：

我们可以基于如下定理通过数值方法来求解加号逆。

定理2 奇异值分解法

设矩阵A∈C^m×n，存在m阶正交阵U和n阶正交阵V使得

我们将其称为矩阵A的奇异值分解。其中

其中λ_i是按照从大到小排列的矩阵A^TA的非零特征值的全体。

求得了矩阵的奇异值分解，则可以知道矩阵的加号逆为：

定理3 Lagrange-Sylvester展开式

设a₁, a₂, …, a_s是A^TA的相异特征值，令：

有

则有

其中

当然，也可以直接采用MATLAB中的pinv命令来求解加号逆。

计算示例

下面给出一个简单的计算示例，介绍如何通过奇异值分解法来求解矩阵的加号逆。考虑矩阵

首先我们需要对矩阵进行奇异值分解。奇异值分解可以看作是对矩阵特征值分解的推广。矩阵特征分解的形式为：

其中Q为特征向量矩阵，Σ为特征值对角阵。

特征值分解的本质是求解给定矩阵的特征值和特征向量，提取出矩阵最重要的特征。但特征值分解仅适用于方阵，很多情况下矩阵A并不一定为方阵，因此我们需要矩阵A右乘其转置矩阵，即A^TA，来获得一个方阵，这样就可以进行特征值分解。可以严格证明，矩阵A^TA的特征向量可以恰好组成矩阵V，而特征值组成矩阵D。证明如下：

将上式进行转置，可以得到：

注意到上式运用到了矩阵U是正交阵的性质。因此有：

最终可以得到：

由此可以看出，矩阵V中包含的列向量实际上就是矩阵A^TA的特征向量，而矩阵D中的对角值（矩阵A的奇异值），即为矩阵A^TA的特征值的根号。同理我们可以证明矩阵AA^T的特征向量组成矩阵U。因此，对于矩阵A的奇异值分解，可以转换为计算矩阵A^TA和AA^T的特征值和特征向量。

下面，我们首先计算矩阵A^TA：

然后计算矩阵A^TA的特征值，由

可得

因此有

解得矩阵A^TA的特征值为：

当λ₁=3时：

因此第一个特征向量为：

注意在计算特征向量时需要将其化为单位向量。

当λ₂=1时：

因此第二个特征向量为：

因此有：

下面还需要确定矩阵U，因此继续计算AA^T：

由于

解得矩阵AA^T的特征值为：

对应的特征向量为：

需要注意的是，虽然矩阵AA^T的特征向量构成了矩阵U，但是其排列顺序并不确定，由于

代入矩阵A、矩阵V和矩阵D可以得到：

比较上式和矩阵AA^T的特征向量，可以推导出矩阵U为：

上述内容展示了对矩阵A进行奇异值分解的过程，在MATLAB中，矩阵的奇异值分解也可以直接通过svd函数实现。

在确定了矩阵A的奇异值分解后，则矩阵A的加号逆可表示为：

下面首先以一个线性拟合问题来说明用广义逆矩阵理论解决最小二乘问题的便捷性。取待拟合的数据点为：

设拟合数据点的表达式为：

为使得拟合函数取得对离散数据点的最佳平方逼近，需要设待拟合的数据点为(x_i, y_i) (i=1,2,3)，并构造误差函数：

根据最小二乘法的定义，拟合表达式应使得误差取得极小值，因此必须满足：

通过移项可以得到：

写为矩阵形式为：

代入数据后可解得：

因此拟合的表达式为：

基于广义逆矩阵理论的求解方法

如果采用广义逆矩阵理论，则求解过程将变得非常简单。我们只需直接写出矛盾方程组：

将其整理为矩阵形式：

通过MATLAB的pinv函数可以直接计算得到矩阵A的加号逆为：

因此拟合函数的系数矩阵为：

则其系数矩阵为：

通过MATLAB不难验证，上述系数矩阵与通过拟合工具箱计算得到的结果完全一致。

本文介绍了广义逆矩阵的一些基本概念及定理，并重点介绍了广义逆矩阵的计算方法及其在线性最小二乘问题中的应用。基于广义逆矩阵理论，可以使得求解线性最小二乘问题的过程规范化，只需建立相应的矛盾方程组，即可直接计算其最小二乘解。