用${R^{m \times n}}$表示$m \times n$实矩阵集合, ${A}\otimes{B}$表示矩阵${A}$与${B}$的Kronecker积, $\overline {{{\hbox{vec}}}}{\rm{(}}{A}{\rm{)}}$表示将矩阵${A}$按行拉直构成的列向量.定义同阶矩阵${A}$与${B}$的内积为$[{A}, {B}] = {\hbox{tr}}({{A}^T}{B})$, 由此导出矩阵的Frobenius范数$\left\|{A} \right\| = \sqrt {\, [{A}, {A}]}$.设$P_i\in R^{n \times n}$为对称正交矩阵, 若$X\in R^{n \times n}$满足$P_1 XP_2=X$, 则称$X$为关于$P_1, P_2$的广义自反矩阵, 记关于$P_1, P_2$的广义自反矩阵集合为$\Omega(P_1, P_2)$.特别地, 当$P_1=P_2$时, 相应的广义自反矩阵是通常的自反矩阵; 当$P_1=P_2$为次单位矩阵时, 相应的广义自反矩阵是中心对称矩阵; 当$P_1=P_2$为单位矩阵时, 相应的广义自反矩阵是一般矩阵.

考虑Riccati矩阵方程

$ AXB+CX^{T}D+\gamma(X)=E_5, $

(1)

其中$A, B, C, D, E_i, X\in R^{n \times n}$,

$ \gamma(X)=XE_1 X+XE_2 X^{T}+X^{T}E_3 X+X^{T}E_4 X^{T}. $

当$E_1=\cdots=E_4=O$时, 方程(1) 是线性矩阵方程, 这类方程在Hamilton力学系统研究和自动控制理论中有重要应用^[1-2].当$E_1=E_3=E_4=O$时, 方程(1) 是一类特殊的Riccati矩阵方程, 这类Riccati矩阵方程在带有测量误差的状态协方差配置问题等工程领域有广泛应用^[3].

目前, 中外学者对一些非线性矩阵方程的求解问题进行了许多研究, 并建立了多种不同的算法.例如, Dehghan等^[4]在求解广义Sylvester矩阵方程组的应用中, 提出了基于牛顿算法求二次矩阵方程特殊解的迭代算法; Long等^[5]建立了求解二次矩阵方程的改进牛顿算法, 并给出了算法的收敛性结论.迄今为止, 关于方程(1) 的广义自反解的迭代算法研究成果尚未见到.本文运用牛顿算法求方程(1) 的广义自反解, 并采用修正共轭梯度法(MCG算法)^[6-9]求导出的线性矩阵方程(LME)的广义自反解或者广义自反最小二乘解(Ls解), 建立求方程(1) 的广义自反解的双迭代算法.

为了书写简洁, 引进记号:

$ \begin{array}{l}\psi(X)=AXB+CX^{T}D+\gamma(X)-E_5,\\ \phi_X(Y)=AYB+CY^{T}D+(XE_1+X^T E_3)Y+Y(E_1 X+E_2 X^T)\\\;\;\;\;\;\;\;\;\;\;\;\;+(XE_2+X^T E_4)Y^T+Y^T(E_3 X+E_4 X^T), \end{array} $

容易导出

$ \psi(X+Y)=\psi(X)+\phi_X(Y)+\gamma(Y), $

(2)

这里$\phi_X(Y)$是$\psi(X)$在“点”$X$沿着“方向”$Y$的导数^[5].

引理1  设$X\in\Omega(P_1, P_2)$是方程(1) 的近似解，那么求方程(1) 的广义自反解等价于求校正值$Y\in\Omega(P_1, P_2)$使得$\psi(X+Y)=O$, 并可线性化为求$Y\in\Omega(P_1, P_2)$使得

$ \phi_X(Y)=-\psi(X). $

(3)

证  已知方程(1) 的近似解$X\in\Omega(P_1, P_2)$时, 令$X^{*}=X+Y$, 那么求方程(1) 的解$X^{*}\in\Omega(P_1, P_2)$, 等价于求校正值$Y\in\Omega(P_1, P_2)$使得$\psi(X+Y)=O$.该式是关于未知矩阵$Y$的二次矩阵方程, 根据牛顿算法的基本原理, 当$Y$的范数较小时, 舍去式(2) 等号右端关于$Y$的二次项$\gamma(Y)$, 即用线性部分近似可得

$ \psi(X+Y)\approx\psi(X)+\phi_X(Y). $

于是求方程(1) 的解$X^{*}\in\Omega(P_1, P_2)$, 可近似的转化为求关于$Y$的线性矩阵方程$\psi(X)+\phi_X(Y)=O$的解$Y\in\Omega(P_1, P_2)$, 即求$Y\in\Omega(P_1, P_2)$使得$\phi_X(Y)=-\psi(X)$.

借鉴文献[4-5]的算法原理, 通过修改有关矩阵的类型, 建立求方程(1) 广义自反解的牛顿算法如下.

第1步  给定初始矩阵$X^{(1)}\in\Omega(P_1, P_2)$, 置$k:\, = 1$;

第2步  如果$\psi(X^{(k)})=O$, 停止; 否则, 求$Y^{(k)}\in\Omega(P_1, P_2)$, 使得

$ \phi_{X^{(k)}}(Y^{(k)})=-\psi(X^{(k)}); $

第3步  计算$X^{(k+1)}=X^{(k)}+Y^{(k)}$, 置$k:\, = k+1$, 转第2步.

对于牛顿算法有以下的收敛性结论^[5]:假设$X^{*}\in R^{n \times n}$是方程(1) 的单根, 且初始矩阵$X^{(1)}$充分接近于$X^{*}$, 那么由牛顿算法确定的矩阵序列$\{X^{(k)}\}$二次收敛于$X^{*}$.

下面建立求LME(3) 广义自反解与广义自反Ls解的MCG算法.考虑LME(3) 的一般形式$(A_i, B_i, C_i, D_i, F, Y\in R^{n \times n})$,

$ \sum\limits_{i=1}^{3}(A_iYB_i+C_iY^TD_i)=F. $

(4)

问题Ⅰ   LME(4) 有广义自反解时, 求$Y\in\Omega(P_1, P_2)$, 使满足LME(4).

问题Ⅱ  LME(4) 无广义自反解时, 求$Y\in\Omega(P_1, P_2)$, 使得

$ \|\sum\limits_{i=1}^{3}(A_iYB_i+C_iY^TD_i)-F\|=\min. $

(5)

若LME(4) 有广义自反解, 称问题Ⅰ相容; 否则, 称问题Ⅰ不相容.

记$w(Y)=\sum\limits_{i=1}^3(A_iYB_i+C_iY^TD_i)$, 那么LME(4) 的正规矩阵方程(指将LME(4) 按行拉直的线性代数方程组的正规方程组的矩阵形式)为

$ \sum\limits_{i=1}^{3}(A_i^Tw(Y)B_i^T+D_i(w(Y))^TC_i)=\sum\limits_{i=1}^{3}(A_i^TFB_i^T+D_iF^TC_i). $

(6)

记LME(4) 和LME(6) 在$Y=Y^{(k)}$处的残量依次为

$ R_k=F-w(Y^{(k)}), \widetilde{R}_k=\sum\limits_{i=1}^{3}(A_i^TR_kB_i^T+D_iR_k^TC_i). $

借鉴文献[7-9]的基本原理, 通过修改有关矩阵的类型或者算法, 建立求解问题Ⅰ的MCG算法(算法1) 如下.

第1步  给定初始矩阵$Y^{(1)}\in\Omega(P_1, P_2)$, 置$k:\, = 1$, 计算

$ R_k=F-w(Y^{(k)}), \quad\widetilde{R}_k=\sum\limits_{i=1}^{3}(A_i^TR_kB_i^T+D_iR_k^TC_i), \quad Z_k=\frac{1}{2}(\widetilde{R}_k+P_1\widetilde{R}_kP_2); $

第2步  若$R_k=O$, 或者$R_k\neq O$而$Z_k=O$, 停止; 否则, 计算

$ \alpha_k=\frac{\|R_{k}\|^2}{\|Z_k\|^2}, \quad Y^{(k+1)}=Y^{(k)}+\alpha_kZ_k; $

第3步  计算

$ \begin{eqnarray*}&R_{k+1}=F-w(Y^{(k+1)}),\quad\widetilde{R}_{k+1}=\sum\limits_{i=1}^{3}(A_i^TR_{k+1}B_i^T+D_iR_{k+1}^TC_i),\\ &\beta_{k+1}=-\frac{[\widetilde{R}_{k+1},Z_k]}{\|Z_k\|^2},\quad Z_{k+1}=\frac{1}{2}(\widetilde{R}_{k+1}+P_1\widetilde{R}_{k+1}P_2)+\beta_{k+1}Z_k;\end{eqnarray*} $

第4步  置$k:\, = k+1$, 转第2步.

可以验证, 算法1中的矩阵满足$Y^{(k)}, Z_k\in\Omega(P_1, P_2)$.对于算法1有以下收敛性定理(证明过程类似文献[7-8]).

定理1  设问题Ⅰ相容, 则对任意初始矩阵$Y^{(1)}\in\Omega(P_1, P_2)$, 算法1可在有限步迭代计算后得到问题Ⅰ的一个解; 若取初始矩阵满足

$ Y^{(1)}=\sum\limits_{i=1}^{3}(A_i^THB_i^T+D_iH^TC_i)+P_1(\sum\limits_{i=1}^{3}(A_i^THB_i^T+D_iH^TC_i))P_2, $

其中任意$H\in R^{n \times n}$, 则算法1可在有限步计算后得到问题Ⅰ的唯一极小范数解, 即LME(4) 的唯一极小范数广义自反解.

定理2  问题Ⅰ不相容的充要条件是存在正整数$k$, 使得由算法1得到的$R_k\neq O$, 而$Z_k=O$.

根据定理2, 在算法1中, 当$R_k\neq O$而$Z_k=O$时, 算法1中断.这表明LME(4) 没有广义自反解.因此, 需要求解问题Ⅱ, 即求LME(4) 的广义自反Ls解.

引理2^[10]  设$X\in R^{n \times n}$, 则存在$n^2$阶对称置换矩阵$T_{n, n}$, 使$\overline{{\rm vec}}(X^{{T}})=T_{n, n}\overline{{\rm vec}}(X)$.

下面通过构造等价的LME, 将求LME(4) 的广义自反Ls解问题, 转化为求等价的LME的广义自反解问题.然后参照算法1, 建立求LME(4) 的广义自反Ls解的迭代算法.约定矩阵的乘积运算优先于矩阵的Kronecker积运算, 引进记号:

$ \begin{array}{l}& y=\overline{{\hbox {vec}}}(Y),\quad f=\overline{{\hbox{vec}}}(F),\quad \widetilde{f}=[f^T,f^T]^T,\\ & N=\left[ {\begin{array}{*{20}{c}} {\sum\limits_{i = 1}^3 (A_i\bigotimes B_i^T+(C_i\bigotimes D_i^T)T_{n,n})} \\ {\sum\limits_{i = 1}^3 (A_iP_1\bigotimes B_i^TP_2+(C_iP_2\bigotimes D_i^TP_1)T_{n,n})} \\ \end{array}} \right],\\ & M_i=A_i^Tw(Y)B_i^T+D_i(w(Y))^TC_i,\quad Q_i=A_i^TFB_i^T+D_iF^TC_i,\\ & \widetilde{M}_i=A_i^Tw(P_1YP_2)B_i^T+D_i(w(P_1YP_2))^TC_i\quad(i=1,2,3),\\ & g(Y)=\sum\limits_{i=1}^3(M_i+P_1\widetilde{M}_iP_2),\quad Q=\sum\limits_{i=1}^3(Q_i+P_1Q_iP_2).\end{array} $

定理3  求解问题Ⅱ等价于求LME

$ g(Y)=Q $

(7)

的广义自反解, 且LME(7) 一定有广义自反解.

证  当$Y\in\Omega(P_1, P_2)$时, $P_1YP_2=Y$.因此, 求解问题Ⅱ等价于求$Y\in\Omega(P_1, P_2)$, 使得

$ \left\|\left[ {\begin{array}{*{20}{c}} {\sum\limits_{i = 1}^3 (A_iYB_i+C_iY^TD_i)} {\sum\limits_{i = 1}^3 (A_i(P_1YP_2) B_i+C_i(P_1YP_2)^T D_i)} \end{array}} \right]-\left[ {\begin{array}{*{20}{c}} {F} {F} \end{array}} \right]\right\|=\min. $

(8)

下面证明求解极小值问题(8) 等价于求解LME(7).将矩阵方程组

$ \left\{{\begin{array}{*{20}{l}} &{\sum\limits_{i = 1}^3 (A_iYB_i+C_iY^TD_i)}=F, \\ &{\sum\limits_{i = 1}^3 (A_i(P_1YP_2) B_i+C_i(P_1YP_2)^T D_i)}=F \end{array}}\right. $

(9)

按行拉直可得线性代数方程组$Ny=\widetilde{f}$, 求其Ls解等价于求LMEs(9) 的Ls解, 即求极小值问题(8) 的解.方程组$Ny=\widetilde{f}$的正规方程组为$N^TNy=N^T\widetilde{f}$, 将其还原为矩阵方程就是LME(7).因为求方程组$Ny=\widetilde{f}$的Ls解等价于求其正规方程组的解, 所以求解问题Ⅱ等价于求LME(7) 的广义自反解.

因为正规方程组$N^TNy=N^T\widetilde{f}$有解, 所以LME(7) 有解.设$\widetilde{Y}$是LME(7) 的一个解(未必是广义自反解), 那么$g(\widetilde{Y})=Q$.令$Y^*=\frac{1}{2}(\widetilde{Y}+P_1\widetilde{Y}P_2)$, 则$Y^*\in\Omega(P_1, P_2)$, 且有$g(Y^*)=Q$, 即LME(7) 有广义自反解.

参照算法1及文献[6], 建立求LME(7) 的广义自反解, 即求解问题Ⅱ的MCG算法(算法2) 如下.

第1步  给定初始矩阵$Y^{(1)}\in\Omega(P_1, P_2)$, 置$k:\, = 1$, 计算

$ R_k=Q-g(Y^{(k)}), \quad\widetilde{R}_k=g(R_k), \quad Z_k=\widetilde{R}_k; $

第2步  若$R_k=O$, 停止; 否则, 计算

$ \alpha_k=\frac{\|R_{k}\|^2}{\|Z_k\|^2}, \quad Y^{(k+1)}=Y^{(k)}+\alpha_kZ_k; $

第3步  计算

$ \begin{array}{l}& R_{k+1}=Q-g(Y^{(k+1)}),\quad\widetilde{R}_{k+1}=g(R_{k+1}),\\ & \beta_{k+1}=-\frac{[\widetilde{R}_{k+1},Z_k]}{\|Z_k\|^2},\quad Z_{k+1}=\widetilde{R}_{k+1}+\beta_{k+1}Z_k;\end{array} $

第4步  置$k:\, = k+1$, 转第2步.

可以验证, 算法2中的矩阵满足$Y^{(k)}, Z_k\in\Omega(P_1, P_2)$, 对于算法2有以下收敛性定理(证明过程类似文献[6]).

定理4   LME(7) 总是有广义自反解的, 对任意初始矩阵$Y^{(1)}\in\Omega(P_1, P_2)$, 算法2可在有限步计算后得到问题Ⅱ的一个解, 即LME(4) 的一个广义自反Ls解; 若取初始矩阵满足

$ Y^{(1)}=g(H)\quad(\mbox{任意}H\in\Omega(P_1, P_2)), $

则算法2可在有限步计算后得到LME(7) 的唯一极小范数广义自反解, 即LME(4) 的唯一极小范数广义自反Ls解.

求方程(1) 的广义自反解, 可采用以下两种计算方案.

方案一

第1步  给定初始矩阵$X^{(1)}\in\Omega(P_1, P_2)$, 置$k:\, = 1$;

第2步  如果$\psi(X^{(k)})=O$, 停止; 否则, 采用算法1求$Y^{(k)}\in\Omega(P_1, P_2)$, 使满足LME(3);当算法1中断(表明LME(3) 无广义自反解)时, 采用算法2求$Y^{(k)}\in\Omega(P_1, P_2)$, 使得

$ \|\phi_{X^{(k)}}(Y^{(k)})+\psi(X^{(k)})\|=\min; $

(10)

第3步  计算$X^{(k+1)}=X^{(k)}+Y^{(k)}$, 置$k:\, = k+1$, 转第2步.

方案二

第1步  给定初始矩阵$X^{(1)}\in\Omega(P_1, P_2)$, 置$k:\, = 1$;

第2步  如果$\psi(X^{(k)})=O$, 停止; 否则, 采用算法2求$Y^{(k)}\in\Omega(P_1, P_2)$, 使满足式(10).当$\phi_{X^{(k)}}(Y^{(k)})=-\psi(X^{(k)})$有广义自反解时, 它的广义自反Ls解就是它的广义自反解;

第3步  计算$X^{(k+1)}=X^{(k)}+Y^{(k)}$, 置$k:\, = k+1$, 转第2步.

例1  采用两种计算方案求方程(1) 的广义自反解.设对称正交矩阵$P_1$为次单位矩阵, $P_2=I-2u^Tu$, 其中$u=(0, 0, 1)$, 系数矩阵和右端项如下:

$ \begin{eqnarray*}& D=\left[ {\begin{array}{*{20}{c}} 1&1&0 \\ 0&1&1 \\ 1&0&-1 \\ \end{array}} \right],\quad E_5=\left[ {\begin{array}{*{20}{c}} -12&-12&4 \\ -12&-12&4 \\ -12&-12&-4 \\ \end{array}} \right],\quad u_1=\left[ {\begin{array}{*{20}{c}} 1 \\ 1 \\ 0 \\ \end{array}} \right],\quad u_2=\left[ {\begin{array}{*{20}{c}} 0 \\ 1 \\ 1 \\ \end{array}} \right],\\ & A=D^T,\quad B=C=I,\quad E_1=E_2=-u_2u_2^T,\quad E_3=-u_1u_1^T,\quad E_4=u_1u_2^T,\end{eqnarray*} $

选取初始矩阵$X^{(1)}=I+P_1IP_2\in\Omega(P_1, P_2)$, $Y^{(1)}=O\in\Omega(P_1, P_2)$, 终止准则$\varepsilon=10^{-9}$, 两种方案的计算结果均为

$ X^{(6)}=\left[ {\begin{array}{*{20}{c}} 2.0000&2.0000&-0.0000 \\ 2.0000&2.0000& 0 \\ 2.0000&2.0000& 0.0000 \\ \end{array}} \right]. $

方案一中, 牛顿算法迭代6次, 算法1和算法2分别迭代36次和35次; 方案二中, 牛顿算法迭代6次, 算法2迭代41次.

计算结果表明:方程(1) 有广义自反解, 由方案一或方案二得到的$X^{(6)}$是它的一个广义自反解; 在方案一中, 采用算法1求解LME(3) 时, 会遇到算法1中断的情况, 此时LME(3) 没有广义自反解$Y^{(k)}$, 需要采用算法2求LME(3) 的广义自反Ls解$Y^{(k)}$.

例2   采用两种计算方案求方程(1) 的广义自反解.取$u=(1, 0, \cdots, 0)$, 对称正交矩阵$P_1=I-2u^Tu$, $P_2$为以次单位矩阵为子矩阵的块对角矩阵, 而

$ \begin{array}{l}& \widetilde{B}_i=\left[ {\begin{array}{*{20}{c}} 5&1&0 \\ 0&6&1 \\ 1&0&-7 \\ \end{array}} \right],\quad \widetilde{D}_i=\left[ {\begin{array}{*{20}{c}} 1&1&0 \\ 0&1&1 \\ 1&0&-1 \\ \end{array}} \right],\quad \widetilde{X}_i=\left[ {\begin{array}{*{20}{c}} 1&0&1 \\ &1&0 \\ & &1 \\ \end{array}} \right],\\ & B=[O,\widetilde{B}_i,I]_{i=1}^N,\quad D=[I,\widetilde{D}_i,O]_{i=1}^N\quad\mbox{(块三对角矩阵-下同), }\\ & A=D^T,\quad C=B^T,\quad E_1=E_4=[I,2I,I]_{i=1}^N,\quad E_2=E_3=[I,{\hbox{ones}}(3),I]_{i=1}^N,\\ & \widehat{X}={\hbox{diag}}(\widetilde{X}_i,\cdots,\widetilde{X}_i)+P_1{\hbox{diag}}(\widetilde{X}_i,\cdots,\widetilde{X}_i)P_2,\\ & E_5=A\widehat{X}B+C\widehat{X}^T D+\widehat{X}E_1\widehat{X}+\widehat{X}E_2\widehat{X}^T+ \widehat{X}^TE_3\widehat{X}+\widehat{X}^TE_4\widehat{X}^T. \end{array} $

选取牛顿算法的初始矩阵

$ X^{(1)}=I+P_1IP_2\in\Omega(P_1, P_2), $

MCG算法的初始矩阵$Y^{(1)}=O\in\Omega(P_1, P_2)$, 终止准则$\varepsilon=10^{-9}$, 两种方案的计算结果见表 1, 表中的外迭代指牛顿算法, 内迭代指MCG算法, 计算时间单位为$s$.

计算结果表明:当LME(3) 总是有广义自反解$Y^{(k)}$时, 方案一比方案二的效率高.

运用牛顿算法求Riccati矩阵方程的广义自反解, 并采用MCG算法求导出的LME的广义自反解或者广义自反最小二乘解, 建立了求Riccati矩阵方程的广义自反解的双迭代算法, 数值算例表明双迭代算法是有效的.修改算法中初始矩阵的类型, 算法1中涉及的矩阵$Z_1$与$Z_{k+1}$的计算公式, 以及算法2中涉及的LME(7) 的构造方式, 还可建立求Riccati矩阵方程的其它特殊解的迭代算法.