超定线性系统最小二乘解的可信验证

The Verification of the Least Square Solution of an Overdetermined Linear System

摘要: 给定一个系数矩阵为病态矩阵的超定线性系统，本文利用区间算法和增广矩阵技术，设计算法输出给定超定线性系统的微小摄动区间系统，算法保证在输出的区间系统中存在系数阵精确秩亏的超定线性系统，算法输出该系数阵精确秩亏的超定线性系统的最小二乘解的高精度近似解和其相应的可信误差界。

Abstract: Given an overdetermined linear system with ill-conditioned coefficient matrix, this paper mainly uses interval algorithm and augmented matrix technology to design an algorithm to output a minimally perturbed interval system for the given overdetermined linear system. The algorithm guarantees that there is an overdetermined linear system with a coefficient matrix of exactly rank deficiency in the output interval system. The algorithm also outputs the high-precision approximate solution of the least square solution of the overdetermined linear system and corresponding verified error bound.

1. 引言

线性系统是计算科学领域中一种很重要的数学模型。线性系统理论在航空航天、化学化工、机械研究等技术领域中都有着广泛的应用实例。超定线性系统是指方程个数大于变量个数的线性系统。研究表明，在数字信号处理和其他面向定量的学科中有许多问题最终都归结为超定线性系统求解的计算问题，最小二乘拟合算法的实质就是对超定线性系统进行求解。

多年来，对于超定线性系统求解问题的研究始终吸引着国内外的专家学者。Levitan等人 [1] 研究了针对特定抽象范数求解超定线性系统最优解的解决方案。Barroda和Roberts [2] 提出了能够解决线性数据拟合问题的有效算法，该算法是对线性规划单纯形法的一种改进，适用于 $l_1$ 问题的初步表述。Cadzow [3] 提出了超定线性系统最小1、2和 $\infty$ 范数近似解的解决方案。Demmel等人 [4] 研究了超定最小二乘问题精确迭代细化的算法、误差范围和数值结果，将线性系统优化算法应用于超定线性系统最小二乘问题的增强线性系统公式。Novoselac和Pavic [5] 研究了 $p$ 范数中超定线性方程组的解性质，从而得到了该问题具有回归、尺度和仿射等变性质。Popa [6] 提出了超定线性系统最小二乘解的数值计算的两种迭代算法，这些算法在Kaczmarz迭代算法的基础上，可以收敛到最小二乘解的序列。Pelckmans等人 [7] 对超定线性系统的噪声观测进行了稀疏估计，提出了一种受高斯噪声干扰的线性方程估计未知参数的方法。Zhijun Zhang等人 [8] 为解决时变的超定问题，提出并设计了一种新颖的变参数收敛微分神经网络，通过研究表明，这种神经网络具有超指数收敛速度的最小二乘解。

近年来，很多学者着手研究超定系统最小二乘解的可信验证。Rump [9] 提出了一种不依赖于正规方程且适用于稀疏矩阵的新算法，该算法用于计算系数阵列满秩的最小二乘问题的最小二乘解和欠定线性系统的最小二范数解。Miyajima和Shinya [10] 在考虑了所有可能的舍入误差的情况下，给出了系数阵为列满秩的欠定系统最小二范数解和超定系统最小二乘解的验证算法，阐明了所获得的误差范围的大小与系数矩阵的条件数相关。Rump [11] 提出的快速算法的误差界是标量，而Miyajima提出的算法的误差界是向量。Miyajima [12] 总结提出了一种用于矩阵方程 $AXB=C$ 的最小二乘解的高效能算法。

本文考虑系数矩阵为病态矩阵的超定线性系统的可信计算。计算给定超定线性系统的微小摄动区间系统，算法保证在输出的区间系统中存在系数阵精确秩亏的超定线性系统，并计算该系数阵精确秩亏的超定线性系统的最小二乘解的可信区间解。

2. 符号和预备知识

令 $ℝ$ 表示实数集合。设 $A\in {ℝ}^{m\times n}$， $A(:)$ 表示由矩阵 $A$ 的列合并得到的长列向量。对于满足条件 $1\le i_1\le i_2\le m$， $1\le i_1\le i_2\le m$ 的正整数 $i_{\text{1}},i_{\text{2}},j_{\text{1}},j_{\text{2}}$， $A_{i_1:i_2,:}$ 表示由矩阵 $A$ 的第 $i_1$ 行到第 $i_2$ 行所构成的子矩阵， $A_{:,j_1:j_2}$ 表示由矩阵 $A$ 的第 $j_1$ 列到第 $j_2$ 列所构成的子矩阵， $A_{i_1:i_2,j_1:j_2}$ 表示由矩阵 $A$ 的第 $i_1$ 行到第 $i_2$ 行，第 $j_1$ 列到第 $j_2$ 列所构成的子矩阵。令 $O_{m,n}$ 表示 $m\times n$ 阶零矩阵， $I_{m\times m}$ 表示 $m\times m$ 阶单位阵。

令 $Iℝ$ 表示区间集合 $\left\{\left[x,\bar{x}\right]:\underset{\_}{x},\bar{x}\in ℝ,x\le \bar{x}\right\}$，用 $mid\left[\underset{\_}{x},\bar{x}\right]$ 表示区间 $\left[\underset{\_}{x},\bar{x}\right]$ 的中点，用length表示区间长度，即 $\text{length}\left(\left[x,\hat{x}\right]\right)=\bar{x}-\underset{\_}{x}$。区间向量 $X={\left(X_1,X_2,\cdots ,X_n\right)}^{\text{T}}\in \mathbb{I}{ℝ}^n$ 是分量为区间的 $n$ 维向量，令 $\mathrm{int}\left(X\right)$ 表示它的内部。区间矩阵 $A$ 是分量为区间的矩阵。给定区间矩阵 $A\in \mathbb{I}{ℝ}^{n\times n}$，若区间矩阵 $A$ 包含的所有实矩阵都是非奇异的，则称区间矩阵 $A$ 非奇异。给定区间矩阵 $A\in \mathbb{I}{ℝ}^{m\times n}$，若区间矩阵 $A$ 包含的所有实矩阵都是满秩的，则称区间矩阵 $A$ 满秩。对于 $A\in \mathbb{I}{ℝ}^{m\times n}$， $mid\left(A\right)$ 是由矩阵 $A$ 中每个分量中点对应的值构成的实矩阵，令 $\text{length}\left(A\right)=\text{max}\left\{\text{length}\left(A_{i,j}\right),1\le i\le m,1\le j\le n\right\}$。

定义1：对于给定的阈值 $\delta >0$，若所有与 $A$ 的2-范数距离不超过 $\delta$ 的 $m\times n$ 阶矩阵的最小秩为 $r$，即

$r=\underset{\begin{array}{l}{\Vert A-B\Vert }_2\le \delta \\ B\in {ℝ}^{m\times n}\end{array}}{\min }\left\{rank\left(B\right)\right\}$

则称 $A$ 的数值 $\delta$ -秩为 $r$，记为 $ran{k}_{\delta }\left(A\right)=r$。

注1：给定矩阵 $A\in {ℝ}^{m\times n}$，其中 $m>n$，对于阈值 $\delta >0$，如果矩阵 $A$ 的奇异值 ${\sigma }_1\left(A\right),\cdots ,{\sigma }_n\left(A\right)$ 满足

${\sigma }_1\left(A\right)\ge \cdots \ge {\sigma }_r\left(A\right)>\delta \ge {\sigma }_{r+1}\left(A\right)\ge \cdots \ge {\sigma }_n\left(A\right)$,

则 $ran{k}_{\delta }\left(A\right)=r$。

引理1：给定矩阵 $A,B\in {ℝ}^{n\times n}$，若 ${\Vert I_n-BA\Vert }_p\le \alpha <1$，其中 $p\in \left\{1,2,\cdots ,\infty \right\}$，则 $A,B$ 均是非奇异矩阵。

定理1：对于给定的区间矩阵 $A\in \mathbb{I}{ℝ}^{n\times n}$ 和区间向量 $b\in \mathbb{I}{ℝ}^n$，向量 $b$ 满足 $0\notin b$，若INTLAB工具箱中verifylss函数成功输出区间向量 $X\subset \mathbb{I}{ℝ}^n$，则区间矩阵 $A$ 非奇异，且 $X$ 满足条件

$\left\{x\in {ℝ}^n:Ax=b\forall A\in A,\forall b\in b\right\}\subseteq X$ (1)

定理2：设矩阵 $A\in \mathbb{I}{ℝ}^{m\times n}$， $m>n$，向量 $b\in \mathbb{I}{ℝ}^m$ 满足 $0\notin b$，当输入区间矩阵 $\left(\begin{array}{cc}A& I_m\\ O_{n,n}& A^{\text{T}}\end{array}\right)$ 和区间向

量 ${\left(b,O_{n,1}\right)}^{\text{T}}$ 时，若INTLAB工具箱中verifylss函数成功地输出区间向量，则区间矩阵 $A$ 满秩。

定理3：设 $f:{ℝ}^n\to {ℝ}^n$ 并且 $f={\left(f_1,f_2,\cdots ,f_n\right)}^{\text{T}}$，其中 $f_1,\cdots ,f_n$ 为连续可微函数。当 $\tilde{x}\in {ℝ}^n$ 且 $X\in I{ℝ}^n$， $0\in X$ 时，令 $f^{\prime }\left(\tilde{x}+X\right)$ 表示系统在区间 $\tilde{x}+X$ 上的雅可比矩阵。设区间矩阵 $M\in I{ℝ}^{n\times n}$ 满足条件 $f^{\prime }\left(\tilde{x}+X\right)\subseteq M$。若对矩阵 $R\in {ℝ}^{n\times n}$，有

$-Rf\left(\tilde{x}\right)+\left(I_n-RM\right)X\subseteq \mathrm{int}\left(X\right),$ (2)

则存在唯一的 $\hat{x}\in \tilde{x}+X$ 使得 $f\left(\hat{x}\right)=0$。

注2：INTLAB工具箱中verifylss函数利用定理3来验证非线性系统的解。

3. 主要结论

设 $A\in {ℝ}^{m\times n}$，其中 $m>n$，令 $A=U\Sigma V^{\text{T}}$ 表示 $A$ 的奇异值分解。若 $ran{k}_{\delta }\left(A\right)=r$，则 $U_{:,r+1:n}$ 和 $V_{:,r+1:m}$ 的列向量分别张成矩阵 $A$ 的近似 $\delta$ -左零空间和近似 $\delta$ -右零空间。定义变量矩阵 $E$

$E=\left(\begin{array}{cccc}{e}_{1,1}& e_{1,2}& \cdots & e_{1,n}\\ e_{2,1}& e_{2,2}& \cdots & e_{2,n}\\ \cdots & \cdots & \cdots & \cdots \\ e_{m,1}& e_{m,2}& \cdots & e_{m,n}\end{array}\right)$ (3)

为了简便起见，用符号 $e$ 表示向量 $E(:)$，用符号 $A\left(e\right)$ 表示矩阵 $A+E$。

引理2：令 $A\in {ℝ}^{m\times n}$，其中 $m>n$，且 $ran{k}_{\delta }\left(A\right)=r$，则矩阵

$\left(\begin{array}{cc}A\left(e\right)& U_{:,r+1:m}\\ V_{:,r+1:n}^{\text{T}}& O_{n-r,n-r}\end{array}\right)$

当 $e=0$ 时非奇异。

引理3：令 $A\in {ℝ}^{m\times n}$， $m>n$，且 $ran{k}_{\delta }\left(A\right)=r$，则矩阵 $H(\; e\; )$

$H\left(e\right)=\left(\begin{array}{cccc}A\left(e\right)& U_{:,r+1:n}& -I_m& O_{m,n-r}\\ V_{:,r+1:n}^{\text{T}}& O_{n-r,n-r}& O_{n-r,m}& -I_{n-r}\\ O_{n,n}& O_{n,n-r}& A{\left(e\right)}^{\text{T}}& V_{:,r+1:n}^{\text{T}}\\ O_{n-r,n}& O_{n-r,n-r}& U_{:,r+1:n}^{\text{T}}& O_{n-r,n-r}\end{array}\right)$ (4)

当 $e=0$ 时非奇异。

引入如下线性系统

$\left(\begin{array}{cccc}A\left(e\right)& U_{:,r+1:n}& -I_m& O_{m,n-r}\\ V_{:,r+1:n}^{\text{T}}& O_{n-r,n-r}& O_{n-r,m}& -I_{n-r}\\ O_{n,n}& O_{n,n-r}& A{\left(e\right)}^{\text{T}}& V_{:,r+1:n}^{}\\ O_{n-r,n}& O_{n-r,n-r}& U_{:,r+1:n}^{\text{T}}& O_{n-r,n-r}\end{array}\right)\left(\begin{array}{l}C\left(e\right)\\ F\left(e\right)\\ Y\left(e\right)\\ Z\left(e\right)\end{array}\right)=\left(\begin{array}{c}{O}_{m,n-r}\\ I_{n-r}\\ O_{n,n-r}\\ O_{n-r,n-r}\end{array}\right).$ (5)

其中矩阵 $C\left(e\right)$， $F\left(e\right)$， $Y\left(e\right)$ 和 $Z\left(e\right)$ 的行维数分别时是 $n$， $n-r$， $m$ 和 $n-r$。

定义矩阵

$\tilde{A}\left(e\right)=A\left(e\right)-Y\left(e\right){\left(C^{\text{T}}\left(e\right)C\left(e\right)\right)}^{-1}{C}^{\text{T}}\left(e\right).$ (6)

命题1：假设向量 $\hat{e}\in {ℝ}^{n^2}$ 满足如下条件：

(1) ${\Vert \hat{e}\Vert }_{\infty }<\frac{1}{\left(m+n\right)\sqrt{m+n}{\Vert H{\left(0\right)}^{-1}\Vert }_2}$，

(2) $F\left(\hat{e}\right)=O_{n-r,n-r}$，

(3) 矩阵 $Z\left(\hat{e}\right)+I_{n-r}$ 非奇异，

(4) 矩阵

$\left(\begin{array}{l}\tilde{A}\left(\hat{e}\right)\\ C{\left(\hat{e}\right)}^{\text{T}}\end{array}\right)$. (7)

满秩，则 $rank\left(\tilde{A}\left(\hat{e}\right)\right)=r$。

证明：易知

$\begin{array}{c}{\Vert I_{m+n}-H{\left(0\right)}^{-1}H\left(\hat{e}\right)\Vert }_2\le {\Vert H{\left(0\right)}^{-1}\Vert }_2{\Vert H\left(0\right)-H\left(\hat{e}\right)\Vert }_2\\ \le \sqrt{m+n}{\Vert H{\left(0\right)}^{-1}\Vert }_2{\Vert H\left(0\right)-H\left(\hat{e}\right)\Vert }_{\infty }\\ \le \left(m+n\right)\sqrt{m+n}{\Vert H{\left(0\right)}^{-1}\Vert }_2{\Vert \hat{e}\Vert }_{\infty }<1.\end{array}$

故矩阵 $H\left(\hat{e}\right)$ 非奇异。若 $F\left(\hat{e}\right)=O_{n-r,n-r}$，则 $A\left(\hat{e}\right)C\left(\hat{e}\right)=Y\left(\hat{e}\right)$。由(6)式可得 $\tilde{A}\left(\hat{e}\right)C\left(\hat{e}\right)=O_{n-r,m}$，断言 $C\left(\hat{e}\right)$ 满秩。事实上，如果 $C\left(\hat{e}\right)$ 的列向量线性相关，则存在一个非零向量 $\alpha \in {ℝ}^{m-r}$ 使得 $C\left(\hat{e}\right)\alpha =0$。由于 $V_{:,r+1:n}^{T}C\left(\hat{e}\right)-Z\left(\hat{e}\right)=I_{n-r}$，所以 $Z\left(\hat{e}\right)\alpha +\alpha =0$，这与矩阵 $Z\left(\hat{e}\right)+I_{n-r}$ 非奇异相矛盾。假设 $rank\left(\tilde{A}\left(\hat{e}\right)\right)<r$，则存在一个向量 $\beta \in null\left(\tilde{A}{\left(\hat{e}\right)}^{\text{T}}\right)$ 使得矩阵 $\left(C\left(\hat{e}\right)\beta \right)$ 满秩。定义非零向量

$\gamma =\left(C\left(\hat{e}\right)\beta \right)\left(\begin{array}{l}{\left(C{\left(\hat{e}\right)}^{\text{T}}C\left(\hat{e}\right)\right)}^{-1}C{\left(\hat{e}\right)}^{\text{T}}\beta \\ -1\end{array}\right),$

通过计算可知

$\left(\begin{array}{l}\tilde{A}\left(\hat{e}\right)\\ C\left(\hat{e}\right){}^{\text{T}}\end{array}\right)\gamma =\text{0}.$

上式与矩阵(7)满秩相矛盾。 □

由引理3可知，存在零向量的一个邻域，使得对这个邻域内的任何向量 $e$，矩阵 $H\left(e\right)$ 均是非奇异的。因此，在这个邻域内，系统(5)的解的每一项对每个变量均可导。对系统(5)的两边关于变量求导，可得如下的线性系统

$H\left(e\right)\left(\begin{array}{l}\frac{\partial C\left(e\right)}{\partial e_{i,j}}\\ \frac{\partial F\left(e\right)}{\partial e_{i,j}}\\ \frac{\partial Y\left(e\right)}{\partial e_{i,j}}\\ \frac{\partial Z\left(e\right)}{\partial e_{i,j}}\end{array}\right)=\left(\begin{array}{l}{O}_{i-1,n-r}\\ -C\left(e\right){}_{j,:}\\ O_{m+n-r-i+j-1,n-r}\\ -Y\left(e\right){}_{i,:}\\ O_{2n-r-j,n-r}\end{array}\right),1\le i\le m,1\le j\le n,$ (8)

$H\left(e\right)\left(\begin{array}{l}\frac{{\partial }^{2}C\left(e\right)}{\partial e_{i_1,j_1}\partial e_{i_2,j_2}}\\ \frac{{\partial }^{2}F\left(e\right)}{\partial e_{i_1,j_1}\partial e_{i_2,j_2}}\\ \frac{{\partial }^{2}Y\left(e\right)}{\partial e_{i_1,j_1}\partial e_{i_2,j_2}}\\ \frac{{\partial }^{2}Z\left(e\right)}{\partial e_{i_1,j_1}\partial e_{i_2,j_2}}\end{array}\right)=\left(\begin{array}{l}{O}_{i_2-1,n-r}\\ {\left(-\frac{\partial C\left(e\right)}{\partial e_{i_1,j_1}}\right)}_{j_2,:}\\ O_{m+n-r-i_2+j_2-1,n-r}\\ {\left(-\frac{\partial Y\left(e\right)}{\partial e_{i_1,j_1}}\right)}_{i_2,:}\\ O_{2n-r-j_2,n-r}\end{array}\right)+\left(\begin{array}{l}{O}_{i_1-1,n-r}\\ {\left(-\frac{\partial C\left(e\right)}{\partial e_{i_2,j_2}}\right)}_{j_1,:}\\ O_{m+n-r-i_1+j_1-1,n-r}\\ {\left(-\frac{\partial Y\left(e\right)}{\partial e_{i_2,j_2}}\right)}_{i_1,:}\\ O_{2n-r-j_1,n-r}\end{array}\right),1\le i\le m,1\le j\le n.$ (9)

考虑如下约束优化问题

$\begin{array}{l}\min \frac{1}{2}{\Vert e\Vert }_2^2\\ s.t.F\left(e\right)=O_{n-r,n-r}\end{array}$ (10)

引入拉格朗日乘子

$\lambda ={\left({\lambda }_{1,1},{\lambda }_{1,2},\cdots ,{\lambda }_{n-r},{\lambda }_{1,2},{\lambda }_{2,2},\cdots ,{\lambda }_{n-r,2},\cdots ,{\lambda }_{1,n-r},{\lambda }_{2,n-r},\cdots ,{\lambda }_{n-r,n-r}\right)}^{\text{T}}$

定义拉格朗日函数

$g\left(e,\lambda \right)=\frac{1}{2}{\Vert e\Vert }_2^2+{\displaystyle \underset{j=1}{\overset{n-r}{\sum }}{\displaystyle \underset{j=1}{\overset{n-r}{\sum }}{\lambda }_{i,j}F{\left(e\right)}_{i,j}}}.$ (11)

经计算可知 $g\left(e,\lambda \right)$ 的梯度为

$g^{\prime }\left(e,\lambda \right)=\left(\begin{array}{l}e+{\displaystyle {\sum }_{i=1}^{n-r}{\displaystyle {\sum }_{j=1}^{n-r}{\lambda }_{i,j}{F}^{\prime }{\left(e\right)}_{i,j}}}\\ F\left(e\right)\end{array}\right).$ (12)

$g\left(e,\lambda \right)$ 的海森矩阵为

$g^{″}\left(e,\lambda \right)=\left(\begin{array}{cc}{I}_{mn}+{{\displaystyle {\sum }_{i=1}^{n-r}{\displaystyle {\sum }_{j=1}^{n-r}{\lambda }_{i,j}{F}^{″}\left(e\right)}}}_{i,j}& F^{\prime }\left(e\right){(:)}^{\text{T}}\\ F^{\prime }\left(e\right)(:)& O_{{\left(n-r\right)}^2,{\left(n-r\right)}^2}\end{array}\right).$ (13)

注3：数值部分利用牛顿迭代法

$\left(\begin{array}{l}\tilde{e}\\ \tilde{\lambda }\end{array}\right)\leftarrow \left(\begin{array}{l}\tilde{e}\\ \tilde{\lambda }\end{array}\right)-g^{″}{\left(\tilde{e},\tilde{\lambda }\right)}^{-1}{g}^{\prime }\left(\tilde{e},\tilde{\lambda }\right)$ (14)

计算 $g^{\prime }\left(e,\lambda \right)=0$ 的高精度近似解 $\left(\tilde{e},\tilde{\lambda }\right)$，其中 $g^{\prime }\left(\tilde{e},\tilde{\lambda }\right)$， $g^{″}\left(\tilde{e},\tilde{\lambda }\right)$ 可由求系统(8) (9)当 $e=\tilde{e}$ 时的解得到。

注4：验证部分利用区间算法计算系统 $g^{\prime }\left(e,\lambda \right)=0$ 近似解的可信误差界。对于任何 $\tilde{E}\in \mathbb{I}{ℝ}^{mn}$，

$\tilde{\Lambda }\in \mathbb{I}{ℝ}^{(n-r)(n-r)}$，利用verifylss函数求解当 $e=\tilde{E}$ 时线性方程组(5) (8) (9)的解，从而得到 $g^{\prime }\left(\tilde{E},\tilde{\Lambda }\right)$ 和

$g^{″}\left(\tilde{E},\tilde{\Lambda }\right)$。利用定理3和verifylss函数的代码，验证部分计算区间向量 $\hat{E}$ 和 $\hat{\Lambda }$，其分别包含实向量 $\hat{e}$ 和 $\hat{\lambda }$，

使得 $g\text{'}\left(\hat{e},\hat{\lambda }\right)=0$。

算法 1

输入 $A\in {ℝ}^{m\times n}$， $m>n$， $b\in {ℝ}^m$， $\delta \in {ℝ}^{+}$ ；

输出 $r\in \mathbb{N}$， $\tilde{A}\left(\hat{e}\right)\in \mathbb{I}{ℝ}^{m\times n}$ 或者“Failure”。

步骤1计算满足条件(1)的正数 $r$。

步骤2计算系统(11)的近似解 $\left(\tilde{e},\tilde{\lambda }\right)$。

步骤3计算包含系统(11)精确解的区间向量 $\left(\hat{E},\hat{\Lambda }\right)$。

步骤4利用(6)式计算区间矩阵 $\tilde{A}\left(\hat{E}\right)$。

步骤5利用定理2验证矩阵 $Z\left(\hat{E}\right)+I_{n-r}$ 是否非奇异，如果失败，则返回“Failure”。

步骤6利用定理3验证矩阵 $\left(\begin{array}{l}\tilde{A}\left(\hat{E}\right)\\ C{\left(\hat{E}\right)}^{\text{T}}\end{array}\right)$ 是否秩亏，如果失败，则返回“Failure”。 □

若算法1成功地输出区间向量 $\hat{E}$，则由命题5可知，存在 $\hat{e}\in \hat{E}$，使得 $rank\left(\tilde{A}\left(\hat{e}\right)\right)=r$。因此，存在

满足条件 $1\le i_1<i_2,\cdots <i_r\le n$ 的正整数序列 $i_1,i_2\cdots ,i_r$，使得区间矩阵 $A{\left(\hat{E}\right)}_{:,i_1,i_2,\cdots ,i_r}$ 满秩，则考虑如下区间线

性系统

$\left(\begin{array}{cc}\tilde{A}{\left(\hat{E}\right)}_{:,i_1,i_2,\cdots ,i_r}& -I_m\\ O_{r,r}& \tilde{A}{\left(\hat{E}\right)}^{\text{T}}{}_{:,i_1,i_2,\cdots ,i_r}\end{array}\right)\left(\begin{array}{l}y\left(\hat{E}\right)\\ \omega \left(\hat{E}\right)\end{array}\right)=\left(\begin{array}{l}b\\ 0\end{array}\right).$ (15)

命题2：设 $r$ 维向量 $x\left(e\right)$ 向量满足如下条件

$x{\left(e\right)}_{:,k}=\{\begin{array}{l}y{\left(e\right)}_{:,j,}当k=i_j,1\le j\le r;\\ 0,其他情况.\end{array}$ (16)

若向量 $\hat{e}\in R^{m\times n}$ 满足 $rank\left(\tilde{A}\left(\hat{e}\right)\right)=r$，则 $x\left(\hat{e}\right)$ 是超定系统

$\tilde{A}\left(\hat{e}\right)x=b$

的最小二乘解。

证明：假设存在向量 $\tilde{x}\in {ℝ}^n$，使得

${\Vert \tilde{A}\left(\hat{e}\right)\tilde{x}-b\Vert }_2\le {\Vert \tilde{A}\left(\hat{e}\right)x\left(\hat{e}\right)-b\Vert }_2.$ (17)

根据 $rank\left(\tilde{A}\left(\hat{e}\right)\right)=r$，可知 $\hat{A}{\left(\hat{e}\right)}_{:r_1,r_2,\cdots ,r_r}$ 满秩。因此，存在实向量 ${\alpha }^{(1)},{\alpha }^{(2)},\cdots ,{\alpha }^{(n-r)}\in {ℝ}^r$，使得

$\tilde{A}{\left(\hat{e}\right)}_{:j_1,j_2,\cdots ,j_{n-r}}=\tilde{A}{\left(\hat{e}\right)}_{:i_1,i_2,\cdots ,i_r}\left({\alpha }^{(1)}{\alpha }^{(2)}\cdots {\alpha }^{(n-r)}\right),$ (18)

其中 $\left\{j_1,j_2,\cdots ,j_{n-r}\right\}=\left\{1,2,\cdots ,n\right\}\backslash \left\{i_1,i_2,\cdots ,i_r\right\}$。显然，区间系统(15)的解包含着下述系统

$\left(\begin{array}{cc}\tilde{A}{\left(\hat{e}\right)}_{:,i_1,i_2,\cdots ,i_r}& -I_m\\ O_{r,r}& \tilde{A}{\left(\hat{e}\right)}^{\text{T}}{}_{:,i_1,i_2,\cdots ,i_r}\end{array}\right)\left(\begin{array}{l}y\left(\hat{e}\right)\\ \omega \left(\hat{e}\right)\end{array}\right)=(\; b\; 0\; )$

的解。由(17) (18)式可知，

${\Vert \tilde{A}\left(\hat{e}\right)\tilde{x}-b\Vert }_2={\Vert \tilde{A}{\left(\hat{e}\right)}_{:i_1,i_2,\cdots ,i_r}\left(\left({\alpha }^{(1)}{\alpha }^{(2)}\cdots {\alpha }^{(n-r)}\right){\tilde{x}}_{j_1,j_2,\cdots ,j_r}+{\tilde{x}}_{i_1,i_2,\cdots ,i_r}\right)-b\Vert }_2\le {\Vert \tilde{A}{\left(\hat{e}\right)}_{:i_1,i_2,\cdots ,i_r}y\left(\hat{e}\right)-b\Vert }_2，$

而 $y\left(\hat{e}\right)$ 是系数矩阵 $\tilde{A}{\left(\hat{e}\right)}_{:i_1,i_2,\cdots ,i_r}$ 和右端向量 $b$ 的超定线性系统的唯一最小二乘解，故

$y\left(\hat{e}\right)={\tilde{x}}_{i_1,i_2,\cdots ,i_r}+\left({\alpha }^{(1)}{\alpha }^{(2)}\cdots {\alpha }^{(n-r)}\right){\tilde{x}}_{j_1,j_2,\cdots ,j_{n-r}}.$

因此 ${\Vert \tilde{A}\left(\hat{e}\right)\tilde{x}-b\Vert }_2={\Vert \tilde{A}\left(\hat{e}\right)x\left(\hat{e}\right)-b\Vert }_2.$ □

算法1

输入 $A\in {ℝ}^{m\times n}$， $m\le n$ ； $b\in {ℝ}^m$， $\delta \in {ℝ}^{+}$ ；

输出 $r\in \mathbb{N}$， $\hat{A}\in \mathbb{I}{ℝ}^{m\times n}$， $\hat{X}\in \mathbb{I}{ℝ}^n$。

步骤1执行算法1。

步骤2利用verifylss函数解线性方程组(15)计算 $y\left(\hat{E}\right)$， $\omega \left(\hat{E}\right)$。

步骤3确定指标集合 $i_1,i_2,\cdots ,i_r$，使得 $\hat{A}{\left(\hat{E}\right)}_{:i_1,i_2,\cdots ,i_r}$ 满秩。

步骤4利用(15) (16)式计算区间向量 $\hat{x}\left(\hat{E}\right)$。

步骤5令 $\hat{A}\leftarrow A\left(\hat{E}\right)$， $\hat{X}\leftarrow \hat{x}\left(\hat{E}\right)$。 □

由命题1和命题2可得如下定理。

定理4：若算法1成功输出 $r\in \mathbb{N}$， $\hat{A}\in \mathbb{I}{ℝ}^{m\times n}$， $\hat{X}\in \mathbb{I}{ℝ}^n$，则存在实矩阵 $\hat{A}\in \hat{A}$，实向量 $\hat{x}\in \hat{X}$ 使得 $\hat{x}$ 是超定线性系统 $\hat{A}x\text{=}b$ 的最小二乘解，且 $rank\left(\hat{A}\right)=r$。

4. 应用实例

例1给定超定线性系统 $Ax=b$，其中，

$A=\left[\begin{array}{ccccc}0.6406& -0.7615& 0.0742& -1.4399& 0.5767\\ -0.5122& 0.9371& -0.9758& -0.0093& -0.0319\\ -0.4263& 0.1856& -0.2452& 0.2856& -0.4537\\ 0.1684& 0.0526& -0.4056& -0.8077& 0.3922\\ 0.4236& -0.5423& 0.3510& -0.4947& 0.2687\\ -0.8300& -0.4306& -0.9047& -1.0113& -1.0731\end{array}\right]$

$b={\left(\begin{array}{cccccc}2.0790& -0.8294& -0.6716& 0.7956& 1.0071& 0.3704\end{array}\right)}^{\text{T}}$

利用算法1可得 $mid\left(X\right)={\left(\begin{array}{cccccc}-1.8809& 0& 0& -0.9969& 0& 1.2501\end{array}\right)}^{\text{T}},$

$mid\left(\hat{A}\right)=\left[\begin{array}{ccccc}0.6405& -0.7614& 0.0742& -1.4398& 0.5767\\ -0.5122& 0.9370& -0.9757& -0.0092& -0.0319\\ -0.2463& 0.1855& -0.2451& 0.2855& -0.4537\\ 0.1683& 0.0526& -0.4055& -0.8076& 0.3922\\ 0.4236& -0.5423& -0.9047& -0.4947& 0.2687\\ -0.8299& -0.4305& -0.2765& -1.0113& -1.0731\end{array}\right]$

$rad\left(\hat{A}\right)=3.8858e-16$， $rad\left(\hat{X}\right)=2.9976e-15$， $length\left(A-\hat{A}\right)=2.7967e-05$。

Table 1. The output results of Example 2

表1. 例2的输出结果

例2 系数矩阵为接近秩亏随机矩阵的超定线性系统 $Ax=b$，表给出了不同的 $\left(m,n,r\right)$ 取值，算法1的输出结果见表1。

NOTES

^*通讯作者。