许多工程问题中都需要求解线性方程组，最小二乘估计是求解线性方程组最常用的方法，它满足最优线性无偏性，因此又称为无偏估计。但当方程组的系数矩阵呈现病态时，由最小二乘法所得的参数估值往往误差很大甚至完全失真，这称为病态最小二乘问题^[1]。例如，考虑一个常见的线性估计模型：

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{A}}\mathit{\boldsymbol{x}} $

(1)

其中，y为m×1维观测向量(含有观测误差)，A为m×n维系数矩阵(列满秩)，x为n×1维未知参数向量。方程(1)两边乘以A^T可得法方程为：

$ z = \mathit{\boldsymbol{Bx}} $

(2a)

$ \mathit{\boldsymbol{B}} = {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}} $

(2b)

$ z = {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{y}} $

(2c)

则由方程(2)可得x的最小二乘估计为：

$ {\mathit{\boldsymbol{x}}_{{\rm{LSE}}}} = {\mathit{\boldsymbol{B}}^{ - 1}}z $

(3)

由于工程问题的复杂性，矩阵B往往呈现病态或严重病态，即矩阵B的条件数很大或矩阵B近似为奇异矩阵，此时，由方程(3)所得的最小二乘估计x_LSE将严重偏离x的真值，必须寻求x的更好估计。

1 岭估计和广义岭估计

为了解决病态最小二乘问题，许多学者开展了岭估计^[2]和广义岭估计^[3]研究。从矩阵方程的角度，岭估计的基本思想是：在矩阵B中增加一个扰动矩阵k I，以降低系数矩阵的条件数，从而获得更加稳定的估值。其计算公式为：

$ {\mathit{\boldsymbol{x}}_{{\rm{RE}}}} = {(\mathit{\boldsymbol{B}} + k\mathit{\boldsymbol{I}})^{ - 1}}z $

(4)

式中，x_RE为x的岭估计，I为与矩阵B同维的单位矩阵，系数k称为岭参数。如果把方程(4)中的扰动矩阵k I更换为一般的对角矩阵Λ，则可得广义岭估计x_GRE的计算公式：

$ {\mathit{\boldsymbol{x}}_{{\rm{GRE}}}} = {(\mathit{\boldsymbol{B}} + \mathit{\boldsymbol{ \boldsymbol{\varLambda} }})^{ - 1}}z $

(5a)

$ \mathit{\boldsymbol{ \boldsymbol{\varLambda} }} = \left[ {\begin{array}{*{20}{c}} {{k_1}}&{}&{}&{}\\ {}&{{k_2}}&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&{{k_n}} \end{array}} \right] $

(5b)

显然，对比式(3)、式(4)和式(5)可知，当k₁=k₂=…=k_n=k时，广义岭估计退化为岭估计；当k₁=k₂=…=k_n=0时，岭估计退化为最小二乘估计。由于方程(4)和方程(5)并非线性方程(2)的无偏估计(例如，方程(4)本质上是z－kIx = Bx的无偏估计，显然并非方程(2a)的无偏估计)，因此，岭估计和广义岭估计均属于有偏估计。目前，岭估计和广义岭估计的研究重点是，如何选择合适的岭参数k，以使得所得的x_RE最接近x的真值，许多学者就此开展了相关研究^[4-9]。

2 基于Neumann级数的有偏估计

为了更好地探究岭估计的本质，统一各种形式的岭估计公式，并进一步统一有偏估计和无偏估计，本文从Neumann级数的角度，推导出x的一系列有偏估计形式，可以清楚地说明现有相关公式的本质联系，为工程中寻求x的最优估计建立理论基础，同时也是Neumann级数的一次新应用。其推导过程如下。

首先，将方程(3)等价变形为：

$ \begin{array}{l} \mathit{\boldsymbol{x}} = {\mathit{\boldsymbol{B}}^{ - 1}} \cdot z = {(\mathit{\boldsymbol{B}} + \mathit{\boldsymbol{K}} - \mathit{\boldsymbol{K}})^{ - 1}}z = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {((\mathit{\boldsymbol{B}} + \mathit{\boldsymbol{K}}) - \mathit{\boldsymbol{K}})^{ - 1}} \cdot z \end{array} $

(6)

式中，矩阵K是扰动矩阵，可以是前述的对角矩阵，也可以是任意的其他矩阵，只要能够降低条件数即可：

$ \rm{cond}(\mathit{\boldsymbol{B}} + \mathit{\boldsymbol{K}}) < cond(\mathit{\boldsymbol{B}}) $

(7)

式中，cond表示条件数。对方程(6)利用Neumann级数^[10-11]展开得到：

$ \begin{array}{l} \mathit{\boldsymbol{x}} = {(\mathit{\boldsymbol{B}} + \mathit{\boldsymbol{K}})^{ - 1}}[\mathit{\boldsymbol{I}} + \mathit{\boldsymbol{K}}{(\mathit{\boldsymbol{B}} + \mathit{\boldsymbol{K}})^{ - 1}} + \\ \mathit{\boldsymbol{K}}{(\mathit{\boldsymbol{B}} + \mathit{\boldsymbol{K}})^{ - 1}}\mathit{\boldsymbol{K}}{(\mathit{\boldsymbol{B}} + \mathit{\boldsymbol{K}})^{ - 1}} + \cdots ] \cdot z \end{array} $

(8)

根据Neumann级数的收敛条件，方程(8)的收敛条件为：

$ \frac{{\left\| \mathit{\boldsymbol{K}} \right\|}}{{\left\| {\mathit{\boldsymbol{B}} + \mathit{\boldsymbol{K}}} \right\|}} < 1 $

(9)

式中，‖·‖表示矩阵的范数。方程(8)即为基于Neumann级数的x估值公式，有着重要意义。它包含现有的最小二乘估计、岭估计和广义岭估计公式，也建立了有偏估计和无偏估计的本质联系，并可以据此引申出一系列新的有偏估计形式。要点如下。

1) 在方程(8)中，若忽略所有高阶项，只保留第一项可得到：

$ \boldsymbol{x} = {\left( {\boldsymbol{B} + \boldsymbol{K}} \right)^{ - 1}}\boldsymbol{z} $

(10)

对比方程(5a)和方程(10)可以发现，如果扰动矩阵K取为对角矩阵Λ，则方程(10)就转化为广义岭估计公式；进一步，如果K取为k I，则方程(10)就转化为岭估计公式。因此，岭估计和广义岭估计本质上就是方程(8)中取第一项的结果。因此，也可以认为岭估计实际为一级有偏估计。

2) 由于方程(8)是由方程(3)等价变形得到的，而方程(3)是x的无偏估计，若方程(8)中只取有限项，则必然是x的有偏估计。这也说明，有偏估计本质上是从无偏估计中取出一部分，或者更加具体一点，岭估计其实是从最小二乘估计中取出了一部分。另外，如果所取的扰动矩阵K满足收敛条件(9)，则有偏估计可以无限逼近无偏估计。

3) 如前所述，方程(8)中如果只取第一项，则可得到常见的岭估计或广义岭估计；如果方程(8)中保留更多的项，则可得到一系列新的有偏估计形式。比如，令Q = K(B+K)^－1，则方程(8)中取前两项和取前3项所得的有偏估计公式分别为：

$ {x^{'}} = {(\mathit{\boldsymbol{B}} + \mathit{\boldsymbol{K}})^{ - 1}} \cdot [\mathit{\boldsymbol{I}} + \mathit{\boldsymbol{Q}} ] \cdot z $

(11)

$ {x^{''}} = {(\mathit{\boldsymbol{B}} + \mathit{\boldsymbol{K}})^{ - 1}} \cdot [\mathit{\boldsymbol{I}} + \mathit{\boldsymbol{Q}} + {\mathit{\boldsymbol{Q}}^2}] \cdot z $

(12)

方程(11)可以称为二级有偏估计，方程(12)可以称为三级有偏估计。接下来的数值算例分析显示，在这一系列的有偏估计形式中，存在着比岭估计或广义岭估计更接近于x真值的解。

3 算例

以文献[12]所用的病态平差模型y = Ax为例，其中系数矩阵A和观测向量y的真值为：

$ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} 2&{ - 5}&1&1&{ - 9.5}\\ { - 2}&4&1&{ - 1.05}&{8.5}\\ { - 2}&1&1&{ - 1}&{2.4}\\ { - 1}&{2.5}&4&{ - 0.5}&7\\ { - 1}&{3.2}&4&{ - 0.5}&{8.4}\\ 1&1&{ - 3}&{0.4}&{0.49}\\ 3&7&{ - 3}&{1.5}&{12.7}\\ 5&{ - 1}&{ - 2}&{2.5}&{ - 3}\\ 4&2&{ - 2}&{2.01}&3\\ 4&3&{ - 2}&2&5 \end{array}} \right], \mathit{\boldsymbol{y}} = \left[ {\begin{array}{*{20}{c}} { - 10.5}\\ {10.45}\\ {1.4}\\ {12}\\ {14.1}\\ { - 0.11}\\ {21.2}\\ {1.5}\\ {9.01}\\ {12} \end{array}} \right] $

(13)

未知参数向量x的真值为x _true=[1 1 1 1 1]^T。显然，该方程组存在严重病态，因为其相应法方程中的系数矩阵B = A^T A的条件数为：cond(B)=1.289 2×10⁵。如果观测向量y无误差，则由方程(13)的最小二乘估计即可获得x的真值。现考虑观测向量y中存在误差的情况，误差的模拟方式为：使用Matlab对观测向量y中的每一个元素添加一个服从正态分布(均值为0，方差为0.01)的随机数。不失一般性，假设含有误差的观测向量为y^c=[－10.178 1, 10.640 4, 1.372 4, 12.050 0, 13.607 7, －0.112 2, 20.942 2, 1.554 0, 9.328 6, 12.112 3]^T。此时，当取不同的扰动矩阵K时，上述方程的最小二乘估计x_LSE、岭估计x_RE(也是一级有偏估计)、二级有偏估计x′和三级有偏估计x″分别列于表 1~5中，其中e_Δx表示x的估计值与真值之间的相对偏差，即

$ {\mathit{\boldsymbol{e}}_{\Delta x}} = \frac{{{{\left\| {{\mathit{\boldsymbol{x}}_{\rm{estimate}}} - {\mathit{\boldsymbol{x}}_{\rm{true}}}} \right\|}_2}}}{{{{\left\| {{x_{\rm{true}}}} \right\|}_2}}} $

(14)

由表 1~5可见，当观测向量y中含有误差时，由最小二乘估计所得结果和真值偏差较大，而各级有偏估计的结果和真值更加接近。就本例而言，二级和三级有偏估计解都比岭估计更接近于真值。上述结果说明，本文基于Neumann级数的有偏估计方法是合理可行的。

接下来分析k值的选取对计算结果的影响规律。不失一般性，由表 1和表 5对比可见，当k=0.2时的各级有偏估计解均比当k=20时的解精确，根据现有文献中选择岭参数的经验可知，当k值越接近于最优岭参数时，其计算结果越准确。因此，本文方法在应用时，可以先用现有的岭参数选择方法来选择一个合适的k值，然后再使用本文所提的各级有偏估计公式，这样能迅速获得更加准确的计算结果。需要指出的是，即使所取的k值远离最优岭参数，比如上述算例中k=20、30、50时，通过增加方程(8)中所保留的级数阶数，也能获得与k=0.2、0.5时精度相当的解。表 6给出了取不同k值时获得精度相当的解所需要的级数阶数。由表 6可见，随着k值的增大，要获得精度相当的解所需要保留的级数阶数也增大，如k=50时，其第760级有偏估计解才与k=0.2时第3级有偏估计解的精度相当。但这也同时说明，即使所取的k值偏大，也总能通过增加级数阶数来获得较好的计算结果。需要注意的是，无论对于哪一个k值，保留的级数阶数并非越多越好，因为随着级数阶数的无限增加，方程(8)将无限逼近最小二乘估计，因此所保留的级数阶数也有一个最优值。表 7给出了当k=0.5时前18级有偏估计解与真值之间的相对误差变化情况，由表 7可见，与真值最接近的解是第11级和第12级有偏估计，因为其对应的相对误差最小(均为0.301 3)。显然，随着级数阶数的增加，解的精度呈现先逐渐提高后又逐渐降低的趋势，当k=0.5时，最优的有偏估计为第11和12级有偏估计。

上述病态最小二乘问题只考虑了观测向量中含有误差的情况，工程实践中，线性模型的系数矩阵里往往也含有误差，相应的病态问题解算中需要同时考虑系数矩阵和观测向量的误差，这种问题称为病态整体最小二乘问题^[1,12-14]。本文方法也可用于病态整体最小二乘问题的解算。为了便于比较，接下来的讨论采用文献[14]所模拟的病态平差模型，其系数矩阵和观测向量由方程(13)中的数据添加随机误差得到，所添加的随机误差均服从均值为0、方差为0.01的正态分布。加入随机误差后的系数矩阵和观测向量具体为：

$ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {1.8812}&{ - 5.1186}&{1.0129}&{1.0806}&{ - 9.5331}\\ { - 2.2202}&{3.8944}&{1.0656}&{ - 1.0268}&{8.4156}\\ { - 1.9014}&{1.1472}&{0.8832}&{ - 1.099}&{2.4498}\\ { - 1.0519}&{2.5056}&{3.9539}&{ - 0.366}&{7.1488}\\ { - 0.9673}&{3.0783}&{3.9738}&{ - 0.471}&{8.3454}\\ {1.0234}&{0.9959}&{ - 3.1213}&{0.5479}&{0.4053}\\ {3.0021}&{6.8872}&{ - 3.1319}&{1.6138}&{12.675}\\ {4.8996}&{ - 1.1349}&{ - 1.9069}&{2.4316}&{ - 2.9337}\\ {3.9053}&{1.9739}&{ - 1.9989}&{1.8808}&{2.9146}\\ {3.9626}&{3.0953}&{ - 2.0645}&{1.9927}&{4.8799} \end{array}} \right], \mathit{\boldsymbol{y}} = \left[ {\begin{array}{*{20}{c}} { - 10.512}\\ {10.443}\\ {1.4485}\\ {11.94}\\ {14.085}\\ { - 0.1535}\\ {21.192}\\ {1.6535}\\ {8.9494}\\ {11.865} \end{array}} \right] $

(15)

该模型法矩阵的条件数为2.083 7×10⁴，病态性严重。文献[14]中给出了各种方法的解算结果，如表 8所示，其中TLS表示整体最小二乘法，LSE表示最小二乘法，RE表示岭估计方法，VOM表示虚拟观测解法，估计值与真值之间的偏差e采用绝对偏差(即e=‖ x_estimate－ x_true‖₂)。

采用本文方法，当k取10时，前5级有偏估计结果如表 9所示。由表 9可见，第3级有偏估计解的精度最好。对比表 8和表 9可见，本文方法第2~5级解的精度略高于表 8中各种方法的解算精度，进一步说明本文方法是合理可行的。

4 结语

本文提出一种基于Neumann级数的有偏估计方法，用于解决病态最小二乘问题。所提方法是Neumann级数的一次新应用，它包含现有的最小二乘估计、岭估计和广义岭估计公式，也建立了有偏估计和无偏估计的本质联系，并可以据此引申出一系列新的有偏估计形式。从本文的数值算例结果来看，在这一系列的有偏估计形式中，存在着比岭估计或广义岭估计更接近于真值的解。算例结果初步说明了本文方法的价值，为解决病态最小二乘问题提供了一条新途径。