最近亏损矩阵的可信验证

The Verification of the Nearest Defective Matrix

摘要: 给定具有n个互异特征值的阶矩阵Ã，以及其近似特征值λ，本文提出的验证算法可计算给定矩阵Ã的微小摄动区间矩阵A，以及给定实数λ的微小摄动区间_∧^Λ，使得在区间矩阵A中存在一个实亏损矩阵。

Abstract: Given an  matrix Ã having n distinct eigenvalues and its approximate eigenvalue λ, we propose a verification algorithm for constructing an interval matrix A near to Ã and an interval _∧^Λ near to λ. The computed interval matrix A is guaranteed to contain a real defective matrix .

1. 引言

如果一个 $n\times n$ 矩阵有n个不同的特征值，则称其为单矩阵。若一个 $n\times n$ 矩阵其线性无关特征向量的个数少于n，则称其为亏损矩阵。单矩阵到其附近亏损矩阵的距离与矩阵特征值的灵敏度分析密切相关。即使矩阵的所有特征值都可以彼此很好地分离，特征值的计算也可能是病态的。Kublanovskaya [1] 通过构造一个子空间的辅助基，提供了一种计算亏损矩阵特征值的方法。Ruhe [2]、Sridhar和Jordan [3]、Kågström和Ruhe [4]、Demmel [5] 分别给出了几种计算亏损矩阵约当标准形和相关不变子空间的稳定算法。Chatelin [6] 基于对牛顿法的修正，设计算法计算亏损特征值的不变子空间。Wilkinson [7] [8] [9] 深入地研究了亏损矩阵特征值问题，并给出了基于摄动分析的全体特征值的显式表达式。Wilkinson [10] [11] 详细地讨论了如何计算与给定单矩阵距离最近的亏损矩阵。Malyshev [12] 给出了一个单矩阵到具有多重特征值矩阵的距离2范数的具体公式。Lippert和Edelman [13] 基于微分几何和奇点理论，计算了一个单矩阵到一个具有二维约当块的矩阵簇的距离。利用矩阵的特征值及特征向量，Alam和Bora [14] 给出了构造最接近给定单矩阵亏损矩阵的数值算法。Alam、Bora、Byers和Overton [15] 通过伪谱分量的合并对文献 [15] 中算法进行了优化。Akinola、Freitag和Spence [16] 提出了基于隐式行列式法的最近亏损矩阵的快速算法。

Rump基于区间算法，开发了Matlab中的INTLAB工具箱，将其用于可信计算，见文献 [17] [18]。本文利用给定单矩阵的近似特征值，设计了基于F范数的单矩阵附近最近亏损矩阵的可信验证算法。

2. 预备部分

令 $ℝ$ 表示全体实数集，设A为 $n\times n$ 矩阵， $\text{tr}\left(A\right)$ 表示A的迹， $A(:)$ 表示由A的每一列合并得到的长列向量。如果A为 $m\times n$ 矩阵， $A_{i_1:i_2,:}$ 表示由矩阵A的第 $i_1$ 行到第 $i_2$ 行所构成的矩阵， $A_{:,j_1:j_2}$ 表示由矩阵A的第 $j_1$ 列到第 $j_2$ 列所构成的矩阵。令 $O_{m,n}$ 表示 $m\times n$ 阶零矩阵， $I_n$ 表示n阶单位阵。

对给定的数 $\stackrel{⌣}{x},\bar{x}$，若 $\stackrel{⌣}{x}<\bar{x}$，则称 $X=\left[\stackrel{⌣}{x},\bar{x}\right]=\left\{x\in ℝ:\stackrel{⌣}{x}\le x\le \bar{x}\right\}$ 为区间，令 $\sup \left(X\right)=\bar{x}$ 和 $\inf \left(X\right)=\stackrel{⌣}{x}$ 分别表示区间 $X$ 的右边界和左边界。设 $\mathbb{I}ℝ$ 表示全体区间的集合，分量为区间的向量和矩阵分别被称为区间向量和区间矩阵。区间向量 $X={\left(X_1,X_2,\cdots ,X_n\right)}^{\text{T}}$ 是一个n元组，其分量 $X_i\in \mathbb{I}ℝ$。定义实向量

$\sup \left(X\right)={\left(\sup \left(X_{\text{1}}\right),\sup \left(X_2\right),\cdots ,\sup \left(X_n\right)\right)}^n$ , $\inf \left(X\right)={\left(\inf \left(X_{\text{1}}\right),\inf \left(X_2\right),\cdots ,\inf \left(X_n\right)\right)}^n$

分别为区间向量 $X$ 的右边界和左边界。

对于区间矩阵 $M\in \mathbb{I}{ℝ}^{n\times n}$，若对满足条件 $M\in M$ 的任意实矩阵M都是对称矩阵，则称区间矩阵 $M$ 为区间对称矩阵。如果 $M$ 内的对称矩阵均为正定矩阵，则区间对称矩阵 $M$ 为区间对称正定矩阵。给定区间对称矩阵 $M$，如果INTLAB工具箱 [17] 中的代码输出1，则 $M$ 为区间正定对称矩阵。

定义1 (见 [18] )：给定矩阵 $A\in {ℝ}^{m\times n}$， $m\ge n$，定义矩阵A的余秩为 $corank\left(A\right)=n-rank\left(A\right)$。对于阈值 $\delta >0$，如果矩阵A的奇异值 ${\sigma }_1\left(A\right),\cdots ,{\sigma }_n\left(A\right)$ 满足

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

则我们说A有数值 $\delta$ 余秩q，记为 $coran{k}_{\delta }\left(A\right)=q$。

引理1 (见 [19] )：给定矩阵 $A\in {ℝ}^{n\times n}$，若存在矩阵 $R\in {ℝ}^{n\times n}$，使得 ${\Vert I_n-RA\Vert }_p<1$，其中 $p\in \left\{1,2,\infty \right\}$，则A是非奇异的。

定理1 (见 [19] )：对于给定的区间矩阵 $A\in \mathbb{I}{ℝ}^{n\times n}$ 和区间向量 $b\in \mathbb{I}{ℝ}^n$，若INTLAB工具箱中verifylss函数成功地输出区间向量 $X\subset \mathbb{I}{ℝ}^n$，则 $X$ 满足条件

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

定理2 (见 [19] )：设 $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 \mathbb{I}{ℝ}^n$，其中 $0\in X$，设 $f^{\prime }\left(\tilde{x}+X\right)$ 为区间 $\tilde{x}+X$ 处 $f$ 的雅可比矩阵。对于 $R\in {ℝ}^{n\times n}$ 和 $M\subset \mathbb{I}{ℝ}^{n\times n}$ 满足 $f^{\prime }\left(\tilde{x}+X\right)\subseteq M$，假设

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

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

3. 主要结论

给定一个具有n个互异特征值的 $n\times n$ 矩阵 $\tilde{A}$。若对某个 $\bar{\lambda }\in ℝ$，若 $coran{k}_{\delta }\left(\tilde{A}-\bar{\lambda }{I}_n\right)=q$ ,则本文的主要工作是计算给定矩阵 $\tilde{A}$ 的微小摄动区间矩阵 $\bar{A}$，以及给定实数 $\bar{\lambda }$ 的微小摄动区间 $\hat{\Lambda }$，使得在区间矩阵 $\bar{A}$ 中存在一个实矩阵 $\hat{A}$，在区间 $\hat{\Lambda }$ 中存在一个实数 $\hat{\lambda }$，算法保证 $\hat{\lambda }$ 为矩阵 $\hat{A}$ 几何重数为q的亏损特征值。

3.1. 数值部分

记 $\tilde{A}-\bar{\lambda }{I}_n$ 的奇异值分解为

$\begin{array}{c}\tilde{A}-\bar{\lambda }{I}_n=U\left(\bar{\lambda }\right)\Sigma \left(\bar{\lambda }\right)V{\left(\bar{\lambda }\right)}^{\text{T}}\\ =\left(u_1{u}_2\cdots u_n\right)diag\left({\sigma }_1,{\sigma }_2,\cdots \text{,}{\sigma }_n\right){\left(v_1{v}_2\cdots v_n\right)}^{\text{T}},\end{array}$

其中 ${\sigma }_1\left(A\right)\ge \cdots \ge {\sigma }_{n-q}\left(A\right)>\delta \ge {\sigma }_{n-q+1}\left(A\right)$。引入未定元矩阵

$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}\\ ⋮& ⋮& \ddots & ⋮\\ e_{n,1}& e_{n,2}& \cdots & e_{n,n}\end{array}\right)$ (1)

为了简便起见，用 $e$ 表示向量 $E(:)$。定义矩阵

$C\left(\lambda ,e\right)=\left(\begin{array}{cc}\tilde{A}-\lambda I_n+E& L\left(\bar{\lambda }\right)\\ R^{\text{T}}\left(\bar{\lambda }\right)& O_{q,q}\end{array}\right)$ (2)

其中 $L\left(\bar{\lambda }\right)=\left(u_{n-q+1}{u}_{n-q+2}\cdots u_n\right)$， $R\left(\bar{\lambda }\right)=\left(v_{n-q+1}{v}_{n-q+2}\cdots v_n\right)$。易知， $\lambda =\bar{\lambda }$， $e=0$ 时，矩阵 $C\left(\bar{\lambda },e\right)$ 非奇异。

引理2：假设 $coran{k}_{\delta }\left(\tilde{A}-\bar{\lambda }{I}_n\right)=q$， $\left|\tilde{\lambda }-\bar{\lambda }\right|<{\zeta }_1$ 且 ${\Vert \tilde{e}\Vert }_{\infty }<{\zeta }_2$。如果

${\zeta }_{\text{1}}+{\zeta }_{\text{2}}<\frac{1}{\sqrt{n+q}{\Vert C{\left(\bar{\lambda },0\right)}^{-1}\Vert }_{\infty }}$ (3)

成立，则当 $\lambda =\tilde{\lambda }$， $e=\tilde{e}$ 时，矩阵 $C\left(\lambda ,e\right)$ 非奇异。

证明：由已知可得

$\begin{array}{c}{\Vert I_{n+1}-C{\left(\bar{\lambda },0\right)}^{-1}C\left(\tilde{\lambda },\tilde{e}\right)\Vert }_2\le \sqrt{n+q}{\Vert I_{n+1}-C{\left(\bar{\lambda },0\right)}^{-1}C\left(\tilde{\lambda },\tilde{e}\right)\Vert }_{\infty }\\ \le \sqrt{n+q}{\Vert C{\left(\bar{\lambda },0\right)}^{-1}\Vert }_{\infty }{\Vert C\left(\bar{\lambda },0\right)-C\left(\tilde{\lambda },\tilde{e}\right)\Vert }_{\infty }\\ <\sqrt{n+q}{\Vert C{\left(\bar{\lambda },0\right)}^{-1}\Vert }_{\infty }\left(\left|\tilde{\lambda }-\bar{\lambda }\right|+{\Vert \tilde{e}\Vert }_{\infty }\right)\\ <1.\end{array}$

则由引理1，可知 $C\left(\tilde{\lambda },\tilde{e}\right)$ 非奇异。 □

令 $X\left(\lambda ,e\right)$ 和 $Y\left(\lambda ,e\right)$ 分别是下列线性系统

$\left(\begin{array}{cc}\tilde{A}-\lambda I_n\text{+}E& L\left(\bar{\lambda }\right)\\ R^{\text{T}}\left(\bar{\lambda }\right)& O_{q,q}\end{array}\right)\left(\begin{array}{c}X\left(\lambda ,e\right)\\ Y\left(\lambda ,e\right)\end{array}\right)=\left(\begin{array}{c}{O}_{n,q}\\ I_q\end{array}\right)$ (4)

解的前n行和后q行。由引理2可知，存在以 $\left(\bar{\lambda },0\right)$ 为中心的邻域，使得对于该邻域内任意的 $\left(\tilde{\lambda },\tilde{e}\right)$，矩阵 $C\left(\tilde{\lambda },\tilde{e}\right)$ 非奇异。进一步，在该邻域内，系统(4)解向量的每个分量关于所有变量都充分可微。对系统(4)两端关于 $\lambda$ 求导可得下列系统

$\left(\begin{array}{cc}\tilde{A}-\lambda I_n\text{+}E& L\left(\bar{\lambda }\right)\\ R^{\text{T}}\left(\bar{\lambda }\right)& O_{q,q}\end{array}\right)\left(\begin{array}{c}\frac{{\partial }^{\left(l\right)}X\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l\right)}}\\ \frac{{\partial }^{\left(l\right)}Y\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l\right)}}\end{array}\right)=\left(\begin{array}{c}l\frac{{\partial }^{\left(l-1\right)}X\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l-1\right)}}\\ O_{q,q}\end{array}\right),l=1,2,\cdots$ (5)

对于 $\left(i,j\right)$，其中 $1\le i\le n,\text{}1\le j\le n$，对系统(4)关于变量 $e_{i,j}$ 求偏导，可得

$\left(\begin{array}{cc}\tilde{A}+E-\lambda I_n& L\left(\bar{\lambda }\right)\\ R^{\text{T}}\left(\bar{\lambda }\right)& O_{q,q}\end{array}\right)\left(\begin{array}{c}\frac{\partial X\left(\lambda ,e\right)}{\partial e_{i,j}}\\ \frac{\partial Y\left(\lambda ,e\right)}{\partial e_{i,j}}\end{array}\right)=-\left(\begin{array}{c}{O}_{i-1,q}\\ X{\left(\lambda ,e\right)}_{j,:}\\ O_{n+q-i,q}\end{array}\right).$ (6)

对于 $\left(i_{\text{1}},j_{\text{1}}\right),\left(i_2,j_2\right)$，其中 $1\le i_{\text{1}}\le n,\text{}1\le j_{\text{1}}\le n,1\le i_2\le n,\text{}1\le j_2\le n$，对系统(6)关于变量 $e_{i_{\text{1}},j_{\text{1}}},e_{i_2,j_2}$ 求偏导，有

$\left(\begin{array}{cc}\tilde{A}+E-\lambda I_n& L\left(\bar{\lambda }\right)\\ R^{\text{T}}\left(\bar{\lambda }\right)& O_{q,q}\end{array}\right)\left(\begin{array}{c}\frac{{\partial }^{\left(\text{2}\right)}X\left(\lambda ,e\right)}{\partial e_{i_{\text{1}},j_{\text{1}}}\partial e_{i_{\text{2}},j_{\text{2}}}}\\ \frac{{\partial }^{\left(\text{2}\right)}Y\left(\lambda ,e\right)}{\partial e_{i_{\text{1}},j_{\text{1}}}\partial e_{i_{\text{2}},j_{\text{2}}}}\end{array}\right)=-\left(\begin{array}{c}{O}_{i_{\text{1}}-1,q}\\ {\left(\frac{\partial X\left(\lambda ,e\right)}{\partial e_{i_{\text{2}},j_{\text{2}}}}\right)}_{j_{\text{1}},:}\\ O_{n+q-i_{\text{1}},q}\end{array}\right)-\left(\begin{array}{c}{O}_{i_{\text{2}}-1,q}\\ {\left(\frac{\partial X\left(\lambda ,e\right)}{\partial e_{i_{\text{1}},j_{\text{1}}}}\right)}_{j_{\text{2}},:}\\ O_{n+q-i_{\text{2}},q}\end{array}\right).$ (7)

对于 $\left(l,i,j\right)$，其中 $l\ge 1,1\le i\le n,\text{}1\le j\le n$，系统(5)两端对 $e_{i,j}$ 求偏导，可得

$\left(\begin{array}{cc}\tilde{A}-\lambda I_n\text{+}E& L\left(\bar{\lambda }\right)\\ R^{\text{T}}\left(\bar{\lambda }\right)& O_{q,q}\end{array}\right)\left(\begin{array}{c}\frac{{\partial }^{\left(l\text{+1}\right)}X\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l\right)}{e}_{i,j}}\\ \frac{{\partial }^{\left(l\text{+1}\right)}Y\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l\right)}{e}_{i,j}}\end{array}\right)=\left(\begin{array}{c}l\frac{{\partial }^{\left(l\right)}X\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l-1\right)}{e}_{i,j}}\\ O_{q,q}\end{array}\right)-\left(\begin{array}{c}{O}_{i-1,q}\\ {\left(\frac{{\partial }^{\left(l\right)}X\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l\right)}}\right)}_{j,:}\\ O_{n+q-i,q}\end{array}\right).$ (8)

对于 $\left(l,i_{\text{1}},j_{\text{1}}\right),\left(l,i_2,j_2\right)$，系统(5)关于变量 $e_{i_{\text{1}},j_{\text{1}}}$ 和 $e_{i_2,j_2}$ 求偏导，可得

$\begin{array}{l}\left(\begin{array}{cc}\tilde{A}+E-\lambda I_n& L\left(\bar{\lambda }\right)\\ R^{\text{T}}\left(\bar{\lambda }\right)& O_{q,q}\end{array}\right)\left(\begin{array}{c}\frac{{\partial }^{\left(l+2\right)}X\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l\right)}\partial e_{i_{\text{1}},j_{\text{1}}}\partial e_{i_{\text{2}},j_{\text{2}}}}\\ \frac{{\partial }^{\left(l+2\right)}Y\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l\right)}\partial e_{i_{\text{1}},j_{\text{1}}}\partial e_{i_{\text{2}},j_{\text{2}}}}\end{array}\right)\\ =\left(\begin{array}{c}l\frac{{\partial }^{\left(l+1\right)}X\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l-1\right)}\partial e_{i_{\text{1}},j_{\text{1}}}\partial e_{i_{\text{2}},j_{\text{2}}}}\\ O_{q,q}\end{array}\right)-\left(\begin{array}{c}{O}_{i_{\text{1}}-1,q}\\ {\left(\frac{{\partial }^{\left(l+1\right)}X\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l\right)}\partial e_{i_{\text{2}},j_{\text{2}}}}\right)}_{j_{\text{1}},:}\\ O_{n+q-i_{\text{1}},q}\end{array}\right)-\left(\begin{array}{c}{O}_{i_{\text{2}}-1,q}\\ {\left(\frac{{\partial }^{\left(l+1\right)}X\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l\right)}\partial e_{i_{\text{1}},j_{\text{1}}}}\right)}_{j_{\text{2}},:}\\ O_{n+q-i_{\text{2}},q}\end{array}\right).\end{array}$ (9)

注1：对于满足条件 $coran{k}_{\delta }\left(\tilde{A}-\bar{\lambda }{I}_n\right)=q$ 的矩阵 $\tilde{A}$ 及实数 $\bar{\lambda }$，数值部分的主要工作是求解以下约束优化问题

$\begin{array}{l}\min \frac{1}{2}{\Vert e\Vert }_2^2\\ \text{s}\text{.t}\text{.}Y\left(\lambda ,e\right)=O_{q,q},\end{array}$ (10)

$\begin{array}{l}\frac{{\partial }^{\left(l\right)}Y\left(\lambda ,e\right)}{\partial {\lambda }^{\left(l-1\right)}}=O_{q,q},1\le l\le k-1,\\ \frac{{\partial }^{\left(k\right)}Y\left(\lambda ,e\right)}{\partial {\lambda }^{\left(k\right)}}\ne O_{q,q}.\end{array}$ (11)

利用牛顿迭代法来计算优化问题(10)的数值稳定点。需要注意的是，条件(11)中出现的正整数k是根据阈值 $\delta$ 计算出来的。初始k值应满足以下条件

${\Vert \frac{{\partial }^{\left(k-1\right)}Y\left(\bar{\lambda },0\right)}{\partial {\lambda }^{\left(k-1\right)}}\Vert }_{\infty }<\delta ,{\Vert \frac{{\partial }^{\left(k\right)}Y\left(\bar{\lambda },0\right)}{\partial {\lambda }^{\left(k\right)}}\Vert }_{\infty }\ge \delta .$

对于更新的 $\left(\tilde{\lambda },\tilde{e}\right)$，如果

${\Vert \frac{{\partial }^{\left(k-1\right)}Y\left(\bar{\lambda },0\right)}{\partial {\lambda }^{\left(k-1\right)}}\Vert }_{\infty }<\delta$

不成立，则需要更新k，使得条件

${\Vert \frac{{\partial }^{\left(k-1\right)}Y\left(\bar{\lambda },0\right)}{\partial {\lambda }^{\left(k-1\right)}}\Vert }_{\infty }<\delta ,{\Vert \frac{{\partial }^{\left(k\right)}Y\left(\bar{\lambda },0\right)}{\partial {\lambda }^{\left(k\right)}}\Vert }_{\infty }\ge \delta .$

成立。 □

注2：令 $Hess\left(\frac{{\partial }^{\left(i\right)}Y\left(\lambda ,e\right)}{\partial {\lambda }^{\left(i\right)}}\right)$ 为 $\frac{{\partial }^{\left(i\right)}Y\left(\lambda ,e\right)}{\partial {\lambda }^{\left(i\right)}}$ 的Hession矩阵，其中 $0\le i\le k-1$，定义非线性系统

$F\left(\lambda ,e\right)=\left(\begin{array}{c}Y\left(\lambda ,e\right)(:)\\ \frac{\partial Y\left(\lambda ,e\right)}{\partial \lambda }(:)\\ ⋮\\ \frac{{\partial }^{\left(k-1\right)}Y\left(\lambda ,e\right)}{\partial {\lambda }^{\left(k-1\right)}}(:)\end{array}\right)=0,$ (12)

则优化问题(10)相应的拉格朗日函数为

$L\left(\lambda ,e,\omega \right)=\frac{1}{2}{\Vert e\Vert }_2^2+{\omega }^{\text{T}}F\left(\lambda ,e\right),$ (13)

其中 $\omega ={\left({\omega }_0,{\omega }_1,\cdots ,{\omega }_{k{q}^2-1}\right)}^{\text{T}}$ 为拉格朗日乘子，函数 $L\left(\lambda ,e,\omega \right)$ 的梯度和Hession矩阵分别为

$L^{\prime }\left(\lambda ,e,\omega \right)=\left(\begin{array}{c}0\\ e\\ O_{k{q}^2,1}\end{array}\right)+\left(\begin{array}{c}{F}^{\prime }{\left(\lambda ,e\right)}^{\text{T}}\omega \\ F\left(\lambda ,e\right)\end{array}\right),$ (14)

$L^{″}\left(\lambda ,e,\omega \right)=\left(\begin{array}{ccc}0& O_{1,n^2}& O_{1,k{q}^2}\\ O_{n^2,1}& I_{n^2}& O_{n^2,k{q}^2}\\ O_{k{q}^2,1}& O_{k{q}^2,n^2}& O_{k{q}^2,k{q}^2}\end{array}\right)+\left(\begin{array}{cc}\text{Hess}\left(F\left(\lambda ,e\right)\right)& F^{\prime }{\left(\lambda ,e\right)}^{\text{T}}\\ F^{\prime }\left(\lambda ,e\right)& O_{k{q}^2,k{q}^2}\end{array}\right).$ (15)

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

${\left(\begin{array}{ccc}\tilde{\lambda }& \tilde{e}& \tilde{\omega }\end{array}\right)}^{\text{T}}\leftarrow {\left(\begin{array}{ccc}\tilde{\lambda }& \tilde{e}& \tilde{\omega }\end{array}\right)}^{\text{T}}-L^{″}{\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)}^{-1}{L}^{\prime }\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)$ (16)

来更新 $\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)$，使得 $L^{\prime }\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)$ 接近零向量。 □

3.2. 验证部分

引理3：若向量 $\left(\hat{\lambda },\hat{e}\right)$ 满足条件 $Y\left(\hat{\lambda },\hat{e}\right)=O_{q,q}$，且

$\left|\hat{\lambda }-\bar{\lambda }\right|+{\Vert \hat{e}\Vert }_{\infty }<\frac{1}{\sqrt{n+q}{\Vert C{\left(\bar{\lambda },0\right)}^{-1}\Vert }_{\infty }}$

成立，则

$corank\left(\tilde{A}+E-\hat{\lambda }{I}_n\right)=q.$ (17)

证明：根据引理2可知， $C\left(\hat{\lambda },\hat{e}\right)$ 非奇异。若 $Y\left(\hat{\lambda },\hat{e}\right)=O_{q,q}$，则有

$\left(\tilde{A}+\hat{E}-\hat{\lambda }{I}_n\right)X\left(\hat{\lambda },\hat{e}\right)=O_{n,q},R^{\text{T}}\left(\hat{\lambda }\right)X\left(\hat{\lambda },\hat{e}\right)=I_q.$

显然， $X\left(\hat{\lambda },\hat{e}\right)$ 是线性无关的。假设 $corank\left(\tilde{A}+\hat{E}-\hat{\lambda }{I}_n\right)>q$，对于某非零向量 $\beta \in \ker \left(\tilde{A}+\hat{E}-\hat{\lambda }{I}_n\right)$ 有

$\left(W\left(\hat{\lambda },\hat{e}\right)\beta \right)$ 满秩，则对于非零实向量 $\gamma$，有

$\left(\begin{array}{c}\tilde{A}+\hat{E}-\hat{\lambda }{I}_n\\ R^T\left(\hat{\lambda }\right)\end{array}\right)\left(W\left(\hat{\lambda },\hat{e}\right)\beta \right)\gamma =0\text{.}$

上式与假设条件矛盾，故(17)成立。 □

假设1：雅可比矩阵 $F^{\prime }\left(0\right)$ 满秩。

命题1：假设 $\left(\hat{\lambda },\hat{e}\right)$ 是系统(12)的解，且满足条件

$\left|\hat{\lambda }-\bar{\lambda }\right|+{\Vert \hat{e}\Vert }_{\infty }<\frac{1}{\sqrt{n+q}{\Vert C{\left(\bar{\lambda },0\right)}^{-1}\Vert }_{\infty }},$

则对每一个正整数l，其中 $\text{1}\le l\le k$，下列向量的序列

$\frac{{\partial }^{\left(t\right)}X{\left(\hat{\lambda },\hat{e}\right)}_{:,s}}{\partial {\lambda }^{\left(t\right)}}1\le s\le q,0\le t\le l-1$ (18)

构成 $\ker {\left(\tilde{A}+\hat{E}-\hat{\lambda }{I}_n\right)}^l$ 的基。

证明：设实矩阵 $\hat{E}\in {ℝ}^{n\times n}$ 满足条件 $\hat{E}(:)=\hat{e}$。当 $l=1$ 时，根据引理3知， $W\left(\hat{\lambda },\hat{e}\right)$ 是 $\ker \left(\tilde{A}+\hat{E}-\hat{\lambda }{I}_n\right)$ 的一个基。假设正整数 $\text{1}\le l\le k-1$，向量组(18)是 $\ker {\left(\tilde{A}+\hat{E}-\hat{\lambda }{I}_n\right)}^l$ 的一组基。下面考虑 $l+\text{1}$ 时的情况。显然

$\ker \left({\left(\tilde{A}+\hat{E}-\hat{\lambda }{I}_n\right)}^{l\text{+1}}\right)\supseteq \text{span}\left\{\frac{{\partial }^{\left(t\right)}W{\left(\hat{\lambda },\hat{e}\right)}_{:,s}}{\partial {\lambda }^{\left(t\right)}}1\le s\le q,0\le t\le l\right\}.$ (19)

接下来证明向量组

$\frac{{\partial }^{\left(t\right)}X{\left(\hat{\lambda },\hat{e}\right)}_{:,s}}{\partial {\lambda }^{\left(t\right)}}1\le s\le q,0\le t\le l$ (20)

线性无关。然而，如果存在一个q维非零实向量序列 ${\alpha }^{\left(\text{0}\right)},{\alpha }^{\left(\text{1}\right)}，\cdots ,{\alpha }^{\left(l\right)}$，使得

${\displaystyle \underset{t=0}{\overset{l}{\sum }}\frac{{\partial }^{\left(t\right)}W\left(\hat{\lambda },\hat{e}\right)}{\partial {\lambda }^{\left(t\right)}}{\alpha }^{\left(t\right)}=0,}$

则

${\displaystyle \underset{t=1}{\overset{l}{\sum }}\frac{{\partial }^{\left(t-1\right)}W\left(\hat{\lambda },\hat{e}\right)}{\partial {\lambda }^{\left(t-1\right)}}{\alpha }^{\left(t\right)}=0,}$

由数学归纳法可知向量组 ${\alpha }^{\left(t\right)},t=1,\cdots ,l$ 全是零向量。因此 $W\left(\hat{\lambda },\hat{e}\right){\alpha }^{\left(0\right)}=0$，而 $W\left(\hat{\lambda },\hat{e}\right)$ 列向量线性无关，故 ${\alpha }^{\left(0\right)}=0$ 成立。综上可得矛盾，所以(20)线性无关。

最后假设存在非零实向量 $z$ 满足条件

$z\in \ker {\left(\tilde{A}+\hat{E}-\hat{\lambda }{I}_n\right)}^{l\text{+1}}\backslash \text{span}\left\{\frac{{\partial }^{\left(t\right)}W{\left(\hat{\lambda },\hat{e}\right)}_{:,s}}{\partial {\lambda }^{\left(t\right)}}1\le s\le q,0\le t\le l\right\},$ (21)

则存在不全为零的向量组 ${\beta }^{\left(t\right)}\in {ℝ}^q$， $t=1,\cdots ,l$，使得

$\left(\tilde{A}+\hat{E}-\hat{\lambda }{I}_n\right)z={\displaystyle \underset{t=0}{\overset{l-1}{\sum }}\frac{{\partial }^{\left(t\right)}W\left(\hat{\lambda },\hat{e}\right)}{\partial {\lambda }^{\left(t\right)}}{\beta }^{\left(t+1\right)}}$

成立。因此，存在一个q维向量 ${\beta }^{\left(\text{0}\right)}$，使得

$z=W\left(\hat{\lambda },\hat{e}\right){\beta }^{\left(0\right)}+{\displaystyle \underset{t=0}{\overset{l-1}{\sum }}\frac{\text{1}}{t+1}\frac{{\partial }^{\left(t+1\right)}W\left(\hat{\lambda },\hat{e}\right)}{\partial {\lambda }^{\left(t+1\right)}}{\beta }^{\left(t+1\right)}},$

这与(21)式矛盾。 □

注4：利用命题1和verifylss函数可知，验证算法计算包含实向量 $\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)$ 的区间向量 $\left(\hat{\Lambda },\hat{E},\hat{W}\right)$，使得区间向量 $\left(\hat{\Lambda },\hat{E},\hat{W}\right)$ 中包含实向量 $\left(\hat{\lambda },\hat{e},\hat{\omega }\right)$，其满足条件 $L^{\prime }\left(\hat{\lambda },\hat{e},\hat{\omega }\right)=0$。对于包含实向量 $\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)$ 的任意区间向量 $\left(\tilde{\Lambda },\tilde{E},\tilde{W}\right)$，利用定理1和verifylss函数，求解当 $\lambda =\tilde{\Lambda },e=\widetilde{{\rm E}},\omega =\tilde{W}$ 时的线性系统(6) (7) (8) (9)，可得到区间梯度向量 $L^{\prime }\left(\tilde{\Lambda },\tilde{E},\tilde{W}\right)$ 和区间雅可比矩阵 $L^{″}\left(\tilde{\Lambda },\tilde{E},\tilde{W}\right)$。 □

4. 主要算法

本节提出一种计算及验证优化问题(10)局部最优解的算法。对于实数x，代码 $intval\left(x\right)$ 输出区间 $\left[x,x\right]$，且对于实数 $a,b$ 满足 $a<b$，代码 $infsub\left(a,b\right)$ 输出区间 $\left[a,b\right]$。对于区间 $X$，代码 $hull\left(X,0\right)$ 返回同时包含 $X$ 和 $0$ 的最小区间。

算法1

输入 $\tilde{A}$ ： $n\times n$ 对称矩阵；

$\delta$ ：数值秩的容差；

$tol$ ：数值牛顿迭代的容差；

$\bar{\lambda }$ ：初始近似特征值；

N：数值迭代的最大次数。

输出 $\hat{\Lambda }\in ℝ$， $\hat{E}\subset \mathbb{I}{ℝ}^{n\times n}$， $q\in \mathbb{N}$ 和 $k\in \mathbb{N}$ 或“失败”。

步骤1 对 $\tilde{A}-\bar{\lambda }{I}_n$ 进行奇异值分解。若 $corank\left(\tilde{A}+E-\hat{\lambda }{I}_n\right)=q$，进行第二步，否则返回失败并停止。

步骤2 数值部分。

步骤2.1 初始化 $iter=0,k=0,\tilde{\lambda }=\bar{\lambda },\tilde{e}=0$。

步骤2.2 当 ${\Vert \frac{{\partial }^{\left(k\right)}Y\left(\tilde{\lambda },\tilde{e}\right)}{\partial {\lambda }^{\left(k\right)}}\Vert }_{\infty }<\delta$ 进行

令 $k\leftarrow k+1$ ；

当 $\lambda =\tilde{\lambda },e=\tilde{e},l=k$，求解线性系统(5)得到 $\frac{{\partial }^{\left(k\right)}Y\left(\tilde{\lambda },\tilde{e}\right)}{\partial {\lambda }^{\left(k\right)}}$。

步骤2.3 如果 $k=1$，返回失败并停止，否则进行下一步。

步骤2.4 求解线性系统(5) (6) (8)得到 $F^{\prime }\left(\lambda ,e\right)$。

步骤2.5 初始化

$\tilde{\omega }\leftarrow \left(F^{\prime }{\left(\tilde{\lambda },\tilde{e}\right)}^{\text{T}}\right)+\left(\begin{array}{c}\tilde{\lambda }\\ \tilde{e}\end{array}\right)$。

步骤2.6 当 $\Vert L^{\prime }\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)\Vert <tol$ 且 $iter<N$，进行

$\lambda =\tilde{\lambda },e=\tilde{e}$，求解线性系统(7) (9)得到 $L^{″}\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)$ ；

利用(16)更新 $\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)$ ；

更新 $iter\leftarrow iter+1$ ；

如果 ${\Vert \frac{{\partial }^{\left(k\right)}Y\left(\tilde{\lambda },\tilde{e}\right)}{\partial {\lambda }^{\left(k\right)}}\Vert }_{\infty }<\delta$，则返回步骤2.2；

求解线性系统(6) (8)计算 $L^{\prime }\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)$。

步骤2.7 如果 $iter=N$，返回失败并停止。

步骤3 验证部分。

步骤3.1 初始化 $R=L^{″}{\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)}^{-1},Z=intval\left(L^{″}{\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)}^{-1}{L}^{\prime }\left(\tilde{\lambda },\tilde{e},\tilde{\omega }\right)\right),X=Z$ 及 $iter=0$。

步骤3.2 当 $iter\le 10$，进行

令 $iter\leftarrow iter+1$ ；

令 $Y=hull\left(Xinfsub\left(0.9,1.1\right)+{10}^{-20}infsup\left(-1,1\right),0\right)$ ；

计算 $L^{″}\left(\tilde{\lambda }+Y_1,\tilde{e}+Y_{2:n^2+1},\tilde{\omega }+Y_{n^2+2:n^2+k{q}^2+1}\right)$ ；

令 $X=Z+\left(I_{n^2+k{q}^2+1}-RM\right)Y$ ；

若 $X\subseteq \mathrm{int}\left(Y\right)$，则令

${\left(\begin{array}{ccc}\hat{\Lambda }& \hat{E}& \hat{W}\end{array}\right)}^{\text{T}}={\left(\begin{array}{ccc}\tilde{\lambda }& \tilde{e}& \tilde{\omega }\end{array}\right)}^{\text{T}}+X$。

步骤3.3 如果 $iter=10$ 或者 $O\notin \frac{{\partial }^{\left(k\right)}\left(\hat{\Lambda },\hat{E}\right)}{\partial {\lambda }^{\left(k\right)}}$，则输出 $\hat{\Lambda },\hat{E},k$ 和q并停止。

定理3：给定矩阵 $\tilde{A}$ 及其近似特征值 $\bar{\lambda }$，如果算法1输出 $\hat{\Lambda },\hat{E},k$ 和q，则存在 $\hat{\lambda }\in \hat{\Lambda }$， $\hat{e}\in \hat{E}$ 使得 $\left(\hat{\lambda },\hat{e}\right)$ 是优化问题(10)的局部最优解。进一步， $\hat{\lambda }$ 是矩阵 $\tilde{A}+\hat{E}$ 几何重数为q的亏损特征值。

证明：由定理2可知，存在 $\left(\hat{\lambda },\hat{e},\hat{\omega }\right)\in \left(\hat{\Lambda },\hat{E},\hat{W}\right)$，使得 $L^{\prime }\left(\hat{\lambda },\hat{e},\hat{\omega }\right)=0$。由引理3和命题1可知，当 $k\ge 2$， $\hat{\lambda }$ 是矩阵 $\tilde{A}+\hat{E}$ 几何重数q的亏损特征值。 □

5. 应用实例

例1给定矩阵

$A=\left[\begin{array}{cccccccc}3.006& 2& 1.005& -1.001& -0.002& -0.001& -0.001& -1\\ 5& 2& 5& -1& -2& -1& -1& 0\\ -5.006& -3& -3.005& 2.001& 3.002& 2.001& 0.001& 2\\ -6& -1& -6& 3& 5& 3& 0& 1\\ -5& -1& -5& 1& 6& 3& 0& 0\\ 1& 0& 1& 0& -1& 1& 0& 0\\ -4& -2& -4& 1& 3& 2& 2& 2\\ 5& 0& 5& -1& -2& -1& -1& 2\end{array}\right]，$

令 $\tilde{A}$ 为

$\tilde{A}=A+{10}^{-n}\left[\begin{array}{cccccccc}0.6294& 0.9150& -0.1565& 0.3575& -0.4462& -0.1225& 0.4187& 0.9195\\ 0.8116& 0.9298& 0.8315& 0.5155& -0.9077& -0.2369& 0.5094& -0.3192\\ -0.7460& -0.6848& 0.5844& 0.4863& -0.8057& 0.5310& -0.4479& 0.1705\\ 0.8268& 0.9412& 0.9190& -0.2155& 0.6469& 0.5904& 0.3594& -0.5524\\ 0.2647& 0.9143& 0.3115& 0.3110& 0.3897& -0.6263& 0.3102& 0.5025\\ -0.8049& -0.0292& -0.9286& -0.6576& -0.3658& -0.0205& -0.6748& -0.4898\\ -0.4430& 0.6006& 0.6983& 0.4121& 0.9004& -0.1088& -0.7620& 0.0119\\ 0.0938& -0.7162& 0.8680& -0.9363& -0.9311& 0.2926& -0.0033& 0.3982\end{array}\right].$

对于不同的n，通过Matlab中的eig代码计算可知 $\tilde{A}$ 为单矩阵。对于不同的n和特征值 $\bar{\lambda }=2.0$ 的矩阵 $\tilde{A}$，表1给出算法1的结果。

Table 1. The performance of Algorithm 1 for Example 1

表1. 例1中算法1的计算结果

例2设 $n\times n$ 矩阵J有一个 $m_1\times m_1$ 的约当块和一个 $\left(m_2-m_1\right)\times \left(m_2-m_1\right)$ 的约当块，两个约当块都与特征值2相关，令P为随机的非奇异矩阵。令 $A=PJ{P}^{-1}$，令

$\begin{array}{l}\tilde{A}=\\ A+{10}^{-4}\left[\begin{array}{cccccccccc}0.3840& -0.7455& -0.0570& 0.9810& -0.1911& 0.7924& 0.2354& -0.9345& 0.7623& 0.0047\\ -0.2393& 0.0260& 0.1411& 0.1012& -0.9535& 0.6953& -0.0520& -0.3028& 0.3942& 0.2712\\ 0.2910& -0.3717& 0.5989& -0.8594& -0.1535& -0.2398& 0.8116& 0.4159& -0.9761& 0.5773\\ -0.8693& -0.1310& 0.1309& -0.8017& -0.6438& -0.4670& 0.3047& -0.6064& 0.1972& 0.7100\\ -0.4320& -0.5132& -0.8087& 0.3289& 0.4160& -0.8258& 0.9311& -0.0740& 0.9080& 0.9349\\ 0.2940& 0.5696& -0.0937& 0.1187& -0.4758& -0.1019& 0.7441& -0.3800& -0.0893& -0.9620\\ -0.2757& 0.7320& -0.6840& 0.7244& 0.0309& -0.9836& -0.9540& 0.5586& 0.7507& -0.2743\\ 0.8089& -0.7438& -0.3218& 0.2609& -0.6594& 0.9109& -0.3479& -0.4232& -0.8774& -0.5354\\ 0.9181& 0.5407& -0.9680& 0.1162& -0.0441& 0.0511& 0.1022& 0.0488& 0.4116& -0.0816\\ 0.1609& -0.5070& 0.5487& 0.4986& -0.1408& 0.9865& -0.8304& -0.6444& 0.7302& -0.6024\end{array}\right]\end{array}$

矩阵A为亏损矩阵，但利用Matlab计算可知矩阵A经过微小摄动后所得的 $\tilde{A}$ 为单矩阵。对于给定的 $\bar{\lambda }=2.0$ 和 $\tilde{A}$，表2给出算法1的计算结果。

Table 2. The performance of Algorithm 1 for Example 2

表2. 例2中算法1的计算结果

基金项目

吉林省自然科学基金(批准号：20180101345JC)。

NOTES

^*通讯作者。