实对称矩阵谱的可信计算

The Verification of the Spectra of the Real Symmetric Matrix

摘要: 本文主要研究实对称矩阵谱的可信计算。给定实对称矩阵，首先利用Matlab中的eig命令求其数值谱，然后利用Kantorovich定理，设计算法计算数值谱的可信误差界。算法保证在该误差界范围内存在一实对称矩阵，该实对称矩阵的精确谱为给定实对称矩阵的数值谱。

Abstract: This paper mainly investigates the verification of the spectra of the real symmetric matrix. Given a real symmetric matrix, we firstly use eig code in Matlab to obtain its numerical spectra. Then by Kantorovich theorem, we provide an algorithm to compute verified error bound such that there exists a perturbed real symmetric matrix within computed error bound, whose exact spectra is the computed numerical spectra of the given matrix.

1. 引言

矩阵的谱为矩阵所有特征值的集合，作为分析矩阵非正规性和特征值摄动的一种有效工具，在化学、控制理论、激光、水动力稳定性等诸多领域都有着重要的应用 [1] [2] [3]。研究表明，在工程计算和实际应用中有许多问题最终都归结为矩阵特征值的计算问题，而且不同的应用会导出一些具有特殊代数结构的矩阵，矩阵代数结构的保留有助于提高特征值计算的精度和效率 [4] [5]。因此，对于具有特殊代数结构的矩阵谱的计算，在实际问题的解决过程中十分重要。

多年来，对于具有特殊代数结构的矩阵特征值求解问题的研究始终吸引着国内外的专家学者。Karow等人 [6] 研究了保持结构化摄动对线性和非线性结构矩阵特征值问题的影响。Graillat [7] [8] 分析了矩阵特征值问题的敏感性，给出了Toeplitz、Hankel、对称、Hermitian和斜Hermitan等结构矩阵的结构化条件数的定义，并证明了其非结构化条件数等于结构化条件数。Adhikari和Alam [9] 对Hamiltonian、反对称、Hermitian等矩阵束的近似特征值的结构化后向误差和结构矩阵束的谱进行分析。基于区间计算，Rump [10] 研究了结构矩阵在结构化摄动下其特征值求解的灵敏度，并证明了Toeplitz矩阵、对称矩阵、Hermitian矩阵等其他类型的结构矩阵，在依范数摄动下，结构化条件数等于非结构化条件数。Noschese和Reichel [11] 提出了一种利用秩-1或投影秩-1摄动来计算结构矩阵近似谱的新方法。

Alon等人 [12] 研究了对角项和上对角线项为独立实随机变量的随机对称矩阵的最大特征值的集中性。Edwards和Jones [13] 利用具有相同均值和方差的高斯概率密度函数，提出了一种直接分析大型对称矩阵谱的方法。给定一个对称矩阵，其每一项依赖于一个参数，HiriartUrruty和Ye [14] 研究了所有特征值的一阶灵敏度。Hladik等人 [15] 考虑了具有扰动区间项的对称矩阵的特征值问题。Hernandez等人 [16] 提出了一种贪婪算法，利用特征向量的局部化特征来计算大型稀疏对称矩阵的特征对。Reid [17] 给出了对称矩阵和广对称矩阵的一些有用的特征值和特征向量的性质。

本文利用Rump区间方法和Kantorovich定理来计算给定实对称矩阵谱的可信误差界，使得在计算的误差范围内存在一个摄动的对称矩阵，其精确谱为给定对称矩阵的数值谱。准确地说，我们将实对称矩阵谱的验证转化为非线性系统根的验证。我们利用Rump区间方法计算非线性系统在Kantorovich定理中出现的常量，然后利用Kantorovich定理计算经验证的零向量误差界作为近似解。

2. 预备知识

令 $ℝ$ 表示实数集合。设 $A\in {ℝ}^{m\times n}$， $A(:)$ 表示由矩阵A的每一列合并得到的长列向量， $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阶单位阵。对于一个 $m\times n$ 的矩阵A，令 $null\left(A\right)$ 表示矩阵A的零空间。

定义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$。

令 $\mathbb{I}ℝ$ 表示全体实区间集合。分量为区间的向量和矩阵分别被称为区间向量和区间矩阵。对于区间矩阵 $A\in \mathbb{I}{ℝ}^{n\times n}$，若对满足条件 $A\in A$ 的任意实矩阵A，A都是非奇异矩阵，则称区间矩阵 $A$ 非奇异。Rump [19] 在Matlab中为区间运算开发了INTLAB工具箱。对于一个非线性系统，若系统的Jacobian矩阵在某个区域上是Lipschitz连续的，则Kantorovich [20] 建立了Kantorovich定理，该定理给出了根据某个区域上的初始近似值信息判断牛顿迭代法是否收敛的充分条件。

引理1 (见 [21] )：给定矩阵 $A,R\in {ℝ}^{n\times n}$，如果矩阵 $I-RA$ 的谱半径小于1，则A是非奇异的。

定理1 (见 [22] )：对于给定的区间矩阵 $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 (见 [20] )：设 $f:{ℝ}^n\to {ℝ}^n$ 并且 $f={\left(f_1,f_2,\cdots ,f_n\right)}^{\text{T}}$，其中 $f_1,\cdots ,f_n$ 为连续可微函数。令 $f^{\prime }\left(x\right)$ 表示系统 $f\left(x\right)=0$ 的雅克比矩阵， $\tilde{x}\in R^n$ 为 $f\left(x\right)=0$ 的近似解，且 $f^{\prime }\left(\tilde{x}\right)$ 可逆。设B为满足条件

$\Vert f^{\prime }{\left(\tilde{x}\right)}^{-1}\Vert \le B$

的常量，令 $\kappa$ 为满足

$\Vert f^{\prime }\left(u\right)-f^{\prime }\left(v\right)\Vert \le \kappa \Vert u-v\Vert ,u,v\in \Omega$,

的Lipschitz常量，其中 $\Omega$ 是包含 $\tilde{x}$ 的一个足够大区域， $\eta$ 为满足条件

$\Vert f^{\prime }{\left(\tilde{x}\right)}^{-1}f\left(\tilde{x}\right)\Vert \le \eta$

的常量。如果对于

$\rho =\frac{1-\sqrt{1-2h}}{h}\eta$ (1)

有 $h=B\kappa \eta \le 1/2$， $\bar{U}\left(\tilde{x},\rho \right)=\left\{x:\Vert x-\tilde{x}\Vert \le \rho \right\}\subset \Omega$，那么存在解 $\hat{x}\in \bar{U}\left(\tilde{x},\rho \right)$，使得 $f\left(\hat{x}\right)=0$。

注1：正如Rall [23] 指出的，定理2中的区域 $\Omega$ 可以取为 $\bar{U}\left(\tilde{x},2\eta \right)$。如果 $\kappa$ 是该区域的Lipschitz常量，则 $h\le 1/2$，当且仅当 $\bar{U}\left(\tilde{x},\rho \right)\subset \bar{U}\left(\tilde{x},2\eta \right)$。

3. 主要结论

对于一个方阵A，如果 $A^{\text{T}}=A$，则称A为对称矩阵，实对称矩阵特征值均为实数。

令 $A^{sym}$ 表示 $n\times n$ 实对称矩阵， $E^{sym}$ 表示相应的摄动矩阵，即

$E^{sym}=\left(\begin{array}{ccccc}{\epsilon }_{1,1}& {\epsilon }_{1,2}& {\epsilon }_{1,3}& \cdots & {\epsilon }_{1,n}\\ {\epsilon }_{1,2}& {\epsilon }_{2,1}& {\epsilon }_{2,2}& \cdots & {\epsilon }_{2,n-1}\\ {\epsilon }_{1,3}& {\epsilon }_{2,2}& {\epsilon }_{3,1}& \cdots & {\epsilon }_{3,n-2}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\epsilon }_{1,n}& {\epsilon }_{2,n-1}& {\epsilon }_{3,n-2}& \cdots & {\epsilon }_{n,1}\end{array}\right)$. (2)

令 $\left\{{\tilde{\lambda }}_1,{\tilde{\lambda }}_2,\cdots ,{\tilde{\lambda }}_k\right\}$ 表示由Matlab中的eig命令计算所得矩阵 $A^{sym}$ 的所有互异特征值的全体。对于 $s=1,2,\cdots ,k$，对 $A^{sym}-{\tilde{\lambda }}_s{I}_n$ 做奇异值分解，有

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

假设对于 $s=1,2,\cdots ,k$， $coran{k}_{\delta }\left(A^{sym}-{\tilde{\lambda }}_s{I}_n\right)=q_s$，其中 $\delta$ 是一个接近0的正实数。

如果 ${\hat{\lambda }}_1,{\hat{\lambda }}_2,\cdots ,{\hat{\lambda }}_k$ 是矩阵 $A^{sym}$ 的所有精确特征值，则对于 $s=1,2,\cdots ,k$，向量组

$v_1\left({\hat{\lambda }}_s\right),v_2\left({\hat{\lambda }}_s\right),\cdots ,v_{n-q_s}\left({\hat{\lambda }}_s\right),u_{n-q_s+1}\left({\hat{\lambda }}_s\right),u_{n-q_s+2}\left({\hat{\lambda }}_s\right),\cdots ,u_n\left({\hat{\lambda }}_s\right)$

是线性无关的。因此，我们可以做出如下合理假设。

假设1：对于 $s=1,2,\cdots ,k$，向量组

$v_1\left({\tilde{\lambda }}_s\right),v_2\left({\tilde{\lambda }}_s\right),\cdots ,v_{n-q_s}\left({\tilde{\lambda }}_s\right),u_{n-q_s+1}\left({\tilde{\lambda }}_s\right),u_{n-q_s+2}\left({\tilde{\lambda }}_s\right),\cdots ,u_n\left({\tilde{\lambda }}_s\right)$

是线性无关的。

对于 $s=1,2,\cdots ,k$，定义矩阵

$C\left(\epsilon \right)=\left(\begin{array}{cc}{A}^{sym}+E^{sym}-{\tilde{\lambda }}_s{I}_n& U_{:,n-q_s+1:n}\left({\tilde{\lambda }}_s\right)\\ U_{:,n-q_s+1:n}{\left({\tilde{\lambda }}_s\right)}^{\text{T}}& O_{q_s,q_s}\end{array}\right)$, (3)

其中

${\epsilon }^{sym}={\left({\epsilon }_{1,1},{\epsilon }_{1,2},\cdots ,{\epsilon }_{1,n},{\epsilon }_{2,1},\cdots ,{\epsilon }_{2,n-1},\cdots ,{\epsilon }_{n-1,1},{\epsilon }_{n-1,2},{\epsilon }_{n,1}\right)}^{\text{T}}$. (4)

引理2：对于 $s=1,2,\cdots ,k$，如果 $coran{k}_{\delta }\left(A^{sym}-{\tilde{\lambda }}_s{I}_n\right)=q_s$，则矩阵

$\left(\begin{array}{cc}{A}^{sym}-{\tilde{\lambda }}_s{I}_n& U_{:,n-q_s+1:n}\left({\tilde{\lambda }}_s\right)\\ U_{:,n-q_s+1:n}{\left({\tilde{\lambda }}_s\right)}^{\text{T}}& O_{q_s,q_s}\end{array}\right)$ (5)

是非奇异的。

证明：

$\begin{array}{l}\left(\begin{array}{cc}{A}^{sym}-{\tilde{\lambda }}_s{I}_n& U_{:,n-q_s+1:n}\left({\tilde{\lambda }}_s\right)\\ U_{:,n-q_s+1:n}{\left({\tilde{\lambda }}_s\right)}^{\text{T}}& O_{q_s,q_s}\end{array}\right)\\ =\left(\begin{array}{cc}U\left({\tilde{\lambda }}_s\right)& O_{q_s,q_s}\\ O_{q_s,n}& I_{q_s}\end{array}\right)\left(\begin{array}{ccc}{\Sigma }_{1:n-q_s,1:n-q_s}\left({\tilde{\lambda }}_s\right)& O_{n-q_s,q_s}& O_{n-q_s,q_s}\\ O_{q_s,n-q_s}& {\Sigma }_{n-q_s+1:n,n-q_s+1:n}\left({\tilde{\lambda }}_s\right)& I_{q_s}\\ O_{q_s,n-q_s}& I_{q_s}& O_{q_s,q_s}\end{array}\right)\left(\begin{array}{cc}{V}_{:,1:n-q_s}{\left({\tilde{\lambda }}_s\right)}^{\text{T}}& O_{n-q_s,q_s}\\ U_{:n-q_s+1:n}\left({\tilde{\lambda }}_s\right)& O_{q_s,q_s}\\ {\Sigma }_{n-q_s+1:n,n-q_s+1:n}\left({\tilde{\lambda }}_s\right)& I_{q_s}\end{array}\right)\end{array}$

证毕。 ¨

定义常量

$K=\min \left\{\frac{1}{n{\sqrt{n+q_s\Vert C{\left(0\right)}^{-1}\Vert }}_{\infty }},s=1,2,\cdots ,k\right\}$. (6)

引理3：若 ${\Vert \epsilon \Vert }_{\infty }<K$，则 $C\left(\epsilon \right),s=1,2,\cdots ,k$，非奇异。

证明：若 ${\Vert \epsilon \Vert }_{\infty }<K$，则对 $s=1,2,\cdots ,k$，

$\begin{array}{c}{\Vert I_{n+1}-C{\left(0\right)}^{-1}C\left(\epsilon \right)\Vert }_2\le \sqrt{n+q_s}{\Vert I_{n+q_s}-C{\left(0\right)}^{-1}C\left(\epsilon \right)\Vert }_{\infty }\\ \le \sqrt{n+q_s}{\Vert C{\left(0\right)}^{-1}\Vert }_{\infty }{\Vert C\left(0\right)-C\left(\epsilon \right)\Vert }_{\infty }\\ <n\sqrt{n+q_s}{\Vert \epsilon \Vert }_{\infty }{\Vert C{\left(0\right)}^{-1}\Vert }_{\infty }\\ <1\end{array}$ (7)

由引理1可知， $C\left(\epsilon \right),s=1,2,\cdots ,k$，非奇异。证毕。 ¨

对于 $s=1,2,\cdots ,k$，令 $W_s\left(\epsilon \right)$ 和 $F_s\left(\epsilon \right)$ 分别是线性系统

$C\left(\epsilon \right)\left(\begin{array}{c}{W}_s\left(\epsilon \right)\\ F_s\left(\epsilon \right)\end{array}\right)=\left(\begin{array}{c}{O}_{n,q_s}\\ I_{q_s}\end{array}\right)$ (8)

解的前n行和后 $q_s$ 行。根据克莱姆法则，我们有如下引理。

引理4：对于 $s=1,2,\cdots ,k$，矩阵 $F_s\left(\epsilon \right)$ 是对称矩阵。

命题1：假设 ${\Vert \tilde{\epsilon }\Vert }_{\infty }<K$。若对 $s=1,2,\cdots ,k$，有 $F_s\left(\tilde{\epsilon }\right)=O_{q_s,q_s}$，则

$corank\left(A^{sym}+{\tilde{E}}^{sym}-{\tilde{\lambda }}_s{I}_n\right)=q_s$.

证明：若 ${\Vert \tilde{\epsilon }\Vert }_{\infty }<K$，则由引理3可知 $C\left(\epsilon \right)$， $s=1,2,\cdots ,k$，非奇异。令 $s=1,2,\cdots ,k$，若 $F_s\left(\tilde{\epsilon }\right)=O_{q_s,q_s}$，则

$\begin{array}{l}\left(A^{sym}+E^{sym}-{\tilde{\lambda }}_s{I}_n\right)W_s\left(\tilde{\epsilon }\right)=O_{n,q_s},\\ U_{:,n-q_s+1:n}{\left({\tilde{\lambda }}_s\right)}^{\text{T}}{W}_s\left(\tilde{\epsilon }\right)=I_{q_s}\end{array}$

若 $W_s\left(\tilde{\epsilon }\right)$ 的列向量线性相关，则 $C_s\left(\tilde{\epsilon }\right)$ 奇异，与条件矛盾。因此 $corank\left(A^{sym}+{\tilde{E}}^{sym}-{\tilde{\lambda }}_s{I}_n\right)\ge q_s$。若 $corank\left(A^{sym}+{\tilde{E}}^{sym}-{\tilde{\lambda }}_s{I}_n\right)>q_s$，则存在非零向量

${\beta }_s\in null\left(A^{sym}+{\tilde{E}}^{sym}-{\tilde{\lambda }}_s{I}_n\right)$,

且 ${\beta }_s$ 与 $W_s\left(\tilde{\epsilon }\right)$ 的列向量线性无关。因此存在非零向量 $b_s\in {ℝ}^{q_s+1}$ 使得

$\left(\begin{array}{c}{A}^{sym}+E^{sym}-{\tilde{\lambda }}_s{I}_n\\ U_{:,n-q_s+1:n}{\left({\tilde{\lambda }}_s\right)}^{\text{T}}\end{array}\right)\left(\begin{array}{cc}{W}_s\left(\tilde{\epsilon }\right)& {\beta }_s\end{array}\right)b_s=0$,

矛盾。故 $corank\left(A^{sym}+{\tilde{E}}^{sym}-{\tilde{\lambda }}_s{I}_n\right)=q_s$。 ¨

引理 5：非线性系统

$G\left(\epsilon \right)=\left(\begin{array}{c}{F}_1{\left(\epsilon \right)}_{1:q_s,1}\\ F_1{\left(\epsilon \right)}_{2:q_s,2}\\ ⋮\\ F_1{\left(\epsilon \right)}_{q_s:q_s,q_s}\end{array}\right)=0$ (9)

是欠定非线性系统。

证明：系统(9)包含 $\frac{n\left(n+1\right)}{2}$ 个变量， ${\displaystyle \underset{i=1}{\overset{s}{\sum }}\frac{q_i\left(q_i+1\right)}{2}}$ 个方程。由于

$\begin{array}{c}\frac{q_1^2+q_2^2+\cdots +q_s^2+q_1+q_2+\cdots +q_s}{2}\le \frac{q_1^2+q_2^2+\cdots +q_s^2+n}{2}\\ =\frac{{\left(q_1+q_2+\cdots +q_s\right)}^2-2{\displaystyle \underset{1\le i\le j\le s}{\sum }{q}_i{q}_j}+n}{2}\\ \le \frac{n^2+n}{2}-{\displaystyle \underset{1\le i\le j\le s}{\sum }{q}_i{q}_j}\\ \le \frac{n\left(n+1\right)}{2}\end{array}$

则引理得证。 ¨

由引理3可知，存在一个以原点为中心的邻域，使得对于该邻域内的任意 $\epsilon$，对 $s=1,2,\cdots ,k$，矩阵 $C\left(\epsilon \right)$ 非奇异。因此在该邻域内，系统(8)解向量的每一项都充分可微。

对于每一对 $\left(i,j\right),1\le i\le n-1,1\le j\le n-i+1$，对系统(8)两端关于 ${\epsilon }_{i,j}$ 求偏导得到

$\left(\begin{array}{cc}{A}^{sym}+E^{sym}-{\tilde{\lambda }}_s{I}_n& U_{:,n-q_s+1:n}\left({\tilde{\lambda }}_s\right)\\ U_{:,n-q_s+1:n}{\left({\tilde{\lambda }}_s\right)}^{\text{T}}& O_{q_s,q_s}\end{array}\right)\left(\begin{array}{c}\frac{\partial W_s\left(\epsilon \right)}{\partial {\epsilon }_{i,j}}\\ \frac{\partial F_s\left(\epsilon \right)}{\partial {\epsilon }_{i,j}}\end{array}\right)=t\left(-\left(\begin{array}{c}{O}_{j+i-2,q_s}\\ W_s{\left(\epsilon \right)}_{i,:}\\ O_{n+q_s-j-i+1,q_s}\end{array}\right)-\left(\begin{array}{c}{O}_{i-1,q_s}\\ W_s{\left(\epsilon \right)}_{j+i-1,:}\\ O_{n+q_s-i,q_s}\end{array}\right)\right)$, (10)

其中

$t=\{\begin{array}{l}\frac{1}{2},j=1\\ 1,2\le j\le n-i+1\end{array}$.

对于每一对 $\left(i_1,j_1\right),\left(i_2,j_2\right),1\le i_1\le n,\text{1}\le j_1\le n-i_1+1,1\le i_2\le n,\text{1}\le j_2\le n-i_2+1$，对系统(10)两端关于 ${\epsilon }_{i_1,j_1},{\epsilon }_{i_2,j_2}$ 求偏导得

$\begin{array}{l}\left(\begin{array}{cc}{A}^{sym}+E^{sym}-{\tilde{\lambda }}_s{I}_n& U_{:,n-q_s+1:n}\left({\tilde{\lambda }}_s\right)\\ U_{:,n-q_s+1:n}{\left({\tilde{\lambda }}_s\right)}^{\text{T}}& O_{q_s,q_s}\end{array}\right)\left(\begin{array}{c}\frac{{\partial }^{(2)}{W}_s\left(\epsilon \right)}{\partial {\epsilon }_{i_1,j_1}\partial {\epsilon }_{i_2,j_2}}\\ \frac{{\partial }^{(2)}{F}_s\left(\epsilon \right)}{\partial {\epsilon }_{i_1,j_1}\partial {\epsilon }_{i_2,j_2}}\end{array}\right)\\ =t_1\left(-\left(\begin{array}{c}{O}_{j_1+i_1-2,q_s}\\ {\left(\frac{\partial W_s\left(\epsilon \right)}{\partial {\epsilon }_{i_2,j_2}}\right)}_{i_1,:}\\ O_{n+q_s-j_1-i_1+1,q_s}\end{array}\right)-\left(\begin{array}{c}{O}_{i_1-1,q_s}\\ {\left(\frac{\partial W_s\left(\epsilon \right)}{\partial {\epsilon }_{i_2,j_2}}\right)}_{j_1+i_1-1,:}\\ O_{n+q_s-i_1,q_s}\end{array}\right)\right)+t_2\left(-\left(\begin{array}{c}{O}_{j_2+i_2-2,q_s}\\ {\left(\frac{\partial W_s\left(\epsilon \right)}{\partial {\epsilon }_{i_1,j_1}}\right)}_{i_2,:}\\ O_{n+q_s-j_2-i_2+1,q_s}\end{array}\right)-\left(\begin{array}{c}{O}_{i_2-1,q_s}\\ {\left(\frac{\partial W_s\left(\epsilon \right)}{\partial {\epsilon }_{i_1,j_1}}\right)}_{j_2+i_2-1,:}\\ O_{n+q_s-i_2,q_s}\end{array}\right)\right)\end{array}$ (11)

其中

$t_1=\{\begin{array}{l}\frac{1}{2},j_1=1\\ 1,2\le j_1\le n-i_1+1\end{array}$, $t_2=\{\begin{array}{l}\frac{1}{2},j_2=1\\ 1,2\le j_2\le n-i_2+1\end{array}$.

当 $\epsilon =0$ 时，求解(10)和(11)，得到雅可比矩阵 $G^{\prime }\left(0\right)$ 和海森矩阵 $G^{″}\left(0\right)$。

假设 2：雅克比矩阵 $G^{\prime }\left(0\right)$ 行满秩。

假设对于指标集 $\mathcal{I}=\left\{\left(i_1,j_1\right),\left(i_2,j_2\right),\cdots ,\left(i_m,j_m\right)\right\}$，矩阵

$\left(\frac{\partial G\left(0\right)}{\partial {\epsilon }_{i_1,j_1}}\frac{\partial G\left(0\right)}{\partial {\epsilon }_{i_2,j_2}}\cdots \frac{\partial G\left(0\right)}{\partial {\epsilon }_{i_m,j_m}}\right)$ (12)

非奇异。由

$\tilde{G}\left({\epsilon }_{i_1,j_1},{\epsilon }_{i_2,j_2},\cdots ,{\epsilon }_{i_m,j_m}\right)={G\left(\epsilon \right)|}_{{\epsilon }_{i,j}=\text{0}\text{for}\left(i,j\right)\notin \mathcal{I}}$ (13)

定义非线性系统

$\tilde{G}\left({\epsilon }_{i_1,j_1},{\epsilon }_{i_2,j_2},\cdots ,{\epsilon }_{i_m,j_m}\right)=0$. (14)

注2：利用定理2，我们计算了当零向量为非线性系统(14)的解时的验证误差界。

对于系统(14)，定义定理2中的常量B， $\eta$，

$B={\Vert \widetilde{G^{\prime }}{\left(0\right)}^{-1}\Vert }_{\infty },\eta ={\Vert {\tilde{G}}^{\prime }{\left(0\right)}^{-1}\tilde{G}\left(0\right)\Vert }_{\infty }$. (15)

定义m维区间向量 $\Omega ={\left({\Omega }_{i_1,j_1},{\Omega }_{i_2,j_2},\cdots ,{\Omega }_{i_m,j_m}\right)}^{\text{T}}$，其每一项均为区间 $\left[-2\eta ,2\eta \right]$。定义同类型的区间摄动矩阵 $\tilde{\Omega }$，其满足当 $\left(i,j\right)\in \mathcal{I}$ 时， ${\tilde{\Omega }}_{i,j}={\Omega }_{i,j}$，当 $\left(i,j\right)\notin \mathcal{I}$ 时， ${\tilde{\Omega }}_{i,j}=0$。

利用verifylss函数，求解当 $E=\tilde{\Omega }$ 时相应的区间线性系统(8) (10) (11)，可得满足下式

$H_{s,t}\supset \frac{{\partial }^{\left(2\right)}G\left(\Omega \right)}{\partial {\epsilon }_{i_s,j_s}\partial {\epsilon }_{i_t,j_t}},1\le s\le t\le m$, (16)

的区间向量集合 $\left\{H_{s,t}:1\le s\le t\le m\right\}$。定义Lipschitz常量

$\kappa =\underset{1\le s\le t\le m}{\max }\max \left\{{\Vert H\Vert }_{\infty }:\forall H\in H_{s,t}\right\}$. (17)

4. 主要算法

算法1

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

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

输出 $\left\{{\tilde{\lambda }}_s\in ℝ:s=1,2,\cdots ,k\right\}$， $\left\{q_s\in ℝ:s=1,2,\cdots ,k\right\}$， $\rho \in ℝ$ 和指标集 $\mathcal{I}$。

步骤1利用Matlab中eig命令计算 $A^{sym}$ 所有互异特征值 $\left\{{\tilde{\lambda }}_s\in ℝ:s=1,2,\cdots ,k\right\}$。

步骤2对 $s=1,2,\cdots ,k$，计算 $A^{sym}-{\tilde{\lambda }}_s{I}_n$ 的奇异值分解。

步骤3求解(10)得到雅可比矩阵 $G^{\prime }\left(0\right)$。

步骤4若 $G^{\prime }\left(0\right)$ 行满秩，选择指标集 $\mathcal{I}$ 使得 $\widetilde{G^{\prime }}\left(0\right)$ 非奇异。

步骤5利用(15)计算常量B， $\eta$。

步骤6利用(17)计算常量 $\kappa$。

步骤7若 $h=B\eta \kappa \le 1/2$，则利用(1)式计算 $\rho$。

步骤8若 $\rho <K$，返回 $\left\{{\tilde{\lambda }}_s:s=1,2,\cdots ,k\right\}$， $\left\{q_s:s=1,2,\cdots ,k\right\}$， $\rho$ 和 $\mathcal{I}$。

定理 3：如果算法1成功，则存在同类型的摄动矩阵 ${\tilde{E}}^{sym}$，其满足条件

$\begin{array}{l}\left|{\tilde{\epsilon }}_{i,j}\right|\le \rho ,\left(i,j\right)\in \mathcal{I};\\ {\tilde{\epsilon }}_{i,j}=0,\left(i,j\right)\notin \mathcal{I}\end{array}$ (18)

$\left\{{\tilde{\lambda }}_s,s=1,2,\cdots ,k\right\}$ 是矩阵 $A^{sym}+{\tilde{E}}^{sym}$ 的精确谱。此外，对于 $s=1,2,\cdots ,k$， $q_s$ 是特征值 ${\tilde{\lambda }}_s$ 的几何重数。

证明若算法1成功，则由定理2可知，对于满足条件(18)的同类型摄动矩阵 ${\tilde{E}}^{sym}$，有 $\tilde{G}\left({\epsilon }_{i_1,j_1},{\epsilon }_{i_2,j_2},\cdots ,{\epsilon }_{i_m,j_m}\right)=0$。因此，存在一个满足条件(18)的区间矩阵 ${\tilde{E}}^{sym}$，使得 $G\left(\tilde{\epsilon }\right)=0$。最后，由命题1可知，对于 $s=1,2,\cdots ,k$， $corank\left(A^{sym}+{\tilde{E}}^{sym}-{\tilde{\lambda }}_s{I}_n\right)=q_s$。 ¨

5. 应用实例

在本节中，演示了Verifylss算法的性能。在Windows 7下，使用Matlab R2012a和INTLAB V5进行了以下实验。

例1给定对称矩阵

$A=\left[\begin{array}{ccccccc}1& 0& 0& 1& 0& 0& 0\\ 0& 1& 1& 0& 0& 0& 0\\ 0& 1& 1& 0& 0& 0& 0\\ 1& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 2\\ 0& 0& 0& 0& 0& 2& 0\\ 0& 0& 0& 0& 2& 0& 1\end{array}\right]$,

由算法1得到 $q_1=q_4=1,q_2=2,q_3=3$， $\left\{{\tilde{\lambda }}_1=-1,{\tilde{\lambda }}_2=0,{\tilde{\lambda }}_3=2,{\tilde{\lambda }}_4=3\right\}$， $\rho =0$。

例2对于对称矩阵

$A=\left[\begin{array}{cccccccc}1& 1& 0& 1& 0& 0& 0& 0\\ 1& 1& 1& 0& 0& 0& 0& 0\\ 0& 1& 1& 0& 0& 0& 0& 0\\ 1& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 2\\ 0& 0& 0& 0& 0& 1& 2& 0\\ 0& 0& 0& 0& 0& 2& 1& 0\\ 0& 0& 0& 0& 2& 0& 0& 1\end{array}\right]$,

应用算法1得到 $q_1=q_6=2,q_2=q_3=q_4=q_5=1$， $\left\{{\tilde{\lambda }}_1=-1.0000,{\tilde{\lambda }}_2=-0.6180,{\tilde{\lambda }}_3=0.3820,{\tilde{\lambda }}_4=1.6180,{\tilde{\lambda }}_5=2.6180,{\tilde{\lambda }}_6=3.0000\right\}$， $\rho =5.5359e-16$。

例3给定一个对称矩阵

$A=\left[\begin{array}{ccccccccc}1& 0& 0& 1& 1& 0& 0& 0& 0\\ 0& 1& 1& 1& 0& 0& 0& 0& 0\\ 0& 1& 1& 0& 0& 0& 0& 0& 0\\ 1& 1& 0& 1& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 1& 1& 0\\ 0& 0& 0& 0& 0& 0& 1& 1& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 1\end{array}\right]$,

应用算法1得到 $q_2=q_4=3,q_1=q_3=q_5=1$， $\left\{{\tilde{\lambda }}_1=-0.7321,{\tilde{\lambda }}_2=0,{\tilde{\lambda }}_3=1.0000,{\tilde{\lambda }}_4=2.0000,{\tilde{\lambda }}_5=2.7321\right\}$， $\rho =5.1876e-16$。

Table 1. The calculation results of ρ in Example 4

表1. 例4中 $\rho$ 的计算结果

例4设 $n\times n$ 的对称矩阵A，其元素为服从区间[0, 1]的均匀分布。表1给出了对于不同n的取值，算法1计算出的可信半径 $\rho$。

基金项目

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

NOTES

^*通讯作者。