本文所用记号: ${\Bbb R}^{m\times n}$为所有$m\times n$实矩阵的集合. $I_k$表示$k$阶单位矩阵, $J_k$表示$k$阶反序单位矩阵. $M^T$和${{M}^{\dagger }}$分别表示矩阵$A$的转置和Moore-Penrose广义逆.矩阵不等式$M\geq N$表示矩阵$M-N$为非负矩阵. $[\alpha]$表示不大于$\alpha$的最大整数.对于矩阵$M=(a_{ij})\in {\Bbb R}^{m\times n}$, $M_+$表示第$ij$位置的元素等于$\max \{0, a_{ij}\}$的矩阵, $M_{-}$表示第$ij$位置的元素等于$\min \{0, a_{ij}\}$的矩阵.在矩阵空间${\Bbb R}^{m\times n}$中定义内积为$\langle M, N\rangle={\rm trace}(N^TM)$, 由此导出的矩阵范数为Frobenius范数, 记为$\|\cdot\|.$

矩阵方程最小二乘问题在结构设计, 振动理论的逆问题, 线性规划等领域有广泛的应用.近年来众多学者更多地关注于应用各类数值算法求解矩阵不等式约束下的矩阵方程最小二乘问题.如Bouhamidi^[1]等利用谱投影梯度算法^[2]研究边界约束$ L\leq X\leq U$下广义Sylvester矩阵方程

$ \begin{equation} \sum\limits_{i=1}^{q}A_iXB_i=C \end{equation} $

(1.1)

最小二乘解问题, 其中$L, U$为给定矩阵; Escalante和Raydan^[3]等利用Dysktra交替投影算法研究对称、半正定、边界约束下矩阵方程$AX=B$最小二乘问题, 文献[4]进而将Dysktra交替投影算法推广应用于双对称、$\varepsilon$正定、非负约束矩阵方程$AX=B$最小二乘问题; 文献[5]应用松弛交替投影算法研究非负约束下的相容矩阵方程$AXB=C$求解问题.但需指出的是, 上述三类求解算法都属于投影类算法, 其缺点是要求容易计算所涉及的不等式约束集合上的投影, 故该类算法只能应用于简单矩阵不等式约束条件.文献[6-7]利用极分解理论分别研究了矩阵不等式约束$CXD\geq E$约束下矩阵方程$AX=B$的一般解和双对称解, 但该方法要求先求出矩阵方程的通解表达式, 对于复杂矩阵方程的最小二乘问题, 由于其通解表达式难以确定, 故极分解方法失效.为了能避免复杂投影矩阵的计算, 同时又能方便地推广到一般的矩阵方程最小二乘问题, 本文将优化理论中求凸集上光滑函数最小值的增广Lagrangian方法^[8-9]应用于求解如下矩阵不等式$EXF\geq D$约束下广义Sylvester矩阵方程(1.1) 最小二乘问题, 数值算例说明该方法是非常有效的.

问题1.1 给定正整数$l, h, m, n, p, t, q$, 线性约束矩阵集合$\Omega$, 矩阵$A_i\in {\Bbb R}^{l\times m}$, $B_i\in {\Bbb R}^{n\times h}$, $(i=1, 2, \cdots, q)$, $C\in {\Bbb R}^{l\times h}$, $E\in {\Bbb R}^{p\times m}$, $F\in {\Bbb R}^{n\times t}$和$D\in {\Bbb R}^{p\times t}$, 求$X\in {\Bbb R}^{m\times n}$, 使得

$ \begin{array}{l} {\rm{minimize}}F(X) = \frac{1}{2}{\left\| {\sum\limits_{i = 1}^q {{A_i}} X{B_i} - C} \right\|^2}\\ {\rm{subject}}\;{\rm{to}}X \in \Omega ,\;\;EXF \ge D. \end{array} $

(1.2)

本文讨论的线性约束矩阵集合$\Omega$为实$(R, S)$对称矩阵集合.实$(R, S)$对称矩阵^[10-11]是一系列线性约束结构矩阵的一般形式, 其定义如下.

定义1.1 若$R\in {\Bbb R}^{m\times m}$, $S\in {\Bbb R}^{n\times n}$且$R={R}^{-1}\neq \pm I_m$, $S=S^{-1}\neq \pm I_n$, 则称$R$和$S$为非平凡实对合矩阵.对于矩阵$X\in {\Bbb R}^{m\times n}$, 若$RXS=X$, 则称$X$为实$(R, S)$对称矩阵.记${\cal S}_{R, S}^{m\times n}$表示全体$m\times n$阶实$(R, S)$对称矩阵构成的集合.

由实$(R, S)$对称矩阵的定义知下列结论成立.

推论1.1 若$N\in {\Bbb R}^{m\times m}$, 则$N+RNS\in {\cal S}_{R, S}^{m\times n}$.若$M\in {\cal S}_{R, S}^{m\times n}$, 则

$ \langle M, N\rangle=\langle M, \frac{1}{2}(N+RNS)\rangle. $

注意到, 若取特定的实对合矩阵$R$和$S$, 由$(R, S)$对称矩阵可得到几类特定的线性约束结构矩阵.如若$R=J_m$, $S=J_n$, 其中$J_s$为$s$阶反序单位矩阵, 则$(R, S)$对称矩阵即为长方形中心对称矩阵, 若还有$m=n$, 则$(R, S)$对称矩阵退化为一般中心对称矩阵; 若$R=S$为反轮换矩阵, 即$R^2=S^2=-I_n$, 则$(R, S)$对称矩阵即为伪中心对称矩阵; 若$R$和$S$为广义反射矩阵, 即$R^H = R$, $R^2 = I$, $S^H = S$, $S^2 = I$, 则$(R, S)$对称矩阵即为自反矩阵:若$R = S$为$(k, p)$镜像矩阵, 则$(R, S)$对称矩阵即为$(k, p)$镜像对称矩阵.

事实上, 若去掉约束条件$X\in \Omega$, 通过拉直算子, 问题(1.2) 可直接转换为如下线性不等式约束二次规划问题

$ \begin{array}{l} {\rm{minimize}}\;\;f(x) = \frac{1}{2}{x^T}Mx + {g^T}x + c,\\ {\rm{subject}}\;{\rm{to}}\;\;Gx \ge f, \end{array} $

进而采用凸规划问题中常用的内点法、积极集法、信赖域法、牛顿法等求解.但是采用拉直算子会成倍数地增加问题的规模, 进而导致计算中大幅度地增加存储量和计算复杂度.易知若问题(1.2) 中$q=1$, $m=n=t=s=100$, 则通过拉直之后线性方程和线性不等式的系数矩阵的阶数均为$10000\times 10000$.

为方便讨论, 定义如下线性算子

$ \begin{eqnarray*} \begin{array}{rrll} {\cal A}:& {\Bbb R}^{m\times n}&\longrightarrow& {\Bbb R}^{l\times h}, \\ &X&\longrightarrow& \sum\limits_{i=1}^{q}A_iXB_i. \end{array} \end{eqnarray*} $

由内积定义$\langle M, N\rangle={\rm trace}(N^TM)$知, $ {\cal A}$关于此内积的转置算子${\cal A}^T$定义为

$ \begin{eqnarray*} \begin{array}{rrll} {\cal A}^T:& {\Bbb R}^{l\times h}&\longrightarrow& {\Bbb R}^{m\times n}, \\ &M&\longrightarrow& \sum\limits_{i=1}^{q}A_i^TMB_i^T. \end{array} \end{eqnarray*} $

先给出实$(R, S)$对称矩阵的结构特征.对于实对合矩阵$R$, 必存在矩阵$P_R\in {\Bbb R}^{m\times r}$和$Q_R\in {\Bbb R}^{m\times (m-r)}$使得

$ \begin{equation} P_R^TP_R=I_r, \ \ Q_R^TQ_R=I_{m-r}\ \ \mbox{且}\ \ RP_R=P_R, \ \ RQ_R=Q_R. \end{equation} $

(2.1)

矩阵$P_R$和矩阵$Q_R$的列分别构成矩阵$R$的属于特征值1和$-1$所对应的特征空间的正交基.对于实对合矩阵$S$, 必存在矩阵$U_S\in {\Bbb R}^{n\times s}$和$V_S\in {\Bbb R}^{n\times (n-s)}$使得

$ \begin{equation} U_S^TU_S=I_s, \ \ V_S^TV_S=I_{n-s}\ \ \mbox{且}\ \ RU_S=U_S, \ \ RV_S=V_S. \end{equation} $

(2.2)

矩阵$U_S$和矩阵$V_S$的列分别构成矩阵$S$的属于特征值1和$-1$所对应的特征空间的正交基.尽管$P_R$, $Q_R$, $U_S$和$V_S$不是唯一的, 但$P_R$和$Q_R$可以分别通过对$I_m+R$和$I_m-R$实施施密特正交化法得到, 而$U_S$和$V_S$可以分别通过对$I_n+S$和$I_n-S$实施施密特正交化法得到.如果$R$和$S$为置换矩阵, 则只需很少的计算量即可得到.特别地, 若$R=J_{2k}$, 即$2k\times 2k$阶的反序单位矩阵, 则可取

$ \begin{equation} \begin{array}{cc} P_R=\frac{1}{\sqrt{2}}\left[\begin{array}{cc} I_k\\J_k \end{array}\right], & Q_R=\frac{1}{\sqrt{2}}\left[\begin{array}{cc} I_k\\-J_k \end{array}\right]. \end{array} \end{equation} $

(2.3)

若$R=J_{2k+1}$, 则可取

$ \begin{equation}\begin{array}{cc} P_R=\frac{1}{\sqrt{2}}\left[\begin{array}{cc} I_k~~&0_{k, 1}\\ 0_{1, k}~~&\sqrt{2}\\ J_k~~&0_{k, 1} \end{array}\right], & Q_R=\frac{1}{\sqrt{2}}\left[\begin{array}{cc} I_k\\ 0_{1, k} \\-J_k \end{array}\right]. \end{array} \end{equation} $

(2.4)

若定义

$ \widehat{P}_R=\frac{P_R^T(I+R)}{2}, ~\widehat{Q}_R=\frac{Q_R^T(I-R)}{2}, ~ \widehat{U}_S=\frac{U_S^T(I+S)}{2}, ~\widehat{V}_S=\frac{V_S^T(I-S)}{2}. $

则有

$ \begin{array}{cccc} {\widehat{P}_R}P_R=I_r, & {\widehat{P}_R}Q_R=0, & \widehat{Q}_RQ_R=I_r, & \widehat{Q}_R{P}_R=0;\\ \widehat{U}_S{U}_S=I_s, & \widehat{U}_S{V}_S=0, & \widehat{V}_S{V}_S=I_s, & \widehat{V}_S{U}_S=0.\\ \end{array} $

由上所列方程即得

$ \begin{array}{cc} \left[\begin{array}{cc} P_R~&Q_R \end{array}\right]^{-1}=\left[\begin{array}{cc} \widehat{P}_R\\ \widehat{Q}_R \end{array}\right], & \left[\begin{array}{cc} U_S&~V_S \end{array}\right]^{-1}=\left[\begin{array}{cc} \widehat{U}_S\\ \widehat{V}_S \end{array}\right]. \end{array} $

若实对合矩阵$R$和$S$是正规的, 则有$\widehat{P}_R={P}_R^T$, $\widehat{Q}_R={Q}_R^T$, $\widehat{U}_S={U}_S^T$和$\widehat{V}_S={V}_S^T$.这时有

$ \begin{array}{cc} \left[\begin{array}{cc} P_R&~Q_R \end{array}\right]^{-1}=\left[\begin{array}{cc} {P}_R^T\\ {Q}_R^T \end{array}\right], & \left[\begin{array}{cc} U_S~&V_S \end{array}\right]^{-1}=\left[\begin{array}{cc} {U}_S^T\\ {V}_S^T \end{array}\right]. \end{array} $

引理2.1^[10] 矩阵$X\in {\cal S}_{R, S}^{m\times n}$的充要条件为$X$可表示为

$ \begin{equation}\label{lemma-001} X=\left[\begin{array}{cc} P_R~&Q_R \end{array}\right]\left[\begin{array}{cc}X_1~&0\\0~&X_2\end{array}\right]\left[\begin{array}{cc} \widehat{U}_S\\ \widehat{V}_S \end{array}\right], \ \ X_1\in {\Bbb R}^{r\times s}, \ X_2\in {\Bbb R}^{(m-r)\times (n-s)}. \end{equation} $

(2.5)

引理2.2 对任意$X\in {\cal S}_{R, S}^{m\times n}$, 若其分解如(2.5) 式所示, 则存在$(R, S)$对称拉直算子矩阵$C_p\in {\Bbb R}^{mn\times (rs+(m-r)(n-s))}$使得

$ \begin{eqnarray*} vec(X)=C_pvec(X_1, X_2). \end{eqnarray*} $

证 由(2.5) 式知

$ \begin{eqnarray*} vec(X)=\left(\left[\begin{array}{cc} \widehat{U}_S\\ \widehat{V}_S \end{array}\right]^T\bigotimes \left[\begin{array}{cc} P_R~&Q_R \end{array}\right]\right)vec\left(\left[\begin{array}{cc}X_1~&0\\0~&X_2\end{array}\right]\right). \end{eqnarray*} $

而

$ \begin{eqnarray*}vec\left(\left[\begin{array}{cc} X_1\\ 0_{m-r, s} \end{array}\right]\right)=diag\left(\left[\begin{array}{cc} I_r\\ 0_{m-r, r} \end{array}\right], \cdots, \left[\begin{array}{cc} I_r\\ 0_{m-r, r} \end{array}\right]\right)vec(X_1), \end{eqnarray*} $

$ \begin{eqnarray*}vec\left(\left[\begin{array}{cc} 0_{r, n-s}\\ X_2 \end{array}\right]\right)=diag\left(\left[\begin{array}{cc} 0_{r, m-r}\\ I_{m-r} \end{array}\right], \cdots, \left[\begin{array}{cc} 0_{r, m-r}\\ I_{m-r} \end{array}\right]\right)vec(X_2).\end{eqnarray*} $

为了不引起混淆, 在上述两式中标明各个零矩阵的阶数.上述第一个等式右边的分块对角矩阵总共有$s$个对角块, 第二个等式右边的分块对角矩阵总共有$n-s$个对角块.令

$ \begin{eqnarray*} W_1=diag\left(\left[\begin{array}{cc} I_r\\ 0_{m-r, r} \end{array}\right], \cdots, \left[\begin{array}{cc} I_r\\ 0_{m-r, r} \end{array}\right]\right)\in {\Bbb R}^{ms\times rs}, \end{eqnarray*} $

$ \begin{eqnarray*} W_2=diag\left(\left[\begin{array}{cc} 0_{r, m-r}\\ I_{m-r} \end{array}\right], \cdots, \left[\begin{array}{cc} 0_{r, m-r}\\ I_{m-r} \end{array}\right]\right)\in {\Bbb R}^{m(n-s)\times (m-r)(n-s)}, \end{eqnarray*} $

则有

$ \begin{eqnarray*} vec\left(\left[\begin{array}{cc}X_1~~&0\\0&X_2\end{array}\right]\right)= \left[\begin{array}{cc}W_1~~&0\\0~~&W_2\end{array}\right] \left[\begin{array}{cc} vec(X_1) \\ vec(X_2) \end{array}\right]. \end{eqnarray*} $

令

$ \begin{equation}\label{C_p} C_p=\left(\left[\begin{array}{cc} \widehat{U}_S\\ \widehat{V}_S \end{array}\right]^T\bigotimes \left[\begin{array}{cc} P_R~~&Q_R \end{array}\right]\right)\left[\begin{array}{cc}W_1~~&0\\0~~&W_2\end{array}\right]\in {\Bbb R}^{mn\times (rs+(m-r)(n-s))}, \end{equation} $

(2.6)

则有

$ vec(X)=C_p\left[\begin{array}{cc} vec(X_1) \\ vec(X_2) \end{array}\right]=C_pvec(X_1, X_2). $

证毕.

引理2.3 令$W\in {\Bbb R}^{m\times n}$, 记

$ \widehat{W}=\left[\begin{array}{cc} \widehat{P}_R\\ \widehat{Q}_R \end{array}\right]W\left[\begin{array}{cc} U_S~~&V_S \end{array}\right]= \left[\begin{array}{cc} \widehat{W}_{11}~~&\widehat{W}_{12}\\ \widehat{W}_{21}~~&\widehat{W}_{22} \end{array}\right], $

则有

$ C_p^\dagger vec(W)=vec(\widehat{W}_{11}, \widehat{W}_{22}). $

证 由$C_p$的定义(2.6) 知

$ \begin{eqnarray*} \begin{array}{rl} C_p^\dagger vec(W)&=\left[\begin{array}{cc}W_1^\dagger~~&0\\ 0~~&W_2^\dagger\end{array}\right]\left(\left[\begin{array}{cc} \widehat{U}_S\\ \widehat{V}_S \end{array}\right]^T\bigotimes \left[\begin{array}{cc} P_R~~&Q_R \end{array}\right]\right)^\dagger vec(W)\\ &=\left[\begin{array}{cc}W_1^\dagger~~&0\\ 0~~&W_2^\dagger\end{array}\right] \left(\left[\begin{array}{cc} U_S~~&V_S \end{array}\right]^T\bigotimes \left[\begin{array}{cc} \widehat{P}_R\\ \widehat{Q}_R \end{array}\right] \right)vec(W)\\ &=\left[\begin{array}{cc}W_1^\dagger~~&0\\ 0~~&W_2^\dagger\end{array}\right] vec\left(\left[\begin{array}{cc} \widehat{P}_R\\ \widehat{Q}_R \end{array}\right]W\left[\begin{array}{cc} U_S ~~&V_S \end{array}\right]\right)\\ &=\left[\begin{array}{cc}W_1^\dagger~~&0\\ 0~~&W_2^\dagger\end{array}\right]vec(\widehat{W}). \end{array} \end{eqnarray*} $

由$W_1$和$W_2$的定义知

$ \left[\begin{array}{cc} W_1^\dagger~~& 0 \end{array}\right]vec(\widehat{W})=vec(\widehat{W}_{11}), \ \ \left[\begin{array}{cc} 0~~& W_2^\dagger \end{array}\right]vec(\widehat{W})=vec(\widehat{W}_{22}), $

故

$ C_p^\dagger vec(W)=\left[\begin{array}{cc} vec(\widehat{W}_{11})\\ vec(\widehat{W}_{22}) \end{array}\right]=vec(\widehat{W}_{11}, \widehat{W}_{22}). $

证毕.

引理2.4 对任意矩阵$W\in {\Bbb R}^{m\times n}$, 若$C_p^Tvec(W)=0$, 则有$W+RWS=0$.

证 由$C_p$的定义(2.6) 知, $C_p$为列满秩矩阵, 有$C_p^\dagger=(C_p^TC_p)^{-1}C_p^T$, 故$C_p^T=(C_p^TC_p)C_p^\dagger$.而$C_p^TC_p$为可逆矩阵, 故

$ C_p^T vec(W)=0\Longleftrightarrow C_p^\dagger vec(W)=0. $

若$W$的分块结构为

$ W= \left[\begin{array}{cc} {P}_R~& {Q}_R \end{array}\right]\left[\begin{array}{cc} \widehat{W}_{11}~&\widehat{W}_{12}\\ \widehat{W}_{21}~&\widehat{W}_{22} \end{array}\right]\left[\begin{array}{cc} \widehat{U}_S\\ \widehat{V}_S \end{array}\right], $

由引理2.3的结论, 若$C_p^\dagger vec(W)=0$, 则有$vec(\widehat{W}_{11})=0$和$vec(\widehat{W}_{22})=0$, 即$\widehat{W}_{11}=0$, $\widehat{W}_{22}=0$, 故

$ \begin{eqnarray*} W+RWS&=&\left[\begin{array}{cc} {P}_R~& {Q}_R \end{array}\right]\left[\begin{array}{cc} 0~&\widehat{W}_{12}\\ \widehat{W}_{21}~&0 \end{array}\right]\left[\begin{array}{cc} \widehat{U}_S\\ \widehat{V}_S \end{array}\right]+R\left[\begin{array}{cc} {P}_R~& {Q}_R \end{array}\right]\left[\begin{array}{cc} 0~&\widehat{W}_{12}\\ \widehat{W}_{21}~&0 \end{array}\right]\left[\begin{array}{cc} \widehat{U}_S\\ \widehat{V}_S \end{array}\right]S\\ &=&\left[\begin{array}{cc} {P}_R~& {Q}_R \end{array}\right]\left[\begin{array}{cc} 0~&\widehat{W}_{12}\\ \widehat{W}_{21}~&0 \end{array}\right]\left[\begin{array}{cc} \widehat{U}_S\\ \widehat{V}_S \end{array}\right]+\left[\begin{array}{cc} {P}_R~&-{Q}_R \end{array}\right]\left[\begin{array}{cc} 0&\widehat{W}_{12}\\ \widehat{W}_{21}~&0 \end{array}\right]\left[\begin{array}{cc} \widehat{U}_S\\-\widehat{V} \end{array}\right]S. \end{eqnarray*} $

证毕.

引理3.1^[12]  假设$x^*$是二次规划问题

$ \min\limits_{x\in R^s}f(x)=\frac{1}{2}x^TMx+g^Tx+c\ \ {\rm subject}\ \ {\rm to}\ \ Gx\geq b $

的一个局部极小点, 则存在向量$y^*$使得

$ Mx^*+g-G^Ty^*=0, \ \ Gx^*\geq b, \ \ \langle y^*, Gx^*-b\rangle=0, \ \ y^*\geq 0. $

定理3.1 矩阵$X^*\in {\cal S}_{R, S}^{m\times n}$是问题(1.2) 的解的充分必要条件是存在非负矩阵$Y^*\in {\Bbb R}^{p\times q}$使得下列条件满足

$ \begin{equation}\label{theorem-shi-1} \left\{\begin{array}{l} {\cal A}^T\left({\cal A}(X^*)-C\right)+R{\cal A}^T\left({\cal A}(X^*)-C\right)S-E^TY^*F^T-RE^TY^*F^TS=0, \\ EX^*F-D\geq 0, \\ \langle Y^*, EX^*F-D\rangle=0. \end{array} \right. \end{equation} $

(3.1)

证 假设存在非负矩阵$Y^*$使得条件(3.1) 满足.定义矩阵函数

$ \widetilde{{\cal F}}(X)={\cal F}(X)+\langle Y^*, D-EXF\rangle, $

则对任意$\widetilde{W}\in {\cal S}_{R, S}^{m\times n}$有

$ \begin{eqnarray*} \widetilde{{\cal F}}(X^*+\widetilde{W}) &=&\frac{1}{2}\|{\cal A}(X^*+\widetilde{W})-C\|^2+\langle Y^*, D-E(X^*+\widetilde{W})F\rangle\\ &=&\frac{1}{2}\|{\cal A}(X^*)-C\|^2+\frac{1}{2}\| {\cal A}(\widetilde{W})\|^2\\ &&+\langle {\cal A}(\widetilde{W}), {\cal A}(X^*)-C\rangle+\langle Y^*, EX^*F\rangle -\langle Y^*, E\widetilde{W}F\rangle\\ &=&\widetilde{{\cal F}}(X^*)+\frac{1}{2}\| {\cal A}(\widetilde{W})\|^2+\langle \widetilde{W}, {\cal A}^T({\cal A}(X^*)-C)-E^TY^*F^T\rangle\\ &=&\widetilde{{\cal F}}(X^*)+\frac{1}{2}\| {\cal A}(\widetilde{W})\|^2+\frac{1}{2}\langle \widetilde{W}, {\cal A}^T({\cal A}(X^*)-C)\\ &&+R{\cal A}^T({\cal A}(X^*)-C)S-E^TY^*F^T-RE^TY^*F^TS\rangle\\ &=&\widetilde{{\cal F}}(X^*)+\frac{1}{2}\| {\cal A}(\widetilde{W})\|^2\\ &\geq& \widetilde{{\cal F}}(X^*). \end{eqnarray*} $

上述不等式说明$X^*$是矩阵函数$\widetilde{ {\cal F}}(X)$的全局最小解。因为$\langle Y^*, EX^*F-D\rangle=0$和$\widetilde{{\cal F}}(X)\geq \widetilde{{\cal F}}(X^*)$对任意$X\in {\cal S}_{R, S}^{m\times n}$都成立, 故

$ \begin{eqnarray*} {\cal F}(X)&\geq& {\cal F}(X^*)+\langle Y^*, D-EX^*F\rangle-\langle Y^*, D-EXF\rangle\\ &=&{\cal F}(X^*)-\langle Y^*, D-EXF\rangle. \end{eqnarray*} $

注意到$Y^*\geq 0$, 故对任意满足矩阵不等式$EXF-D\geq 0$的所有$X\in {\cal S}_{R, S}^{m\times n}$都有

$ {\cal F}(X)\geq {\cal F}(X^*), $

即$X^*$是问题(1.2) 的解.

反过来, 假设$X^*$是问题(1.2) 的解, 下面说明存在非负矩阵$Y^*$使得条件(3.1) 成立.事实上, 由拉直映射, 问题(1.2) 等价于

$ \begin{eqnarray} \rm{minimize}&&{\cal F}(X)= \frac{1}{2}vec(\widetilde{X})^TC_p^T\sum\limits_{i=1}^{q}\left(B_i\otimes A_i^T\right)\sum\limits_{i=1}^{q}\left(B_i^T\otimes A_i\right)C_pvec(\widetilde{X})\\ && -vec(C)^T(B^T\otimes A)C_pvec(\widetilde{X})+\frac{1}{2}vec(C)^Tvec(C)\\ \rm{subject}\ \rm{to} && (F^T\otimes E)C_pvec(\widetilde{X})\geq vec(D), \end{eqnarray} $

(3.2)

其中$C_p$为第2节中讨论的$(R, S)$对称拉直算子矩阵, 为方便讨论, 记$\widetilde{X}$为引理2.2中与$X\in {\cal S}_{R, S}^{m\times n}$所对应的矩阵$(X_1, X_2)$.因此, 若$X^*$是问题1.2的解, 即$\widetilde{X}^*=(X_1^*, X_2^*)$是问题3.2的解, 由引理3.1, 存在非负矩阵$Y^*\in {\Bbb R}^{r\times t}$使得下列条件满足

$ \begin{eqnarray*} \left\{\begin{array}{l} C_p^T\bigg[ \sum\limits_{i=1}^{q}\left(B_i\otimes A_i^T\right)\sum\limits_{i=1}^{q}\left(B_i^T\otimes A_i\right)C_pvec(\widetilde{X}^*)\\[4mm] -\sum\limits_{i=1}^{q}\left(B_i\otimes A_i^T\right)vec(C)-(F\otimes E^T)vec(Y^*)\bigg] =0, \\[3mm] (F^T\otimes E)C_pvec(\widetilde{X}^*)\geq vec(D), \\ \langle vec(Y^*), (F^T\otimes E)C_pvec(\widetilde{X}^*)- vec(D)\rangle=0. \end{array} \right. \end{eqnarray*} $

由$vec(X^*)=C_pvec(\widetilde{X}^*)$及引理2.4知, 上述条件等价于

$ \begin{eqnarray*} \left\{\begin{array}{l} \sum\limits_{i=1}^{q}A_i^T\left(\sum\limits_{i=1}^{q}A_iX^*B_i-C\right)B_i^T +R\sum\limits_{i=1}^{q}A_i^T\left(\sum\limits_{i=1}^{q}A_iX^*B_i-C\right)B_i^TS\\[3mm] -E^TY^*F^T-RE^TY^*F^TS=0, \\ EX^*F-D\geq 0, \\ \langle Y^*, EX^*F-D\rangle=0, \end{array} \right. \end{eqnarray*} $

即(3.1) 式成立, 定理得证.

定义增广Lagrangian函数

$ \begin{equation}\label{lagrangian-1} L_\rho(X, Z)=\frac{1}{2}\left\| {\cal A}(X)-C\right\|^2+\frac{\rho}{2}\left\|\left(D-EXF+\frac{Z}{\rho}\right)_+\right\|^2, \end{equation} $

(3.3)

其中$Z\geq 0$称为Lagrangian乘子矩阵, $\rho>0$为罚函数.易知

$ \nabla_X L_\rho(X, Z)={\cal A}^T\left( {\cal A}(X)-C\right)-\rho E^T\left(D-EXF+\frac{Z}{\rho}\right)_+F^T. $

我们提出求解矩阵最优化问题(1.2) 的迭代格式如下.

算法3.1  (增广Lagrangian方法求解问题(1.2))

第1步:输入矩阵$A, B, C, E, F, D$和充分大的参数矩阵$Z_{\max}>0$.输入标量$\gamma>1$, $r\in (0, 1)$, $\rho_1>0$, 误差容许值$\varepsilon>0$和近似误差容许值$\varepsilon_k\downarrow 0$.选取初始矩阵$\overline{Z}_1$使得$0\leq \overline{Z}_1\leq Z_{\max}$.令$k\leftarrow 1$.

第2步:计算如下问题的近似估计矩阵$X_k$

$ \begin{equation} \begin{array}{rl} \rm{minimize}& L_{\rho_k}(X, \overline{Z}_k)\ \ {\rm subject}\ {\rm to}\ \ X\in {\cal S}_{R, S}^{m\times n}. \end{array} \end{equation} $

(3.4)

即计算近似矩阵$X_k$使得$\|\nabla L_{\rho_k}(X_k, \overline{Z}_k)+R\nabla L_{\rho_k}(X_k, \overline{Z}_k)S\| < \varepsilon_k$.

第3步:定义$Z_k=(\overline{Z}_k+\rho_k(D-EX_kF))_+.$

第4步:若$\|\nabla L_{\rho_k}(X_k, \overline{Z}_k)+R\nabla L_{\rho_k}(X_k, \overline{Z}_k)S\| < \varepsilon$, $\|\{EX_kF-D, Z_k\}_{-}\| < \varepsilon$, 则终止运行, 输出结果.否则, 转第5步.

第5步:如果$k=1$或

$ \begin{equation} \|\{EX_kF-D, Z_k\}_{-}\|\leq r\|\{EX_{k-1}F-D, Z_{k-1}\}_{-}\|, \end{equation} $

(3.5)

定义$\rho_{k+1}=\rho_k$.否则, 定义$\rho_{k+1}=\gamma\rho_k.$

第6步:按如下方法更新$\overline{Z}_{k+1}$使得$0\leq \overline{Z}_{k+1}\leq Z_{max}:$当$0\leq ({Z}_k)_{ij}\leq (Z_{max})_{ij}$, $i=1, 2, \cdots, p$, $j=1, 2, \cdots, q$时, 令$(\overline{Z}_{k+1})_{ij}=(Z_k)_{ij}$.

第7步:令$k\leftarrow k+1$并转第2步.

在运行算法3.1时, 选取参数$\gamma=5$, $r=0.5$, $\rho_1=1$; 选取矩阵$Z_{\max}$的所有元素等于$10^{20}$, 选取初始矩阵$\overline{Z}_1$为零矩阵以及误差容许值$\varepsilon=10^{-8}$, 迭代精度$\varepsilon_k$按如下方式选取

$ \begin{equation}\varepsilon_0=10^{-4}, \ \ \varepsilon_k=\left\{\begin{array}{ll} 0.1\varepsilon_{k-1}&\mbox{若}\ \varepsilon_{k-1}>\varepsilon, \\ \varepsilon_{k-1}&\mbox{若}\ \varepsilon_{k-1} <\varepsilon. \end{array}\right.\end{equation} $

(3.6)

算法3.1中的问题(3.4) 是一个结构约束矩阵最优化问题, 对任意给定的$A, B, C, D, E$, $F, \overline{Z}_k$和标量$\rho_k$, 该矩阵最优化问题总是有解的.本文利用谱投影梯度(SPG)方法^[2]求解问题(3.4).所谓谱投影梯度是求解凸集上光滑函数最小值的一种非单调投影梯度方法.在求解问题(3.4) 中, 谱投影梯度方法格式简单、易编程并且不需要利用矩阵分解, 克服了传统的谱步长非单调全局搜索梯度方法收敛慢的不足, 是求解大规模凸约束问题的有效方法.本文将SPG算法做适当改变, 用于求解问题(3.4) 的近似估计值.其迭代格式可以描述为如下步骤.

算法3.2 (SPG方法计算子问题(3.4) 的近似估计矩阵$X_k$)

第1步:输入矩阵$A, B, C, D, E$和$F$; 输入整数$M>1$, 充分小的参数$\alpha_{\min}>0$, $\widetilde{\gamma}\in (0, 1)$, 充分大的参数$\alpha_{\max}>\alpha_{\min}$, 和安全参数$0 < \sigma_1 < \sigma_2 < 1$.选取初始参数$\alpha_1\in [\alpha_{\min}, \alpha_{\max}]$, 初始矩阵$X_1\in {\cal S}_{R, S}^{m\times n}$和迭代精度$\varepsilon_k\downarrow 0$.令$i\leftarrow 1$.

第2步:若$\|\nabla L_{\rho_k}(X_k, \overline{Z}_k)+R\nabla L_{\rho_k}(X_k, \overline{Z}_k)S\| < \varepsilon_k$, 程序运行终止(这时$X_i$是问题(3.4) 的一个近似估计矩阵).

第3步:计算$dX_i=-\frac{1}{2}\alpha_i(\nabla L_{\rho_k}(X_i, \overline{Z}_k)+R\nabla L_{\rho_k}(X_i, \overline{Z}_k)S)$.令$\lambda=1.$

第4步:计算$X_+=X_i+\lambda dX_i$.

第5步:若

$ \begin{equation} L_{\rho_k}(X_+, \overline{Z}_k)\leq \max\limits_{1\leq j\leq \min\{i, M\}}L_{\rho_k}(X_{i-j}, \overline{Z}_k)+\widetilde{\gamma}\lambda\langle dX_i, \nabla L_{\rho_k}(X_i, \overline{Z}_k)\rangle, \end{equation} $

(3.7)

定义$\lambda_i=\lambda$, $X_{i+1}=X_+$, $s_i=X_{i+1}-X_{i}$, $y_i=\nabla L_{\rho_k}(X_{i+1}, \overline{Z}_k)-\nabla L_{\rho_k}(X_i, \overline{Z}_k)$, 转第6步.

若(3.7) 式不成立, 定义$\lambda_{new}\in [\sigma_1\lambda, \sigma_2\lambda]$, 令$\lambda=\lambda_{new}$转第4步.

第6步:计算$b_i=\langle s_i, y_i\rangle$.若$b_i\leq 0$, 令$\alpha_i=\alpha_{\max}$.否则, 计算

$ a_i=\langle s_i, s_i\rangle, \ \alpha_i=\min\{\alpha_{\max}, \max\{\alpha_{\min}, a_i/b_i\}\}. $

第7步:令$i\leftarrow i+1$转第2步.

在运行算法3.2时, 选取整数$M=15$, $\widetilde{\gamma}=10^{-4}$, $\alpha_{\min}=10^{-30}$, $\alpha_{\max}=10^{30}$, $\sigma_1=0.1$, $\sigma_2=0.9$, $\lambda_{new}=(\sigma_1\lambda+\sigma_2\lambda)/2$和$\alpha_0=1$.迭代精度$\varepsilon_k$的选取如(3.6) 式所示.第一次运行算法3.2时的初始矩阵$X_1$取为零矩阵, 第$k$次($k\neq 1$)运行的初始矩阵$X_1$取为算法3.1的第$k-1$次近似解.

定理3.2 若由算法3.1产生的数列$\{\rho_k\}$有界, 则由算法3.1产生的矩阵序列$\{X_k\}$的极限$X^*$是问题(1.2) 的解.

证 注意到$\varepsilon_k\downarrow 0$以及$\|\nabla L_{\rho_k}(X_k, \overline{Z}_k)+R\nabla L_{\rho_k}(X_k, \overline{Z}_k)S\| < \varepsilon_k$, 则有

$ \begin{eqnarray*} && \lim\limits_{k\rightarrow \infty}\|{\cal A}^T({\cal A}( X_k)-C)-\rho_k E^T\left(D-EX_kF+\frac{\overline{Z}_k}{\rho_k}\right)_+F^T\\ &&+R{\cal A}^T({\cal A}(X_k)-C)S-\rho_k RE^T\left(D-EX_kF+\frac{\overline{Z}_k}{\rho_k}\right)_+F^TS\|=0. \end{eqnarray*} $

记

$ \lim\limits_{k\rightarrow \infty}X_k=X^*, \ \ \lim\limits_{k\rightarrow \infty}\rho_k=\rho^*, \ \ \lim\limits_{k\rightarrow \infty}\overline{Z}_k=\overline{Z}^*, $

则有

$ \begin{eqnarray*} && {\cal A}^T({\cal A}(X^*)-C)-\rho^* E^T\left(D-EX^*F+\frac{\overline{Z}^*}{\rho^*}\right)_+F^T\\ &&+R{\cal A}^T({\cal A}(X^*)-C)S-\rho^* RE^T\left(D-EX^*F+\frac{\overline{Z}^*}{\rho^*}\right)_+F^TS=0. \end{eqnarray*} $

令$Y^*=\rho^*\left(D-EX^*F+\frac{\overline{Z}^*}{\rho^*}\right)_+$, 则有

$ Y^*\geq 0, \ \ {\cal A}^T({\cal A}(X^*)-C)+R{\cal A}^T( {\cal A}(X^*)-C)S-E^TY^*F^T-RE^TY^*F^TS=0. $

因为$\{\rho_k\}$有界, 则存在$k_0\in {\Bbb N}$使得不等式(3.5) 对一切$k\geq k_0$成立.因此, 有

$ \lim\limits_{k\rightarrow \infty}\|\{EX_kF-D, \ Z_k\}_{-}\|=0. $

注意到$Z_k\geq 0$对任意整数$k\in {\Bbb N}$成立, 则有

$ \lim\limits_{k\rightarrow \infty}\|(D-EX_kF)_+\|=0, \ \ \lim\limits_{k\rightarrow \infty}\langle Z_k, \ EX_kF-D\rangle=0. $

即有

$ EX^*F-D\geq 0, \ \ \langle Z^*, \ EX^*F-D\rangle=0. $

由$Z_k$和$Y^*$的定义知$(Z^*)_{ij}>0$当且仅当$(Y^*)_{ij}>0$ $(i=1, 2, \cdots, p, j=1, 2, \cdots, q)$.因此

$ \langle Y^*, EX^*F-D\rangle=0. $

综上所述$X^*$满足条件(3.1).由定理3.1知$X^*$是问题(1.2) 的一个解.

本节所有数据均是通过MATLAB 7.0, 机器精度$10^{-16}$, Intel$^{\circledR}$ Core i3处理器, 2.13GHz的PC中获得.为方便计算, 本节例题中的实对合矩阵均取为反序单位矩阵, 即$R=J_m$, $S=J_n$.

例4.1 本例通过考虑问题(1.2) 的如下简化形式

$ \begin{array}{*{20}{l}} {{\rm{minimize}}\;\frac{{\rm{1}}}{{\rm{2}}}{{\left\| {AXB{\rm{ - }}C} \right\|}^{\rm{2}}}}\\ {{\rm{subject}}\;{\rm{to}}X \in S_{R,S}^{m \times n},\;\;EXF \ge D} \end{array} $

(4.1)

来说明算法3.1是可行且高效的.令$l=7$, $h=m=n=8$, $p=6$, $t=7$.按如下Matlab方式给定矩阵$A, B, C, D, E$和$F$:

$ \begin{eqnarray*} A={\left[\begin{array}{cccccccc} 1.5& 0.5 & 2.3 &~~-1.4~~& 1.5& ~~-1~~&-1.9& ~~ 0.3\\ 1.1& ~~-1.4~~ & 0.8& 1.4& -0.2 & 1& -0.8 & ~~1.7\\ -3.3& 4.2& -2.4& -4.2 & 0.6& -3& 2.4 &~~ -5.1\\ 3 & 1 & 4.6& -2.8& 3 & -3& -3.8 & ~~ 0.6\\ 1.1& -1.4& 0.8 & 1.4& -0.2& 1& -0.8 & ~~1.7\\ -5.5& 7 & -4 & -7 & 1 & -5& 4 &~~ -8.5\\ -1.5& -0.5 & -2.3& 1.4& -1.5 & 1& 1.9 & ~~-0.3 \end{array}\right]}, \end{eqnarray*} $

$ \begin{eqnarray*} E={\left[\begin{array}{cccccccc} 0.1& 1.4 & 0.2& ~~ 0.2~~& 0.1& ~~ 0.0~~ & 1.6 & ~~-1.4\\ 0.2& ~~-0.9~~& 1.3& ~~-0.2~~&-0.1 & 1.1& 0.7 & ~~0.4\\ -0.6& -0.5& 1.3& -0.9& -0.0 & -0.0& 0.0 & ~~ -1.4\\ -1.0& 1.4& -1.2& 0.0& -0.6& ~~ -0.0~~& 1.7& ~~-0.9\\ 0.9& 1.4& 1.6& -0.2& 2.0& -0.5& 0.0 & ~~0.9\\ -0.2& -1.0& 0.6& 2.3& 0.6& -0.6& 0.2 & ~~0.3 \end{array}\right]}, \end{eqnarray*} $

$ \begin{eqnarray*} F={\left[\begin{array}{cccccccc} -0.8 & ~~-1.7~~ & 1.9 &~~-0.6 ~~&-0.6& ~~ 1.2~~ & 0.4\\ -0.3 & 0.7& -0.5& 0.3 & 2.0 & ~~-0.9~~& -0.4\\ 0.2 & -1.0& 1.7& -0.8 & -1.7 & -0.1& 0.2\\ 0.2 & -1.1& 0.8& 0.1 & -0.6 & -0.4 & 0.9\\ -1.7 & 1.7& 0 & 1.5 & 0.7 & -0.1 & 0.9\\ -1.2& 2.1& -0.3 & 0.5 & -0.8& 0.3 & 0.6\\ -1.7 & -0.8& 1.6 & -0.7 & 0 & -0.6 & 0\\ 0 & 0& 0.7& 0.3 & 0.6 & 0.8 & -0.6 \end{array}\right]}. \end{eqnarray*} $

$ \begin{eqnarray*} D={\left[\begin{array}{cccccccc} -0.5783& ~~-0.1974~~&-1.9147&~~-3.4000~~& 1.3566& ~~ -2.3000~~& -1.2521\\ -0.4000& -1.5000& -11.5811 & -1.1148& 0.6088 & -0.3737& -3.7877\\ -0.3000& 0& -9.9539& -1.5447& -0.2000 & -0.6899& -4.3020\\ 0.8220& -0.1394& -0.7963 & -0.3000& -0.9000 & 0.7052 & -0.1000\\ 1.4000& -0.0775& -9.4103& -0.2459& -1.0787 & 0.0435 & -1.8165\\ 0.5000& 0.1307& -4.4220& -0.6175& 0.3000 & -0.3000 & -1.8059 \end{array}\right]}, \end{eqnarray*} $

$B=eye(8, 8)$和$C=ones(7, 8)$, 其中$eye(8, 8)$为8阶的单位矩阵, $ones(7, 8)$为$7\times 8$的元素全为1的矩阵.利用算法3.1通过15次外迭代步, 耗时6.7381s得到近似解

$ \begin{eqnarray*} X^{*}=\left[\begin{array}{cccccccc} 2.8328 & 0.5349 & 3.6081& ~-3.8602~& 5.3070 &~-2.4549 ~& 0.3851 & ~-1.4725\\ -0.2654 & ~-0.3707~ & -0.4740& 0.3678& -0.4180 & 0.0457& -0.3579& ~ 0.1037\\ -1.0000 & -1.0000 & -1.0000 & -1.0000 & -1.0000 & -1.0000& -1.0000& ~-1.0000\\ 0.3640 & -0.3544& 0.6625 & -1.9444 & 1.1986 & -1.4163& -0.4058 & ~-1.1121\\ -1.1121& -0.4058& -1.4163 & 1.1986 & -1.9444 & 0.6625& -0.3544 & ~0.3640\\ -1.0000& -1.0000& -1.0000 & -1.0000 & -1.0000 & -1.0000& -1.0000 & ~-1.0000\\ 0.1037& -0.3579& 0.0457 & -0.4180 & 0.3678 & -0.4740& -0.3707& ~-0.2654\\ -1.4725& 0.3851 & -2.4549 & 5.3070 & -3.8602 & 3.6081& 0.5349 & ~2.8328 \end{array}\right], \end{eqnarray*} $

由$X^*$的元素排列知$X^*$是$(R, S)$对称矩阵(中心对称), 经计算得

(a) $\|AX^{*}B-C\|=6.3246, $

(b) $\|A^T(AX^{*}B-C)B^T+RA^T(AX^{*}B-C)B^TS\|=5.9005\times 10^{-7}, $

(c) $\left\|{A}^T(AX^*B-C)B^T+R{A}^T(AX^*B-C)B^TS-E^TY^*F^T-RE^TY^*F^TS\right\|=2.0351\times 10^{-7}, $

(d) $\|(EX^{*}F-{D})_{-}\|=3.8527\times 10^{-10}, $

(e) $\langle Y^*, EX^*F-D\rangle=1.3522\times 10^{-13}.$

由(a)(b)可知$X^*$是矩阵方程$AXB=C$的最小二乘$(R, S)$对称解(因为$\nabla {\cal F}(X^*)$在${\cal S}_{R, S}^{m\times n}$上的投影为0), 由(c)-(e)可知$X^*$满足定理3.1中的对应于问题(4.1) 的条件(3.1).

例4.2 本例考虑问题(1.2) 的如下形式

$ \begin{equation} \begin{array}{rl} \rm{minimize}& \frac{1}{2}\left\|A_1XB_1+A_2XB_2-C\right\|^2\\ \rm{subject}\ \rm{to} & X\in {\cal S}_{R, S}^{m\times n}, \ \ EXF\geq D. \end{array} \end{equation} $

(4.2)

按如下Matlab方式生成矩阵

$ \begin{eqnarray*}\begin{array}{lll} A_1=randn(l, m), &B_1=randn(n, h), & C=A\overline{X}B, \\ E=randn(p, m), &\ F=randn(n, t), & D=E\overline{X}F-ones(p, t), \end{array} \end{eqnarray*} $

其中$randn(m, n)$表示元素服从标准分布的$m\times n$随机矩阵, $\overline{X}=\overline{Y}+R\overline{Y}S\in {\cal S}_{R, S}^{m\times n}$, $\overline{Y}=rand(m, n)$.矩阵$C$的选取保证矩阵方程$AXB=C$在${\cal S}_{R, S}^{m\times n}$内是相容的, 矩阵$D$的选取保证问题(4.2) 是相容的.

当$l\geq m$且$h\geq n$时, 问题(4.1) 只有唯一解$\overline{X}$.对于不同系统维数产生的随机系数矩阵$A, B, E, F$, 表 1给出了算法3.1收敛的迭代时间和相对误差, 这里相对误差定义为

$ Re_{error}=\frac{X_k-\overline{X}}{\overline{X}}, $

其中$X_k$为算法收敛所得的近似解. 表 2给出了当系数维数取为$l=80$, $h=90$, $m=80$, $n=80$, $p=80$和$t=80$时的收敛数据, 其中$k$表示外迭代(算法3.1) 次序, $Error_{Algorithm 3.1}$表示外迭代精度, 即算法3.1中对应于$X_k$和$Z_k$的指标值

$ (\|\nabla L_{\rho_k}(X_k, \overline{Z}_k)+R\nabla L_{\rho_k}(X_k, \overline{Z}_k)S\|^2+\|\{EX_kF-D, Z_k\}_{-}\|^2)^{1/2}, $

$Inner_{Steps}$表示对应于外迭代第$k$步所需的内迭代(算法3.2) 次数, $Error_{Algorithm3.2}$表示内迭代精度, 即算法3.2的终止精度指标值

$ \|\nabla L_{\rho_k}(X_k, \overline{Z}_k)+R\nabla L_{\rho_k}(X_k, \overline{Z}_k)S\|. $

当$l < m$或$h < n$时, 问题(4.1) 有无穷多解, 由算法得到的序列$\{X_k\}$不一定会收敛到$\overline{X}$.在这种情况下, 对于不同系统维数产生的随机系数矩阵$A, B, E, F$, 表 3给出了算法3.1收敛的迭代时间、剩余范数$\|A_1X_kB_1+A_2X_kB_2-C\|$和$\|(EX_kF-D)_-\|$. 表 4给出了当系数维数取为$l=80$, $h=90$, $m=80$, $n=80$, $p=80$和$t=80$时的收敛数据, 其中$k$表示外迭代(算法3.1) 次序, $Error_{Algorithm 3.1}$表示外迭代精度, $Inner_{Steps}$表示对应于外迭代第$k$步所需的内迭代(算法3.2) 次数, $Error_{Algorithm 3.2}$表示内迭代精度.

本文讨论了矩阵不等式$EXF\geq D$约束下广义Sylvester矩阵方程$\sum_{i=1}^{q}A_iXB_i=C$的$(R, S)$对称最小二乘问题, 其中$EXF\geq D$表示矩阵$EXF-D$非负.结合$(R, S)$对称拉直算子矩阵, 本文将矩阵优化问题看作为带不等式约束的凸二次规划问题, 给出问题有解的充分必要条件, 进而将优化理论中求凸集上光滑函数最小值的增广Lagrangian方法应用于求解该问题, 并给出矩阵形式的迭代格式, 避免拉直算子引起的问题规模大幅增加的缺点.数值例子中以该问题的简化形式-问题(4.1) 为例说明了矩阵形式的增广Lagrangian方法是可行且高效的.特别是当问题只有唯一解时, 算法的迭代效率是非常高的.