在控制理论、随机滤波和超分辨率图像恢复等领域, 会遇到形如
的含有未知矩阵逆幂的非线性矩阵方程[ 1-6], 其中$E_i, F_i, G, X\in{ \bf{R}}^{n\times n}$(实矩阵集合). 针对方程(1)的特例$X+F^{\mathrm{T}}_1X^{-1}F_1=G$, 陈小山等[ 1]给出了方程存在唯一对称正定解的充分条件, 以及对称正定解和最大解的存在区间, 并利用微分方法给出了最大解的一阶扰动界; 针对方程(1)的特例$X+F^{\mathrm{T}}_2X^{-2}F_2=I$, Zhao 等[ 2]研究了$F_2$可逆时, 方程存在对称正定解的充要条件; Peng 等[ 3]给出了当$\alpha\geq1$时求方程$X+A^{\mathrm{T}}X^{-\alpha}A=Q$的最小对称正定解的不动点迭代算法, 并讨论了算法的收敛性问题; 李静等[ 4-5]给出了方程$X-A^{\mathrm{T}}X^{-k}A=Q (k\geq1)$存在对称正定解的充分条件和解的性质, 并建立了两种迭代算法; Ran 等[ 6]研究了一般非线性矩阵方程$X+A^{\mathrm{T}}f(X)A=Q$存在唯一半正定解的充分条件, 并给出了半正定解的扰动分析. 将方程(1)中的$X^{-1}$替换为$X$可得
本文建立求方程(2)的对称解的双迭代算法.
首先运用牛顿算法求方程(2)的对称解, 然后采用修正共轭梯度法(MCG算法)求由牛顿算法每一步迭代计算 导出的线性矩阵方程(LME)的对称解或者对称最小二乘解(Ls解), 建立求方程(2)的对称解的双迭代算法. 通常的牛顿算法要求每一步的LME总有唯一对称解, 因此需要对方程(2)的系数矩阵和常数项矩阵附加较强的限定, 否则算法就会中断. 文中采用的MCG算法不同于通常的共轭梯度法, 它不要求涉及的线性代数方程组的系数矩阵和常数项矩阵 对称正定、可逆或者列满秩, 因此总是可行的. 文中建立的双迭代算法仅要求方程(2)有对称解, 不要求它的对称解唯一, 也不对它的系数矩阵和常数项矩阵做附加限定.
用${\bf{SR}}^{n\times n}$表示$n$阶实对称矩阵集合, 引入以下矩阵函数:
当$X^{-1}Y$的谱半径小于1时($Y$的范数较小时能够满足), 可以导出
这里$\phi_{X}(Y)$是$\psi(X)$在“点”$X$沿着“方向”$Y$的Fréchet导数.
引理 设$X\in{ \bf{SR}}^{n\times n}$是方程(2)的近似解, 那么, 求方程(2)的对称解等价于求校正值 $Y\in{ \bf{SR}}^{n\times n}$使得$\psi(X+Y)=O$, 并可以线性化为求$Y\in{ \bf{SR}}^{n\times n}$使得$\phi_{X}(Y)=-\psi(X)$.
证 设$X^*\in{ \bf{SR}}^{n\times n}$是方程(2)的精确解, 已知$X\in{ \bf{SR}}^{n\times n}$是方程(2)的近似解, 令$X^*=X+Y$, 那么, 求方程(2)的解$X^*\in{ \bf{SR}}^{n\times n}$等价于求校正值 $Y\in{ \bf{SR}}^{n\times n}$使得$\psi(X+Y)=O$. 这是关于$Y$的非线性矩阵方程, 根据牛顿算法的基本原理, 当$Y$的范数较小时, 舍去式(3)右端关于$Y$的高次项$\gamma_{X}(Y)$, 即用线性部分近似可得 $\psi(X+Y)\approx\psi(X)+\phi_{X}(Y).$ 于是, 求方程(2)的解$X^*\in{ \bf{SR}}^{n\times n}$, 可近似的转化为求关于$Y$的线性矩阵方程$\phi_{X}(Y)=-\psi(X)$的对称解.
上述引理中, $\phi_{X}(Y)=-\psi(X)$的解$Y\in{ \bf{SR}}^{n\times n}$一般不是$\psi(X+Y)=O$的精确解, 从而由$X^*=X+Y$确定的$X^*\in{ \bf{SR}}^{n\times n}$也不是$\psi(X)=O$的精确解, 但是可以看作一个近似解. 借鉴文献[7], 通过修改某些矩阵的类型, 建立求方程(2)的对称解的牛顿算法如下.
第1步: 给定初始矩阵$X^{(1)}\in{ \bf{SR}}^{n\times n}$, 置$k:=1$.
第2步: 如果$\psi(X^{(k)})=O$, 停止; 否则, 求$Y^{(k)}\in{ \bf{SR}}^{n\times n}$, 使得
第3步: 计算$X^{(k+1)}=X^{(k)}+Y^{(k)}$, 置$k:=k+1$, 转第2步.
需要指出: 对于牛顿算法中的某个$k$, 当LME(4)没有对称解$Y$时, 可用它的对称Ls解来代替, 这也是本文研究的双迭代算法的一个特点. 对于牛顿算法有如下的收敛性结论[7]: 假设$X^*$是方程(2)的单根, 且初始矩阵$X^{(1)}$充分接近于$X^*$, 那么由牛顿算法确定的矩阵序列$\{X^{(k)}\}$收敛于$X^*$.
用$A\otimes B$表示矩阵$A$与$B$的Kronecker积, $\overline{\mathrm{vec}}(A)$表示将矩阵$A$按行拉直构成的列向量, $\|A\|$表示矩阵$A$的Frobenius范数. 下面建立求LME(4)的对称解与对称Ls解的MCG算法. 考虑LME(4)的一般形式
其中$A_i, B_i, F\in{ \bf{R}}^{n\times n}$. 研究以下两个问题:
问题 Ⅰ 设LME(5)有对称解, 求$Y\in{ \bf{SR}}^{n\times n}$, 使满足LME(5).
问题 Ⅱ 设LME(5)无对称解, 求$Y\in{ \bf{SR}}^{n\times n}$, 使得$\|\sum\limits_{i=1}^7A_{i}YB_{i}-F\|=\mathrm{min}$.
3.1 求解问题Ⅰ的MCG算法
借鉴文献[8]的算法原理, 通过修改有关矩阵的类型或者算法, 建立求解问题Ⅰ的MCG算法. 引进记号: $u(Y)=\sum\limits_{i=1}^7A_{i}YB_{i},$ $w(R)=\sum\limits_{i=1}^7A_{i}^{\mathrm{T}}RB_{i}^{\mathrm{T}}.$
算法 1(问题Ⅰ的MCG算法)
第1步 任意给定初始矩阵$Y^{(1)}\in{ \bf{SR}}^{n\times n}$, 置$k:=1$, 计算 $R_{k}=F-u(Y^{(k)}),$ $ \widetilde{R}_{k}=w(R_{k}),$ $Z_{k}=\frac{1}{2}(\widetilde{R}_{k}+{\widetilde{R}_{k}}^{\mathrm{T}});$
第2步 如果$R_{k}=O$, 或者$R_{k}\neq O$而$Z_{k}=O$, 停止计算; 否则, 计算 $Y^{(k+1)}=Y^{(k)}+\frac{\|R_{k}\|^2}{\|Z_{k}\|^2}Z_{k}; $
第3步 计算
第4步 置$k:=k+1$, 转第2步.
可以验证, 算法1中的矩阵满足$Y^{(k)}\in{ \bf{SR}}^{n\times n}$, 且有如下的收敛定理(证明过程类似文献[ 8]).
定理 1 设问题Ⅰ相容(指LME(5)有对称解), 则对任意初始矩阵$Y^{(1)}\in{ \bf{SR}}^{n\times n}$, 算法1可在有限步计算后求得问题Ⅰ的一个解, 即LME(5)的一个对称解; 若取初始矩阵$Y^{(1)}$满足 $Y^{(1)}=\frac{1}{2}[w(H)+(w(H))^{\mathrm{T}}]$(任意$H\in{ \bf{R}}^{n\times n}$), 则算法1可在有限步计算后得到问题Ⅰ的唯一极小范数解, 即LME(5)的唯一极小范数对称解; 问题Ⅰ不相容的充要条件是存在正整数$k$, 使得由算法1得到的$R_{k}\neq O$而$Z_{k}=O$.
3.2 求解问题Ⅱ的MCG算法
根据定理1, 在算法1中, 当$R_{k}\neq O$而$Z_{k}=O$时, 算法1中断, 这表明LME(5)无对称解. 因此, 需要求解问题Ⅱ, 即求LME(5)的对称Ls解. 下面通过构造约束正规矩阵方程, 将求LME(5)的对称Ls解问题转化为求约束正规矩阵方程的对称解的问题. 然后参照算法1, 建立求LME(5)的对称Ls解的MCG算法. 引进记号 $Q=w(F)+(w(F))^{\mathrm{T}},$ $g(Y)=w(u(Y))+[w(u(Y^{\mathrm{T}}))]^{\mathrm{T}}.$
定理 2 求解问题Ⅱ等价于求LME(约束正规矩阵方程)
的对称解, 并且LME(6)一定有对称解.
证 求解问题Ⅱ等价于求$Y\in{ \bf{SR}}^{n\times n}$, 使得
下面证明求解极小值问题(7)等价于求LME(6)的对称解. 令 $M=\left( \begin{array}{c} {\sum\limits_{i=1}^7A_{i}\otimes B_{i}^{\mathrm{T}}} \\ {(\sum\limits_{i=1}^7A_{i}\otimes B_{i}^{\mathrm{T}})T_{n,n}} \end{array}\right),$ $f=\overline{\mathrm{vec}}\left( \begin{array}{c} F\\ F \end{array}\right),$ $y=\overline{\mathrm{vec}}(Y),$ 其中$T_{m,n}$表示满足$\overline{\mathrm{vec}}(A^{\mathrm{T}}_{m\times n})=T_{m,n}\overline{\mathrm{vec}}(A)$的$mn$阶置换矩阵. 将线性矩阵方程组
按行拉直可得线性代数方程组$My=f$, 求它的Ls解等价于求线性矩阵方程组(8)的Ls解, 即求极小值问题(7)的解. 方程组$My=f$的正规方程组为$M^{\mathrm{T}}My=M^{\mathrm{T}}f$, 还原为矩阵形式就是LME(6). 因为求方程组$My=f$的Ls解, 就是求它的正规方程组$M^{\mathrm{T}}My=M^{\mathrm{T}}f$的解, 即求LME(6)的解, 所以求解问题Ⅱ等价于求LME(6)的对称解.
下面证明LME(6)有对称解. 因为正规方程组$M^{\mathrm{T}}My=M^{\mathrm{T}}f$有解, 所以LME(6)有解. 设$\widetilde{Y}$是LME(6)的一个解(未必是对称解), 那么$g(\widetilde{Y})=Q$. 令$\widetilde{Y}^*=\frac{1}{2}(\widetilde{Y}+\widetilde{Y}^{\mathrm{T}})$, 则 $\widetilde{Y}^*\in{ \bf{SR}}^{n\times n}$, 且有$g(\widetilde{Y}^*)=Q$, 故LME(6)有对称解.
参照算法1以及文献[ 8-9]的算法原理, 建立求LME(6)的对称解, 即求解问题Ⅱ的MCG算法如下.
算法 2(问题Ⅱ的MCG算法)
第1步 任意给定初始矩阵$Y^{(1)}\in{ \bf{SR}}^{n\times n}$, 置$k:=1$, 计算 $R_{k}=Q-g(Y^{(k)}),$ $\widetilde{R}_{k}=g(R_{k}),$ $Z_{k}=\widetilde{R}_{k};$
第2步 如果$R_{k}=O$, 停止计算; 否则, 计算 $Y^{(k+1)}=Y^{(k)}+\frac{\|R_{k}\|^2}{\|Z_{k}\|^2}Z_{k};$
第3步 计算 $R_{k+1}=Q-g(Y^{(k+1)}),$ $\widetilde{R}_{k+1}=g(R_{k+1}),$ $Z_{k+1}=\widetilde{R}_{k+1}+\frac{\|R_{k+1}\|^2}{\|R_{k}\|^2}Z_{k};$
可以验证, 算法2中的矩阵满足$Y^{(k)}\in{ \bf{SR}}^{n\times n}$, 且有如下的收敛定理(证明过程类似文献[ 9]).
定理 3 对任意的初始矩阵$Y^{(1)}\in{ \bf{SR}}^{n\times n}$, 算法2可在有限步计算后求得问题Ⅱ的一个解, 即LME(5)的一个对称Ls解; 若取初始矩阵$Y^{(1)}$满足 $Y^{(1)}=g(H)$ (任意$H\in{ \bf{SR}}^{n\times n}$), 则算法2可在有限步计算后得到LME(6)的唯一极小范数对称解, 也就是LME(5)的唯一极小范数对称Ls解.
求方程(2)的对称解, 可采用以下两种计算方案.
方案一
第1步 给定初始矩阵$X^{(1)}\in{ \bf{SR}}^{n\times n}$, 置$k:=1$;
第2步 如果$\psi(X^{(k)})=O$, 停止; 否则, 采用算法1求$Y^{(k)}\in{ \bf{SR}}^{n\times n}$, 使满足LME(4); 若算法1中断(此时LME(4)没有对称解), 采用算法2求$Y^{(k)}\in{ \bf{SR}}^{n\times n}$, 使得
第3步 计算$X^{(k+1)}=X^{(k)}+Y^{(k)}$, 置$k:=k+1$, 转第2步.
方案二
第2步 如果$\psi(X^{(k)})=O$, 停止; 否则, 采用算法2求$Y^{(k)}\in{ \bf{SR}}^{n\times n}$, 使满足式(9). 当LME(4)有对称解时, 它的对称Ls解就是它的对称解;
下面的方案一(a)是指算法2的初始矩阵取零矩阵, 方案一(b)是指算法2的初始矩阵取算法1中断前得到的矩阵, 用 $k_{12}$表示方案一中算法1的中断次数, $k_0$、$k_1$及$k_2$分别表示牛顿算法、算法1及算法2的总迭代次数, $t$表示计算时间, $n$表示未知矩阵的阶数. 算例使用Matlab7.0软件-CPU3.00GHz微机计算, 取牛顿算法的终止准则为 $10^{-7}$, MCG算法的终止准则为$10^{-8}$, 没有特别约定时MCG算法的初始矩阵为$Y^{(1)}=O$.
例 1 采用双迭代算法和Li's算法(文献[ 4]算法)求方程(2)的特例$X^{-1}-F_3^{\mathrm{T}}X^3F_3=G$的对称解.
(1) 系数矩阵和常数项矩阵如下(取自文献[ 4]):
根据文献[ 4], 取初始对称矩阵为$X^{(1)}=\frac{5}{6}I$, 迭代次数和计算时间见表 1.
由方案一和方案二求得方程(2)的对称解均为
由Li's算法求得对应的方程(1)的对称解$X^{(7)}$就是$(X^{(5)})^{-1}$.
(2) 系数矩阵和常数项矩阵如下: $F_3=\left( {\begin{array}{ccc} 0.3&0.1&0.7\\ 0.1&0.2&0.5\\ 0.3&0.1&0.4 \end{array}} \right), G=\mathrm{ones}(3).$ 取初始对称矩阵为$X^{(1)}=\frac{2}{3}I$, 迭代次数和计算时间见表 2.
由方案一和方案二求得方程(2)的对称解均为 $X^{(8)}=\left( {\begin{array}{ccc} 1.7668&-0.6312&-0.6261\\ -0.6312&1.9587&-0.6773\\ -0.6261&-0.6773&1.1429 \end{array}} \right),$ Li's算法失效的原因是常数项矩阵不满足文献[ 4]中算法的要求.
需要指出, 双迭代算法不对方程(2)的系数矩阵和常数项矩阵做附加限定, 不必选取特定的初始对称矩阵, 但求出的只是方程(2)的对称解, 不一定具备正定性. 而Li's算法对方程(2)的常数项矩阵有要求, 选取特定的初始对称矩阵时, 可求出方程(2)的对称正定解.
例 2 采用两种计算方案求方程(2)的对称解, 其系数矩阵和常数项矩阵等如下:
令$X^{(0)}=\mathrm{diag}(\widetilde{X}^{(0)},\cdots,\widetilde{X}^{(0)}), X^{(1)}=\mathrm{diag}(\widetilde{X}^{(1)},\cdots,\widetilde{X}^{(1)}),$ 构造矩阵
那么方程(2)有对称解, 从而方程(1)有对称解, 取牛顿算法的初始矩阵为$X^{(1)}$.
(1) 选取子矩阵$\widetilde{X}^{(1)}=U_1$, 两种方案的迭代次数和计算时间见 表 3.
两种方案求出方程(2)的对称解均为已知的$X^{(0)}$, 对应方程(1)的对称解为$(X^{(0)})^{-1}$.
(2) 选取子矩阵$\widetilde{X}^{(1)}=U_2$, 两种方案的迭代次数和计算时间见表 4.
两种方案求出方程(2)的对称解(未知矩阵阶数$n=6$的情形)均为
对应方程(1)的对称解为$(X^{(5)})^{-1}$. 对比(1)中的计算结果可见, 当牛顿算法的初始矩阵不同时, 求得方程(2)的对称解也不同(比较(1)中的$X^{(0)}$与(2)中的$X^{(5)}$), 这说明方程(2)的对称解不唯一, 从而方程(1)的对称解不唯一.
多个算例的计算结果表明 当某一步的LME(4)有对称解时, 方案一一般比方案二的效率高, 其原因是当LME(4)有对称解时, 方案二总是采用求LME(4)的对称Ls解的算法2进行计算, 而算法2比算法1的计算耗时多.
作变量替换可将方程(1)转化为方程(2), 运用牛顿算法求方程(2)的对称解, 并采用MCG算法求由牛顿算法每一步迭代计算 导出的LME的对称解或者对称最小二乘解, 建立了求方程(2)的对称解的双迭代算法, 算例表明双迭代算法是有效的. 通过修改牛顿算法、算法1和算法2中初始矩阵的类型, 以及算法1中$Z_1$与$Z_{k+1}$的计算公式和算法2中涉及的LME(6)的构造方式, 可以建立求方程(2)的其它特殊解的双迭代算法. 需要指出, 本文建立的双迭代算法求出的只是方程(2)的一个对称解, 不一定具备正定性, 如何修改算法使求出的对称解具备正定性, 尚需进一步探讨.