本文考虑如下带有列正交约束的广义Sylvester方程极小化问题.

问题1.1  给定矩阵$ A_1, A_2, \ldots A_N\in {\Bbb R}^{l\times n} $, $ B_1, B_2, \ldots B_N\in {\Bbb R}^{p\times s} $$ (n\gg p) $和$ C\in {\Bbb R}^{l\times s} $. 求$ X\in {\Bbb R}^{n\times p} $使得

(1.1) $ \begin{equation} \begin{array}{rc} \min& { } \frac{1}{2}\Big\| \sum\limits_{i = 1}^{N}A_iXB_i-C\Big\|^{2}\\ \rm{s.t.}& \ X^TX = I_p, \end{array} \end{equation} $

其中$ \|\cdot\| $表示Frobenius范数.

带约束条件的广义Sylvester方程求解问题或最小二乘问题因在图像处理、控制论等领域的广泛应用, 近年来已被充分研究. 文献[1-4] 利用矩阵形式的Krylov子空间法给出了结构约束下广义Sylvester方程迭代解法. 文献[5]利用递阶辨识原理给出了无约束下广义Sylvester方程递阶辨识迭代法, 其实质为最速下降法的推广. 文献[6-7]利用谱投影梯度法和条件梯度法研究极小化问题

(1.2) $ \begin{equation} \rm{min}\ \ \frac{1}{2}\Big\|\sum\limits_{i = 1}^{q}A_iXB_i-C\Big\|^2, \ \ \ \rm{s.t.}\ X\in \Omega, \end{equation} $

其中$ \Omega_1 = \{X: L\leq X\leq U\} $, $ \Omega_2 = \{X: \|X\|\leq \delta\} $和$ \Omega = \Omega_1\cap\Omega_2 $, 并将该方法应用于一类带正则化项的图像恢复问题. 文献[8]利用非精确的交替方向法研究了问题(1.2), 其中$ \Omega = \{X: L\leq EXF\leq U\} $和$ X\in {\cal S} $, 且$ {\cal S} $为具有线性结构的线性子空间, 并将模型和方法应用于带正则化项的具有特殊对称性的图像恢复问题. 更多不同约束条件下广义Sylvester方程及其特殊方程形式的求解参见论文[9-13], 若问题1.1中的$ A_i $和$ B_i $均为单位阵, 则问题1.1退化为约束矩阵最佳逼近问题, 相关的研究也很多, 这里不再细述. 但是正交约束下的广义Sylvester方程或其特殊方程形式的求解问题相关研究成果并不多. 特别地, 文献[14]利用$ k $阶反对称主子阵的$ n(n>k) $阶正交矩阵的C-S分解研究了反对称正交约束矩阵方程$ AX = B $的求解. 但是该方法不能应用于正交约束下具有一般形式的广义Sylvester方程求解问题.

注意到问题(1.1)的可行域$ \{X\in {\Bbb R}^{n\times p}| X^TX = I_p, \; n\gg p\} $通常被称为Stiefel流形, 记为$ {\rm{St}}(n, p) $. 本文从黎曼流形的角度给出问题1.1的有效数值求解算法. 目前, 经过国内外学者的不断探讨, Stiefel流形优化问题的研究已经取得一系列的成果. 现有的方法主要有可行方法和不可行方法. 不可行方法的特点是迭代点不一定保持列正交约束, 但能逐渐减少迭代点的不可行度, 其主要方法有Majorization算法^[15], 交替方向法^[16], 分裂正交约束法^[17]和邻近点交替极小增广Lagrangian算法^[18]等. 可行方法的主要特点是迭代更新点始终保持列正交约束性质, 并使目标函数逐渐下降直至局部极小点. 可行方法主要有基于流形的黎曼梯度下降法、黎曼共轭梯度法、黎曼牛顿法和信頼域法等, 而其中黎曼共轭梯度法因是欧氏空间上非线性共轭梯度法的流形推广, 近年来应用较广^[19-24]. 如文献[19]构造黎曼非线性共轭梯度法求解Stiefel流形上的矩阵奇异值计算问题. 文献[21] 构造黎曼共轭梯度算法求解定秩矩阵流形上的低秩矩阵完全化问题. 文献[22-24] 分别将欧氏空间上的修正Polak-Ribière-Polyak(MPRP)非线性共轭梯度法^[25]和改进Fletcher-Reeves非线性共轭梯度法推广至黎曼乘积流形, 解决随机矩阵逆特征值问题和双随机逆特征值问题; 文献[20] 基于Cayley变换和差分收缩, 将欧氏空间上的Dai非单调共轭梯度法推广至黎曼流形.

基于欧氏空间上的MPRP非线性共轭梯度法, 文献[25]和^[22]给出的黎曼乘积流形上的推广形式, 本文构造一类黎曼MPRP共轭梯度算法求解列正交约束下的广义Sylvester方程极小化问题: 问题1.1, 并给出算法全局收敛性. 该迭代格式的主要特点是搜索方向恒能保证该目标函数下降. 数值实验表明, 相比于基于流形的黎曼梯度下降法和几类常用的不可行方法, 该算法对于问题模型(1.1)是高效可行的. 本文余下内容组织如下: 第2节给出黎曼共轭梯度法的基本框架, 并结合Stiefel流形的几何性质, 给出求解问题模型(1.1)的黎曼MPRP共轭梯度法的迭代格式; 第3节给出算法的全局收敛性分析; 第4节给出数值算例和数值比较验证所提算法的可行性和高效性. 最后第5节给出结论.

欧氏空间上的MPRP非线性共轭梯度法^[25]主要用于求解如下无约束优化问题

$ \begin{eqnarray*} \min\limits_{x\in R^n} f(x). \end{eqnarray*} $

该算法的迭代格式如下

(2.1) $ \begin{equation} x_{k+1} = x_{k}+\alpha_k\eta_k, \ \ k = 0, 1, \cdots, \end{equation} $

其中搜索方向

(2.2) $ \begin{equation} \eta_k = \left\{ \begin{array}{lr} -\nabla f(x_k), & {\rm{if}}\ k = 0;\\ -\nabla f(x_k)+\beta_k^{PRP} \eta_{k-1}-\theta_k y_{k-1}, & {\rm{if}}\ k>0, \end{array} \right. \end{equation} $

其中

(2.3) $ \begin{equation} \beta_k^{PRP} = \frac{\nabla f(x_k)^Ty_{k-1}}{\|\nabla f(x_{k-1})\|^2}, \; \; y_{k-1} = \nabla f(x_k)-\nabla f(x_{k-1}), \; \; \theta_k = \frac{\nabla f(x_k)^T\eta_{k-1}}{\|\nabla f(x_{k-1})\|^2}. \end{equation} $

由(2.2)和(2.3)式易得$ \eta_k\nabla f(x)^T = -\|\nabla f(x)\|^2 $, 这说明搜索方向$ \eta_k $是目标函数$ f $在$ x_k $处的下降方向. 迭代格式中步长$ \alpha_k = \overline{\alpha}_k\lambda^t $, 其中$ t $是满足如下不等式的最小非负整数

(2.4) $ \begin{equation} f(x_{k}+\alpha_k\eta_k)\leq f(x_k)-\delta\alpha_k^2 \|\nabla f(x_k)\|^2, \end{equation} $

其中$ \overline{\alpha}_k $表示第$ k $步的初始步长, 具体表示形式在下文中(2.27)式. $ \delta, \lambda\in (0, 1) $.

为了解决问题(1.1), 需将欧氏空间上的MPRP非线性共轭梯度法推广至黎曼流形$ {\cal M} $. 在欧氏空间上, 新的迭代点$ x_{k+1} $可通过平移得到, 但在黎曼流形是不能实现的. 为了能使新的迭代点属于流形$ {\cal M} $, 利用黎曼优化$ \min\limits_{x\in {\cal M}} f(x) $中的一般迭代格式

(2.5) $ \begin{equation} x_{k+1} = {\cal R}_{x_k}(\alpha_k\eta_k), \end{equation} $

其中$ \eta_k $为位于流形$ {\cal M} $上点$ {x_k} $处切平面$ T_{x_k}{\cal M} $内的搜索方向, $ {\cal R} $表示$ T_{x_k}{\cal M}\rightarrow {\cal M} $的收缩映射(retraction)^[26]. 在黎曼流形$ {\cal M} $上, 黎曼梯度$ \nabla f(x_k)\in T_{x_k}{\cal M} $, 但搜索方向$ \eta_{k-1}\in T_{x_{k-1}}{\cal M} $, 而位于不同切平面的向量不能直接相加. 由此可知(2.1)-(2.3)式无法直接推广到黎曼流形$ {\cal M} $. 为了能够实现加法运算, 需要通过向量转移算子(vector transport)^[26]将$ \eta_{k-1} $和$ \nabla f(x_{k-1}) $映射到$ x_k $的切平面$ T_{x_k}{\cal M} $上, 分别记为

$ {\cal T}_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})\; \; \; \; \rm{和}\; \; \; \; {\cal T}_{\alpha_{k-1}\eta_{k-1}}(\nabla f(x_{k-1})). $

通过上述变化, 可得到搜索方向

(2.6) $ \begin{equation} \eta_k = \left\{ \begin{array}{lr} -\nabla f(x_k), & {\rm{if}}\ k = 0;\\ -\nabla f(x_k)+\beta_k {\cal T}_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})-\theta_k y_{k-1}, & {\rm{if}}\ k>0, \end{array} \right. \end{equation} $

这里迭代参数

(2.7) $ \begin{equation} y_{k-1} = \nabla f(x_k)-{\cal T}_{\alpha_{k-1}\eta_{k-1}}(\nabla f(x_{k-1})), \end{equation} $

(2.8) $ \begin{equation} \beta_k = \frac{\nabla f(x_k)^Ty_{k-1}}{\|\nabla f(x_{k-1})\|^2}, \; \; \theta_k = \frac{\nabla f(x_k)^T{\cal T}_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})}{\|\nabla f(x_{k-1})\|^2}. \end{equation} $

给定第$ k $步初始步长$ \overline{\alpha}_k $, 参数$ \delta, \lambda\in (0, 1) $以及搜索方向$ \eta_k $, 步长$ {\alpha}_k = \overline{\alpha}_k\lambda^t $满足如下Armijo法则

(2.9) $ \begin{equation} f({\cal R}_{x_k}(\alpha_k\eta_k))\leq f(x_k)-\delta\alpha_k^2 \|\nabla f(x_k)\|^2. \end{equation} $

下面针对问题(1.1)的可行域介绍Stiefel流形$ {\rm{St}}(n, p) $的基础知识, 更多性质可参考文献[26]. $ \forall X\in {\rm{St}}(n, p) $, 点$ X $处的切空间可表示为

$ {\bf T}_X{\rm{St}}(n, p) = \{\xi\in {{\Bbb R}} ^{n\times p}:X^{T}\xi+\xi^{T}X = 0.\} $

显然点$ X $的切空间$ {\bf T}_X{\rm{St}}(n, p) $为与$ X $相关的线性子空间, 故对$ \forall \xi, \zeta\in {\bf T}_X{\rm{St}}(n, p) $, 其内积可定义为$ \langle \xi, \zeta\rangle_X = {\rm{trace}}(\xi^T\zeta) $, 诱导范数$ \|.\|_X $即为欧氏空间中的Frobenius范数. 为方便内积和范数可简记为$ \langle\cdot, \cdot\rangle $和$ \|\cdot\| $. $ \forall M\in R^{n\times p} $, $ M $在点$ X\in {\rm{St}}(n, p) $的切空间$ {\bf T}_X{\rm{St}}(n, p) $的正交投影可表示为

(2.10) $ \begin{equation} {\bf P}_{X}(M) = (I-XX^T)M+X {\rm{skew}}(X^TM) = M-X{\rm{sym}}(X^TM), \end{equation} $

其中对$ \forall N\in R^{p\times p} $, $ {\rm{skew}}(N) = \frac{1}{2}(N-N^T), \ {\rm{sym}}(N) = \frac{1}{2}(N+N^T). $由迭代格式(2.5)-(2.8) 可知需要用到Stielfel流形上的收缩算子和向量转移算子, 这里应用Stielfel流形上常用的两类算子: 通过QR分解构建收缩算子和正交投影定义向量转移^[26]. $ \forall \xi\in {\bf T}_X{\rm{St}}(n, p) $, 基于QR分解的收缩算子定义为

(2.11) $ \begin{equation} {\cal R}_X(\xi) = Qf(X+\xi). \end{equation} $

这里$ Qf(A) $表示非奇异矩阵$ A\in{\Bbb R}^{n\times p} $的QR分解$ A = \widetilde{Q}\widetilde{R} $中的Q因子, 其中$ \widetilde{Q}\in {\rm{St}}(n, p) $, $ \widetilde{R}\in{\Bbb R}^{p\times p} $为对角元全为正的上三角阵. $ \forall \eta, \xi\in {\bf T}_X{\rm{St}}(n, p) $, 基于正交投影的向量转移定义为

$ {\cal T}_{\xi}(\zeta) = {\bf P}_{{\cal R}_X(\xi)}(\zeta), $

即$ \zeta $在点$ {\cal R}_X(\xi) $处的切平面$ {\bf T}_{{\cal R}_X(\xi)}{\rm{St}}(n, p) $内的正交投影. 结合切平面正交投影公式(2.10), 可得该向量转移的具体计算公式为

(2.12) $ \begin{equation} {\cal T}_{\xi}(\zeta) = \zeta-Qf(X+\xi){\rm{sym}}\left(Qf(X+\xi)^T\zeta\right). \end{equation} $

下面讨论黎曼梯度的具体计算公式. 为方便记$ {f}(X) $为问题(1.1)中的目标函数, 即

$ f(X): = \frac{1}{2}\Big\| \sum\limits_{i = 1}^{N}A_iXB_i-C\Big\|^{2}, \; \; \; X\in {\rm{St}}(n, p), $

记$ \overline{f}(X) $为定义在$ {{\Bbb R}} ^{n\times p} $上与$ {f}(X) $具有相同表达式的函数, 即

$ \overline{f}(X): = \frac{1}{2}\Big\| \sum\limits_{i = 1}^{N}A_iXB_i-C\Big\|^{2}, \; \; \; X\in {{\Bbb R}}^{n\times p}. $

同时定义

$ F(X): = \sum\limits_{i = 1}^{N}A_iXB_i-C, \; \; \; X\in {\rm{St}}(n, p) $

和

$ \overline{F}(X): = \sum\limits_{i = 1}^{N}A_iXB_i-C, \; \; \; X\in {{\Bbb R}}^{n\times p}, $

则$ {f}(X) $和$ {F}(X) $分别可视为$ \overline{{f}}(X) $和$ \overline{{F}}(X) $在流形$ {\rm{St}}(n, p) $上的限制, 即$ {f}(X) = \overline{{f}}(X)|_{{\rm{St}}(n, p)} $和$ {F}(X) = \overline{{F}}(X)|_{{\rm{St}}(n, p)} $. 因为$ {\rm{St}}(n, p) $可视为欧氏空间$ {{\Bbb R}} ^{n\times p} $的一类嵌入子流形, 由文献[26, p48] 可知, $ f $关于$ X\in {\rm{St}}(n, p) $的黎曼梯度$ {\rm{grad}} f(X) $可$ f $关于$ X\in {{\Bbb R}}^{n\times p} $的欧氏梯度$ \nabla \overline{f}(X) $在点$ X $切平面的正交投影, 即

(2.13) $ \begin{equation} {\rm{grad}} {f}(X) = {\bf P}_{X}\left(\nabla \overline{{f}}(X)\right) = \nabla \overline{{f}}(X)-X {\rm{sym}}\left(\nabla \overline{{f}}(X)^{T}X\right), \end{equation} $

其中$ \overline{f} $在$ X\in {{\Bbb R}} ^{n\times p} $的欧氏梯度可通过简单计算得

(2.14) $ \begin{equation} \nabla\overline{f}(X) = \frac{\partial\overline{f}(X)}{\partial X} = \sum\limits_{j = 1}^{N}A_j^T\left(\sum\limits_{i = 1}^{N}A_iXB_i-C\right)B_j^T. \end{equation} $

下面计算$ F(X) $的微分公式及下文将会用到的回调映射. $ \forall X\in {\rm{St}}(n, p), \xi\in {\bf T}_X{\rm{St}}(n, p) $有

$ \begin{eqnarray*} \overline{F}(X+t\xi)-\overline{F}(X) = \left(\sum\limits_{i = 1}^{N}A_i(X+t\xi)B_i-C\right)-\left(\sum\limits_{i = 1}^{N}A_iXB_i-C\right) = t\sum\limits_{i = 1}^{N}A_i\xi B_i. \end{eqnarray*} $

由于$ \overline{F} $是$ {\Bbb R}^{n\times p}\rightarrow {\Bbb R}^{l\times s} $的光滑映射, 而$ F $是$ \overline{F} $在$ {\Bbb R}^{n\times p} $的嵌入黎曼子流形$ {\rm{St}}(n, p) $上的限制, 故可知函数$ F $在点$ X\in {\rm{St}}(n, p) $上的微分$ {\cal D}F(X):{\bf T}_X{\rm{St}}(n, p)\longrightarrow {\bf T}_{F(X)}{\Bbb R}^{l\times l} \simeq {\Bbb R}^{l\times l} $可表示为^[26]

(2.15) $ \begin{equation} {\cal D}F(X)[\xi] = {\cal D}\overline{F}(X)[\xi] = \sum\limits_{i = 1}^{N}A_i\xi B_i, \; \; \; \; \xi\in {\bf T}_X{\rm{St}}(n, p). \end{equation} $

此外$ f $的黎曼梯度与$ F(X) $的微分存在如下关系^[26]

(2.16) $ \begin{equation} {\rm{grad}} f(X) = {\cal D}F(X)^{\ast}[F(X)], \end{equation} $

这里$ {\cal D}F(X)^{\ast} $表示$ {\cal D}F(X) $的伴随算子. 记$ {\bf T}{\rm{St}}(n, p) $为流形$ {\rm{St}}(n, p) $上的切向量丛, 定义回调映射$ \widehat{f}:{\bf T} {\rm{St}}(n, p)\longrightarrow {\Bbb R} $和$ \widehat{F}:{\bf T} {\rm{St}}(n, p)\longrightarrow {\Bbb R}^{l\times s} $如下

(2.17) $ \begin{equation} \widehat{f}(\xi) = f\left({\cal R}(\xi)\right), \; \; \; \; \widehat{F}(\xi) = F\left({\cal R}(\xi)\right), \ \ \ \ \xi\in {\bf T} {\rm{St}}(n, p). \end{equation} $

记$ \hat{f}_{X} $和$ \hat{F}_{X} $分别表示回调映射$ \hat{f} $和$ \hat{F} $在点$ X $切平面$ {\bf T}_X{\rm{St}}(n, p) $上的限制, 则有

(2.18) $ \begin{equation} {\cal D}F(X) = {\cal D}\widehat{F}_X(0_X), \ \ \ \ X\in {\rm{St}}(n, p), \end{equation} $

其中$ 0_{X} $表示切平面$ {\bf T}_X{\rm{St}}(n, p) $上的原点^[26]. 且$ f $的黎曼梯度与其回调函数$ \hat{f} $的欧氏梯度具有如下关系^[26]

(2.19) $ \begin{equation} {\rm{grad}} f(X) = \nabla \hat{f}_X(0_{X}). \end{equation} $

根据以上讨论, 求解问题1.1的黎曼MPRP共轭梯度算法描述如算法1所示.

算法1  求解问题1.1的黎曼MPRP共轭梯度法.

给定初值$ X_0\in {\rm{St}}(n, p) $. 令$ \varepsilon, \rho, \delta\in(0, 1) $, $ \eta_0 = -{\rm{grad}}f(X_0) $, $ k = 0 $. {$ \| {\rm{grad}} f(X_k)\| >\varepsilon $ } 令

(2.20) $ \begin{equation} \eta_k = -{\rm{grad}}f(X_k)+\beta_{k}{\cal T}_{\alpha_{k-1}\eta_{k-1}}\eta_{k-1}-\theta_kY_{k-1}, \end{equation} $

其中

(2.21) $ \begin{equation} Y_{k-1} = {\rm{grad}} f(X_k)-{\cal T}_{\alpha_{k-1}\eta_{k-1}}{\rm{grad}} f(X_{k-1}), \end{equation} $

(2.22) $ \begin{equation} \beta_{k} = \frac{\langle{\rm{grad}} f(X_k), Y_{k-1}\rangle}{\|{\rm{grad}} f(X_{k-1})\|^2}, \end{equation} $

(2.23) $ \begin{equation} \theta_k = \frac{\langle {\rm{grad}} f(X_k), {\cal T}_{\alpha_{k-1}\eta_{k-1}}\eta_{k-1}\rangle}{\|{\rm{grad}}f(X_{k-1})\|^2}. \end{equation} $

选取步长$ \alpha_k = \max\{\overline{\alpha}\rho^{j}, j = 0, 1, 2\cdots \} $使其满足

(2.24) $ \begin{equation} f\left({\cal R}_{X_k}(\alpha_k\eta_k)\right)\leq f(X_k)-\delta\alpha_{k}^2\|\eta_k\|^2. \end{equation} $

令$ X_{k+1} = {\cal R}_{X_k}(\alpha_k\eta_k) $, 且$ k\leftarrow k+1 $.

注2.1  类似于欧氏空间中的无约束优化问题, 我们可以用$ \|{\rm{grad}}f(X_k)\|\leq\varepsilon $ (这里常数$ \varepsilon>0 $)来作为黎曼优化的停止条件. 事实上, 问题1.1的拉格朗日函数为

$ {\cal L}(X, \Lambda) = \overline{f}(X)-\frac{1}{2}{\rm{trace}}\left(\Lambda(X^TX-I_p)\right), $

其中$ \Lambda $为对应于对称约束条件$ X^TX = I_p $的对称拉格朗日乘子矩阵. $ {\cal L}(X, \Lambda) $分别对$ X $和$ \Lambda $求偏导可得

$ \begin{eqnarray*} \frac{\partial {\cal L}( {X, \Lambda} )}{\partial X} = \nabla \overline{f}(X)-X\Lambda = 0, \ \ \ \frac{\partial {\cal L}( {X, \Lambda} )}{\partial \Lambda} = X^TX-I_p = 0. \end{eqnarray*} $

对上式第一个等式左乘$ X^T $, 可得拉格朗日乘子矩阵$ \Lambda = G^{T}X $. 进而可得问题(1.1) 的一阶优化条件为$ \nabla \overline{f}(X)-X\nabla \overline{f}(X)^TX = 0 $和$ X^TX = I_p $. 事实上, 可证明在$ X^TX = I_p $条件下$ \nabla \overline{f}(X)-X\nabla \overline{f}(X)^TX = 0 $等价于$ {\rm{grad}} f(X) = 0 $, 因为由公式(2.13)

$ {\rm{grad}} f(X) = \nabla \overline{f}(X)-X {\rm{sym}}\left(\nabla \overline{f}(X)^{T}X\right) = \left(I_p-\frac{1}{2}XX^T\right)\left(\nabla \overline{f}(X)-X\nabla \overline{f}(X)^TX\right). $

注2.2  若算法1中每次迭代的初始步长均取为$ \overline{\alpha} $, 收敛效果并不好. 下面参考文献[22, 25]给出一个获取合理初始步长的方法. 对定义在两个切空间$ {\bf T}_{X_k}{\rm{St}}(n, p) $和$ {\bf T}_{F(X_k)}{\Bbb R}^{l\times s}\simeq {\Bbb R}^{l\times s} $的回调映射$ \widehat{F}_{X_k}(\xi) $, 其在点$ 0_{X_k} $处沿着搜索方向$ \eta_k $的一阶近似为

(2.25) $ \begin{equation} \widehat{F}_{X_k}(t\xi_k)\approx \widehat{F}_{X_k}(0_{X_k})+t{\cal D}\widehat{F}_{X_k}(0_{X_k})[\eta_k]. \end{equation} $

当$ t\rightarrow 0 $时, 由(2.17)和(2.18)式知上式等价于

$ F(R_{X_k}(t\eta_k))\approx F(X_k)+t{\cal D}F(X_k)[\eta_k]. $

且当$ t\rightarrow0 $时, 由(2.25)式我们有

$ \begin{eqnarray*} \widehat{f}_{X_k}(t\eta_k)& = &\frac{1}{2}\|\widehat{F}_{X_k}(t\eta_k)\|^2\\ &\approx& \frac{1}{2}\|\widehat{F}_{X_k}(0_{X_k})\|^2+t\langle {\cal D}\widehat{F}_{X_k}(0_{X_k})^{*}[\widehat{F}_{X_k}(0_{X_k})], \; \eta_k\rangle+\frac{1}{2}t^2\|{\cal D}\widehat{F}_{X_k}(0_{X_k})[\eta_k]\|^2. \end{eqnarray*} $

由(2.17)和(2.18)式知当$ t\rightarrow0 $时上式变为

(2.26) $ \begin{equation} f\left({\cal R}_{X_k}(t\eta_k)\right)\approx f(X_k)+t\langle ({\cal D}F(X_k))^{*}[F(X_k)], \eta_k\rangle+\frac{1}{2}t^2\|{\cal D}F(X_k)[\eta_k]\|^2 . \end{equation} $

由$ (2.16) $式知当$ t\rightarrow0 $时$ f $在点$ X_k $处沿着方向$ \eta_k $的局部近似为

$ w_k(t) = f(X_k)+t\langle {\rm{grad}} f(X_k), \eta_k\rangle+\frac{1}{2}t^2\|{\cal D}F(X_k)[\eta_k]\|^2. $

上式两边对$ t $进行求导得

$ w'_k(t) = \langle {\rm{grad}} f(X_k), \eta_k\rangle +t\|{\cal D}F(X_k)[\eta_k]\|^2. $

令$ w'_k(t) = 0 $便有

$ t = -\frac{\langle {\rm{grad}} f(X_k), \eta_k\rangle }{\|{\cal D}F(X_k)[\eta_k]\|^2}. $

结合(2.15)式, 故令初始步长为

(2.27) $ \begin{equation} \overline{\alpha}_k: = \frac{|\langle {\rm{grad}} f(X_k), \; \eta_k\rangle| }{\|\sum\limits_{i = 1}^{N}A_i\eta_k B_i\|^2}. \end{equation} $

首先给出算法1的如下基本性质.

(1) 由(2.20)-(2.23)式可知对$ \forall k\geq 1 $, 有

(3.1) $ \begin{eqnarray} \langle \eta_k, {\rm{grad}}f(X_k)\rangle & = &-\langle{\rm{grad}} f(X_k), \; {\rm{grad}} f(X_k)\rangle-\theta_k \langle Y_{k-1}, \; {\rm{grad}} f(X_k)\rangle{}\\ &&+\beta_{k}\langle {\cal T}_{\alpha_{k-1}\eta_{k-1}}\eta_{k-1}, \; {\rm{grad}} f(X_k)\rangle{}\\ & = &-\|{\rm{grad}} f(X_k)\|^2, \end{eqnarray} $

由此可知搜索方向$ \eta_k $是目标函数$ f $在$ X_k $切平面$ {\bf T}_{X_k}{\rm{St}}(n, p) $内的一个下降方向.

(2) 由于$ f $有下界, 结合(2.24)式可知由算法1产生的函数序列$ \{ f(X_k)\} $是单调下降的, 故为收敛的. 于是有$ \sum\limits_{k = 0}^\infty \alpha_k^2\|\eta_k\|^2 <\infty $, 即

(3.2) $ \begin{equation} \lim\limits_{k\rightarrow \infty}\alpha_k\|\eta_k\| = 0. \end{equation} $

接下来对算法$ 1 $进行全局收敛性分析. 欧氏空间上MPRP共轭梯度法的全局收敛性分析见文献[25], 黎曼乘积流形上黎曼MPRP共轭梯度法的全局收敛性分析见文献[22]. 为论文完整性, 本文给出收敛性分析完整框架.

引理3.1  $ \forall X_0\in {\rm{St}}(n, p) $, 集合

(3.3) $ \begin{equation} {\cal Z} = \{X\in {\rm{St}}(n, p) \mid f(X)\leq f(X_0)\} \end{equation} $

是一个紧集.

证  对$ \forall X\in {\cal Z} $, 有$ \|X\| = \sqrt{p} $, 即集合$ {\cal Z} $是有界的. 又因为目标函数$ f $是一个连续函数, 可知$ {\cal Z} $是闭集. 因此$ {\cal Z} $是一个紧集.

命题3.1  对算法1产生的序列$ \{Y_{k-1}\} $和$ \{\eta_{k-1}\} $, 存在一个常数$ \lambda>0 $, 使得当$ k $足够大时有

$ \|Y_{k-1}\|\leq \lambda \alpha_{k-1}\|\eta_{k-1}\|, $

这里$ \alpha_{k-1} $为算法$ 1 $中定义的步长.

证  由$ (2.12) $和$ (2.21) $式可知

(3.4) $ \begin{eqnarray} Y_{k-1}& = &{\rm{grad}} f(X_k)-{\cal T}_{\alpha_{k-1}\eta_{k-1}}\left({\rm{grad}}f(X_{k-1})\right){}\\ & = &{\rm{grad}}f(X_k)-{\bf P}_{{\cal R}_{X_{k-1}}(\alpha_{k-1}\eta_{k-1})}\left({\rm{grad}} f(X_{k-1})\right){}\\ & = &{\rm{grad}} f(X_k)-{\bf P}_{X_k}\left({\rm{grad}} f(X_{k-1})\right){}\\ & = &{\bf P}_{X_k}\left({\rm{grad}} f(X_k)-{\rm{grad}} f(X_{k-1})\right) \end{eqnarray} $

由于$ {\rm{St}}(n, p) $是$ {\Bbb R}^{n\times p} $的一个嵌入子流形, 且对$ \forall X\in {\rm{St}}(n, p) $都有$ {\bf T}_X{\rm{St}}(n, p)\subset{\Bbb R}^{n\times p} $, 故(2.13) 式p中的黎曼梯度可以间接看作为$ {\rm{St}}(n, p) $和$ {\Bbb R}^{n\times p} $间的连续映射. 于是对$ \forall X_1, X_2\in {\rm{St}}(n, p) $, $ {\rm{grad}} f(X_1), {\rm{grad}} f(X_2) $可看作是$ {\Bbb R}^{n\times p} $中的向量, 故$ {\rm{grad}} f(X_2)-{\rm{grad}}f(X_1) $是有意义的. 我们说$ {\rm{grad}}f(X) $在紧集$ {\cal Z} $上是Lipschitz连续的, 即对任意的$ X_1, X_2\in {\cal Z}\subset {\rm{St}}(n, p) $, 存在一个常数$ \nu_{L_1}>0 $, 使得

(3.5) $ \begin{equation} \|{\rm{grad}} f(X_2)-{\rm{grad}} f(X_1)\|\leq \nu_{L_1} \rm{dist}(X_1, X_2) \end{equation} $

这里的"dist"表示$ {\rm{St}}(n, p) $上的黎曼距离. 由于$ {\rm{St}}(n, p) $是一个紧流形, $ {\cal R} $为定义在$ {\rm{St}}(n, p) $上的收缩算子, 则存在两个常数$ \mu>0, \delta_\mu>0 $, 使得对$ \forall X\in {\rm{St}}(n, p) $, $ \eta_X\in {\bf T}_X{\rm{St}}(n, p) $及$ \|\eta_X\|\leq\delta_\mu $, 有

(3.6) $ \begin{equation} \mu\|\eta_X\|\geq \rm{dist}({\cal R}_X(\eta_X), X). \end{equation} $

由(3.1)式知对$ k $足够大有

(3.7) $ \begin{equation} \|\alpha_{k-1}\eta_{k-1}\|\leq\delta_\mu. \end{equation} $

于是由(3.4)-(3.7) 式知对$ k $足够大有

$ \begin{eqnarray*} \|Y_{k-1}\|& = &\|{\bf P}_{X_k}({\rm{grad}} f(X_k)-{\rm{grad}}f(X_{k-1}))\|\\ &\leq&\|{\rm{grad}} f(X_k)-{\rm{grad}} f(X_{k-1})\|\\ &\leq&\nu_{L_1} \rm{dist}(X_k, X_{k-1}) = \nu_{L_1} \rm{dist}\left({\cal R}_{X_{k-1}}(\alpha_{k-1}\eta_{k-1}), X_{k-1}\right)\\ &\leq&\nu_{L_1}\mu\|\alpha_{k-1}\eta_{k-1}\| \equiv \lambda\alpha_{k-1}\|\eta_{k-1}\|. \end{eqnarray*} $

令$ \lambda = \nu_{L_1}\mu $则命题得证.

在算法$ 1 $中, 对$ \forall k\geq1 $, 有

(3.8) $ \begin{equation} \|{\cal T}_{\alpha_{k-1}\eta_{k-1}} \eta_{k-1}\| = \|{\bf P}_{{\cal R}_{X_{k-1}}(\alpha_{k-1}\eta_{k-1})} (\eta_{k-1})\| = \|{\bf P}_{X_k}(\eta_{k-1})\|\leq\|\eta_{k-1}\|. \end{equation} $

由上式, (2.22)及(2.23)式可知, 对$ \forall k\geq1 $, 有如下不等式成立

(3.9) $ \begin{eqnarray} \beta_k& = &\frac{\langle {\rm{grad}} f(X_k), Y_{k-1}\rangle}{\|{\rm{grad}} f(X_{k-1})\|^2}\leq \frac{\|{\rm{grad}} f(X_k)\| \|Y_{k-1}\| }{\|{\rm{grad}} f(X_{k-1})\|^2}, {}\\ \theta_k & = &\frac{\langle {\rm{grad}} f(X_k), {\cal T}_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})\rangle}{\|{\rm{grad}} f(X_{k-1})\|^2}\leq \frac{\|{\rm{grad}} f(X_k)\| \|{\cal T}_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})\|}{\|{\rm{grad}} f(X_{k-1})\|^2}\\ &\leq&\frac{\|{\rm{grad}} f(X_k)\| \|\eta_{k-1}\|}{\|{\rm{grad}} f(X_{k-1})\|^2}.{} \end{eqnarray} $

由引理3.1知$ {\cal Z} $是紧集, 又因目标函数$ f $是连续可微的且$ {\rm{grad}}f(X) $可作为定义在$ {\rm{St}}(n, p) $与$ {\Bbb R}^{n\times p} $间的一个连续非线性映射, 于是存在常数$ \gamma>0 $使得

(3.10) $ \begin{equation} \|{\rm{grad}} f(X)\|\leq\gamma , \ \ X\in {\cal Z}. \end{equation} $

又因目标函数序列$ \{f(X_k)\} $是单调递减的, 故由算法1中产生的迭代序列$ \{X_k\} $属于紧集$ {\cal Z} $.

引理3.2  若存在常数$ \varepsilon>0 $, 使得对$ \forall k $有

(3.11) $ \begin{equation} \| {\rm{grad}}f(X_k)\|\geq\varepsilon, \end{equation} $

则存在常数$ \Theta>0 $, 使得对$ \forall k $有

(3.12) $ \begin{equation} \|\eta_k\|\leq\Theta. \end{equation} $

证  由命题3.1, 公式(2.20)及(3.8)-(3.12)可知, 存在整数$ k_1>0 $, 使得对$ \forall k>k_1 $有

(3.13) $ \begin{eqnarray} \|\eta_k\| &\leq & \|{\rm{grad}} f(X_k)\|+\|\beta_k{\cal T}_{\alpha_{k-1}\eta_{k-1}}(\eta_{k-1})\|+\|\theta_kY_{k-1}\|{}\\ &\leq &\gamma+\frac{\|{\rm{grad}} f(X_k)\| \|Y_{k-1}\|}{\|{\rm{grad}} f(X_{k-1})\|^2}\|\eta_{k-1}\|+\frac{\|{\rm{grad}} f(X_k)\|\|\eta_{k-1}\| }{\|{\rm{grad}} f(X_{k-1})\|^2}\|Y_{k-1}\|{}\\ &\leq & \gamma+2\frac{\|{\rm{grad}}f(X_k)\| \|Y_{k-1}\|}{\|{\rm{grad}}f(X_{k-1})\|^2}\|\eta_{k-1}\|{}\\ &\leq & \gamma+2\frac{\gamma\lambda\alpha_{k-1}\|\eta_{k-1}\|}{\varepsilon^2}\|\eta_{k-1}\|. \end{eqnarray} $

由$ (3.2) $式可得

$ \lim\limits_{k\rightarrow \infty}\frac{2\gamma\lambda}{\varepsilon^2}\alpha_{k-1}\|\eta_{k-1}\| = 0. $

于是存在常数$ r\in(0, 1) $及整数$ k_2>0 $, 使得对$ \forall k>k_2 $都有

(3.14) $ \begin{equation} \frac{2\gamma\lambda}{\varepsilon^2}\alpha_{k-1}\|\eta_{k-1}\|\leq r. \end{equation} $

再由上式及(3.13)式可知对$ \forall k>k_{*} = \max\{k_1, k_2\} $有

$ \begin{eqnarray*} \|\eta_k\|&\leq&\gamma+r\|\eta_{k-1}\|\leq\gamma(1+r+r^2+, \cdots , r^{k-k_{*}-1})+r^{k-k_{*}}\|\eta_{k_{*}}\|\leq \frac{\gamma}{1-r}+\|\eta_{k_{*}}\|. \end{eqnarray*} $

最后令$ \Theta_1 = \max\{\|\eta_1\|, \|\eta_2\|, \cdots , \|\eta_{k_{*}}\|\} $可知对$ \forall k $, 有

$ \|\eta_k\| \leq \max\{\Theta_1, \frac{\gamma}{1-r}+\|\eta_{k_{*}}\| \}\equiv \Theta. $

证毕.

定理3.1  令$ \{X_k\} $是由算法$ 1 $产生的迭代序列, 则有

$ \lim\limits_{k\rightarrow \infty} \inf \| {\rm{grad}} f(X_k)\| = 0. $

证  反证法. 假设存在常数$ \epsilon>0 $, 使得对$ \forall k $有

(3.15) $ \begin{equation} \|{\rm{grad}} f(X_k)\|\geq\epsilon. \end{equation} $

则由(3.1)式知对$ \forall k\geq1 $有

$ \alpha_k\| {\rm{grad}} f(X_k)\|^2 = -\alpha_k\langle \eta_k, {\rm{grad}} f(X_k)\rangle \leq \alpha_k\|\eta_k\|\cdot \| {\rm{grad}} f(X_k)\|, $

即对$ \forall k\geq1 $有

$ 0\leq \alpha_k\|{\rm{grad}}f(X_k)\|\leq \alpha_k \|\eta_k\|. $

由上式及(3.2)式可知

$ \lim\limits_{k\rightarrow \infty} \alpha_k\| {\rm{grad}}f(X_k)\| = 0. $

因此, 若$ { } \mathop {\lim \inf }\limits_{k \to \infty } \alpha_k>0 $, 则$ { } \mathop {\lim \inf }\limits_{k \to \infty } \|{\rm{grad}}f(X_k)\| = 0 $而这与$ (3.15) $式矛盾.

现假设$ { } \mathop {\lim \inf }\limits_{k \to \infty } \alpha_k = 0 $, 可先设

(3.16) $ \begin{equation} \lim\limits_{k\rightarrow \infty}\alpha_k = 0. \end{equation} $

根据算法1的(3.24)式知, 当$ k $足够大时有$ \rho^{-1}\alpha_k $不满足(2.24)式, 即

(3.17) $ \begin{equation} f\left({\cal R}_{X_k}(\rho^{-1}\alpha_k\eta_k)\right)-f(X_k)\geq -\delta\rho^{-2}\alpha_k^2\|\eta_k\|^2. \end{equation} $

由目标函数$ f $及收缩算子$ {\cal R} $都是连续可微的可知回调函数$ \widehat{f} $也是连续可微的, 于是存在常数$ \delta_2>0, \nu_{L_2}>0 $使得

(3.18) $ \begin{equation} \|{\rm{grad}}\widehat{f}_X(\eta_x)-{\rm{grad}} \widehat{f}_X(\eta_x)\|\leq \nu_{L_2}\| \eta_x-\eta_x\|, \; \; \; \; X\in{\cal Z}, \; \eta_x, \eta_x \in {\bf T}_X{\rm{St}}(n, p), \end{equation} $

这里$ \|\eta_x-\eta_x\|<\delta_2 $. 再由中值定理, 引理3.2, 公式(3.1) 及(3.17)知, 存在常数$ \varphi_k\in(0, 1) $, 使得对$ k $足够大有

(3.19) $ \begin{eqnarray} f\left({\cal R}_{X_k}(\rho^{-1}\alpha_k\eta_k)\right)-f(X_k) & = &f\left({\cal R}_{X_k}(\rho^{-1}\alpha_k\eta_k)\right)-f\left({\cal R}_{X_k}(0_{X_k})\right){}\\ & = &\widehat{f}_{X_k}(\rho^{-1}\alpha_k\eta_k)-\widehat{f}_{X_k}(0_{X_k}){}\\ & = &\rho^{-1}\alpha_k \langle {\rm{grad}}\widehat{f}_{X_k}(\varphi_k\rho^{-1}\alpha_k\eta_k), \eta_k\rangle{}\\ & = &\rho^{-1}\alpha_k \langle {\rm{grad}}\widehat{f}_{X_k}(0_{X_k}), \eta_k\rangle+\rho^{-1}\alpha_k \langle {\rm{grad}}\widehat{f}_{X_k}(\varphi_k\rho^{-1}\alpha_k\eta_k), \eta_k\rangle{}\\ &&-\rho^{-1}\alpha_k \langle {\rm{grad}}\widehat{f}_{X_k}(0_{X_k}), \eta_k\rangle{}\\ & \leq& \rho^{-1}\alpha_k \langle {\rm{grad}}\widehat{f}_{X_k}(0_{X_k}), \eta_k\rangle+\nu_{L_2}\varphi_k^2\rho^{-2}\alpha_k^2\|\eta_k\|^2{}\\ &\leq& -\rho^{-1}\alpha_k \| {\rm{grad}} f(X_k)\|^2 +\nu_{L_2}\rho^{-2}\alpha_k^2\|\eta_k\|^2. \end{eqnarray} $

将上式与$ (3.17) $式联立知, 对$ \forall k $足够大时有

$ \| {\rm{grad}} f(X_k)\|^2\leq(\nu_{L_2}+\delta)\rho^{-1}\alpha_k\|\eta_k\|^2. $

最后由引理3.2知搜索方向序列$ \{\eta_k\} $是有界的, 结合(3.16)式便有

$ \lim\limits_{k\rightarrow \infty} \|{\rm{grad}} f(X_k)\| = 0. $

这也与$ (3.15) $式矛盾, 故可知定理3.1成立, 即产生的迭代序列$ \{X_k\} $是收敛的.

本节利用数值实验和数值比较验证算法1对于求解问题(1.1)是高效可行的. 所有的数据均是通过MATLAB(R2014a), Inte$ l^{\circledR} $ Core i5处理器, 3.20GHz的PC中获得. 为说明算法1的高效性, 本节将与两类不可行算法-交替方向法^[16]和分裂正交法^[17], 以及三类可行算法- Majorization算法^[15]以及两类基于Stiefel流形的梯度下降法进行数值比较.

首先给出各比较算法的基本迭代思路和迭代格式主要步骤.

交替方向法(ADMM)^[16]. 引入辅助矩阵$ V\in {{\Bbb R}}^{n\times p} $, 问题(1.1) 可转化为

(4.1) $ \begin{equation} \min\limits_{X, V\in {{\Bbb R}}^{n\times p}}\; \frac{1}{2}\big\| \sum\limits_{i = 1}^{N}A_iXB_i-C\big\|^{2} \; \; \; {\rm{subject\ to}}\; X-V = 0, \; \; V^TV-I_p = 0. \end{equation} $

问题(4.1)的拉格朗日函数为

$ \begin{eqnarray*} {\cal L}_\alpha(X, V, \Lambda) = \frac{1}{2}\big\| \sum\limits_{i = 1}^{N}A_iXB_i-C\big\|^{2} -\langle \Lambda, X-V\rangle+\frac{\alpha}{2}\|X-V\|^2, \end{eqnarray*} $

其中$ \Lambda\in {{\Bbb R}}^{n\times p} $为拉格朗日乘子矩阵, $ \alpha $为罚参数. ADMM的迭代思想是交替极小化参数矩阵$ V $和$ X $. 其主要迭代格式如下

(4.2) $ \begin{eqnarray} X_{k+1}:& = &\mathop{\rm argmin}\limits_ {X\in {{\Bbb R}}^{n\times p}} {\cal L}_\alpha(X, V_k, \Lambda_k){}\\ & = &\mathop{\rm argmin}\limits_{V\in {{\Bbb R}}^{n\times p}}\frac{1}{2}\big\| \sum\limits_{i = 1}^{N}A_iXB_i-C\big\|^{2} +\frac{\alpha}{2}\|X-V_k-\frac{1}{\alpha}\Lambda_k\|^2, \\ V_{k+1}:& = &\mathop{\rm argmin}\limits_{V:V^TV = I_p} {\cal L}_\alpha(X_{k+1}, V, \Lambda_k) = \mathop{\rm argmin}\limits_{V:V^TV = I_p} \|V-(X_{k+1}-\frac{1}{\alpha}\Lambda_k)\|^2, {}\\ \Lambda_{k+1}:& = &\Lambda_{k}-(X_{k+1}-V_{k+1}).{} \end{eqnarray} $

(4.2)式中的$ X $子问题是一个无约束矩阵优化问题, 可采用带BB步长的梯度下降法求解. (4.2)式中的$ V $子问题实质上是Stiefel流形上的正交投影问题, 其解析解可表示为$ V_{k+1} = UI_{n\times p}W^T $, 其中$ U\in {{\Bbb R}}^{n\times n} $, $ W\in {{\Bbb R}}^{p\times p} $为奇异值分解$ X_{k+1}-\Lambda_k = UDW^T $中的正交奇异因子, $ D\in {{\Bbb R}}^{n\times p} $为奇异值分解中的对角矩阵.

正交分裂法(SOC)^[17]. 正交分裂法的迭代思想与交替方向法相似, 即引入变量分裂正交约束, 然后交替极小两个变量矩阵. 引入变量矩阵$ P = X\in {{\Bbb R}}^{n\times p} $分裂正交约束条件, 问题(1.1)可转化为

(4.3) $ \begin{equation} \min\limits_{X, P\in {{\Bbb R}}^{n\times p}}\; \frac{1}{2}\big\| \sum\limits_{i = 1}^{N}A_iXB_i-C\big\|^{2}\; \; \; {\rm{subject\ to}}\; X = P, \; \; P^TP = I_p. \end{equation} $

因为(4.3)式中第一个约束条件是线性约束, 故可以采用Bregman迭代处理, 即迭代求解

(4.4) $ \begin{eqnarray} &&{ } (X_k, P_k) = \mathop{\rm argmin}\limits_ {X, P\in {{\Bbb R}}^{n\times p}}\frac{1}{2}\big\| \sum\limits_{i = 1}^{N}A_iXB_i-C\big\|^{2}+\frac{\alpha}{2}\|X-P+\Lambda_{k-1}\|^2, \; \; \rm{subject\ to}\; P^TP = I_p, {}\\ && \Lambda_k = \Lambda_{k-1}+X_k-P_k. \end{eqnarray} $

类似于交替方向乘子法和分裂Bregman迭代, (4.4)式的第一个极值问题采用交替极小变量矩阵$ X $和$ P $的方式求解, 即得到如下迭代格式

(4.5) $ \begin{eqnarray} && X_{k}: = \mathop{\rm argmin}\limits_{X\in {{\Bbb R}}^{n\times p}} \frac{1}{2}\big\| \sum\limits_{i = 1}^{N}A_iXB_i-C\big\|^{2}+\frac{\alpha}{2}\|X-P_{k-1}+\Lambda_{k-1}\|^2, {}\\ &&P_{k}: = \mathop{\rm argmin}\limits_{P: P^TP = I_p}\frac{\alpha}{2}\|P-(X_{k}+\Lambda_{k-1})\|^2, \\ &&\Lambda_{k}: = \Lambda_{k-1}-X_{k}+P_{k}.{} \end{eqnarray} $

Majorization算法(Major)^[15]. Majorization算法的基本迭代思想是对于$ {{f}}(X) $, 构造一个约束极值点容易求的辅助函数$ {g}(X) $且对$ \forall X\in {\rm{St}}(n, p) $满足$ {g}(X)\geq f(X) $. 每一次迭代都利用$ {g}(X) $的极小值点来得到$ {{f}}(X) $的可行下降点. 下面针对问题(1.1)的目标函数$ {{f}}(X) $, 简述如何构造有效辅助函数$ g(X) $. 任意给定初始点$ \widetilde{X}\in {\rm{St}}(n, p) $, 定义$ E = X-\widetilde{X} $, 则有

(4.6) $ \begin{equation} {f}(X) = {f}(\widetilde{X})+\langle \sum\limits_{i = 1}^N A_i^T\left(\sum\limits_{i = 1}^N A_i\widetilde{X}B_i-C\right)B_i^T, \; E \rangle+\frac{1}{2}\big\|\sum\limits_{i = 1}^N A_iEB_i\big\|^2. \end{equation} $

又因为下面不等式成立

(4.7) $ \begin{equation} {\rm{trace}}(B_j^TE^TA_j^TA_iEB_i) = {\rm{trace}}(E^TA_j^TA_iEB_i B_j^T)\leq \lambda_{ij}\|E\|^2, \; \; i\neq j \end{equation} $

和

(4.8) $ \begin{equation} \|A_iEB_i\|^2 = {\rm{trace}}(E^TA_i^TA_iEB_i B_i^T)\leq \lambda_{i}\|E\|^2, \end{equation} $

其中$ \lambda_{ij} $为$ \frac{1}{2}(A_i^TA_j) \otimes (B_iB_j^T)+\frac{1}{2}(A_j^TA_i) \otimes (B_jB_i^T) $的最大特征值, $ \lambda_{i} $为$ (A_i^TA_i) \otimes (B_iB_i^T) $的最大特征值. 结合(3.6)-(3.8) 式有

(4.9) $ \begin{equation} {f}(X)\leq {f}(\widetilde{X})+\langle \sum\limits_{i = 1}^N A_i^T\left(\sum\limits_{i = 1}^N A_i\widetilde{X}B_i-C\right)B_i^T, \; E \rangle+\left(\frac{1}{2}\sum\limits_{i = 1}^N\lambda_i+\sum\limits_{i<j}\lambda_{ij}\right)\|E\|^2. \end{equation} $

将$ E = X-\widetilde{X} $代入(4.9)式右端表达式, 并记得到的关于$ X $的函数为$ g(X) $, 即

$ \begin{eqnarray} {g}(X)&: = &{f}(\widetilde{X})+\big\langle \sum\limits_{i = 1}^N A_i^T\big(\sum\limits_{i = 1}^N A_i\widetilde{X}B_i-C\big)B_i^T, \; X-\widetilde{X} \big\rangle+\big(\frac{1}{2}\sum\limits_{i = 1}^N\lambda_i+\sum\limits_{i<j}\lambda_{ij}\big)\|X-\widetilde{X}\|^2\\ &: = &c+\big\langle \sum\limits_{i = 1}^N A_i^T\big(\sum\limits_{i = 1}^N A_i\widetilde{X}B_i-C\big)B_i^T-\big(\sum\limits_{i = 1}^N\lambda_i+2\sum\limits_{i<j}\lambda_{ij}\big)\widetilde{X}, \; X\big\rangle, \end{eqnarray} $

其中$ c $为与$ X $无关的常量, 且上式中用到$ \|X\|^2 = p $, $ \forall X\in {\rm{St}}(n, p) $. 易知对$ \forall X\in {\rm{St}}(n, p) $有$ {g}(X) $$ \geq f(X) $, 且有$ {g}(\widetilde{X}) = f(\widetilde{X}) $. 若记

$ M^T = \sum\limits_{i = 1}^N A_i^T\left(\sum\limits_{i = 1}^N A_i\widetilde{X}B_i-C\right)B_i^T-\left(\sum\limits_{i = 1}^N\lambda_i+2\sum\limits_{i<j}\lambda_{ij}\right)\widetilde{X}, $

则辅助函数$ g(X) $可表示为

$ g(X) = c+{\rm{trace}}(MX). $

现若$ X_*\in {\rm{St}}(n, p) $是$ {g}(X) $在约束集合$ {\rm{St}}(n, p) $中的极小值点, 则对$ \forall X\in {\rm{St}}(n, p) $, $ {g}(X_*)\leq g(X) $. 特别地, 有$ {g}(X_*)\leq g(\widetilde{X}) $. 同时又有$ g(\widetilde{X}) = f(\widetilde{X}) $. 因此对于目标函数$ f(X) $和任意初始点$ \widetilde{X}\in {\rm{St}}(n, p) $, 有$ f(X_*)\leq f(\widetilde{X}) $, 也即找到了一个使得目标函数值下降的可行点$ X_* $. 重复该过程直到目标函数$ {f}(X) $不再下降, 则可得到$ {f}(X) $的局部极小值点. 同时, $ {g}(X) $在约束集合$ {\rm{St}}(n, p) $中的极小值点容易求得. 事实上, 令$ -M\in {{\Bbb R}}^{p\times n} $的简约版奇异值分解为$ -M = P\Lambda Q^T $, 其中$ P\in {{\Bbb R}}^{p\times p}, Q\in {{\Bbb R}}^{n\times p} $, 且$ P^TP = PP^T = Q^TQ = I_p $, $ \Lambda\in {{\Bbb R}}^{p\times p} $为$ M $的奇异值组成的对角阵. 则可知极值问题$ \min_{X^TX = I_p}{g}(X) $的极小值点可表示为$ X_* = QP^T $.

基于BB步长非单调搜索技术的梯度下降法(OptStiefel)^[27]. 文献[27]提出的基于Cayley变换和结合BB步长非单调搜索的梯度下降法, 其迭代更新格式为

$ \begin{eqnarray*} X_{k+1} = \left(I_n+\frac{\alpha_k}{2}\tilde{A}_k\right)^{-1}\left(I_n-\frac{\alpha_k}{2}\tilde{A}_k\right)X_k. \end{eqnarray*} $

其中$ \alpha_k $表示搜索步长, 反对称矩阵$ \tilde{A}_k = G_kX_k^T-X_kG_k^T $, $ G_k $为目标函数$ f(X) $关于$ X\in {{\Bbb R}}^{n\times p} $的欧氏梯度(2.14)在$ X_k $的函数值. 当$ p\ll n $时, 为降低算法复杂度, 结合Sherman-Morrison-Woodbury公式, 迭代格式简化为

$ \begin{eqnarray*} X_{k+1} = X_k - \alpha_k U_k\left(I_{2p} + \frac{\alpha_k}{2}V_k^TU_k\right)^{-1}V_k^TX_k, \end{eqnarray*} $

其中$ U_k = [G_k, X_k] $, $ V_k = [X_k, -G_k] $. 对于给定参数$ \delta, \rho, \varrho\in (0, 1) $以及搜索方向$ \eta_k $, 步长$ {\alpha}_k = \overline{\alpha}_k\rho^s $满足如下非单调Armijo法则

(4.10) $ \begin{equation} f({\cal R}_{x_k}(\overline{\alpha}_k\rho^s\eta_k))\leq C_k+\delta\overline{\alpha}_k\rho^s \langle \nabla f(x_k), \eta_k\rangle. \end{equation} $

其中$ \overline{\alpha}_k $表示第$ k $的初始步长, 且由如下Barzilai-Borwein(BB)方式确定

$ \overline{\alpha}_k^{1} = \frac{{\rm{trace}}(S_{k-1}^TS_{k-1})}{|{\rm{trace}}(Y_{k-1}^TS_{k-1})|}, \ \ \ \rm{或}\ \ \overline{\alpha}_k^{2} = \frac{|{\rm{trace}}(S_{k-1}^TY_{k-1})|}{{\rm{trace}}(Y_{k-1}^TY_{k-1})}, $

其中$ S_{k-1} = X_k-X_{k_1} $, $ Y_{k-1} = \nabla f(X_k)-\nabla f(X_{k-1}) $. (3.10) 中参考函数$ C_{k+1} $为$ C_k $和$ f(X_k) $的凸组合, 即$ C_{k+1} = (\varrho Q_kC_k+f(X_{k+1}))/Q_{k+1} $, 而$ C_0 = f(X_0) $, $ Q_{k+1} = \varrho Q_k+1 $和$ Q_0 = 1 $.

基于BB步长非单调搜索技术的混合方向下降法(Mixed)^[28]. 文献[28] 提出的基于混合方向和非单调线搜索的梯度类方法. 其迭代更新方案为

$ X_{k+1} = {\bf \pi}({X_k-\alpha_kH_k}), $

其中$ \pi $表示$ {{\Bbb R}}^{n\times p}\rightarrow {\rm{St}}(n, p) $的正交投影算子, $ H_k $为位于$ X_k $切平面上两个方向$ G_k-X_kG_kX_k^T $和$ I_n-X_kX_k^T $的线性组合, 即$ H_k = \theta_1(G_k-X_kG_kX_k^T)+\theta_2 (I_n-X_kX_k^T) $. 步长$ \alpha_k $采用文献[27]中结合BB步长的非单调线搜索方式确定(4.10)式.

(1) 测试模型. 测试问题(1.1)的如下具体形式

$ \min\limits_{X\in {\rm{St}}(n, p)} \frac{1}{2}\|A_1XB_1+A_2XB_2\|^2. $

(2) 终止标准. 算法的终止标准参见文献[27]: $ \|{\rm{grad}}f(X_{k})\|\leq \varepsilon $或$ X_{\rm diff}\leq \varepsilon_x $且$ f_{\rm diff}\leq \varepsilon_f $, 其中

$ X_{\rm diff} = \frac{\|X_{k}-X_{k-1}\|}{\sqrt{n}}, \ \ \ {f}_{\rm diff} = \frac{|f(X_{k})-f(X_{k-1})|}{|f(X_{k-1})|+1}. $

(3) 参数选定. 算法1中参数$ \varepsilon = 10^{-3} $, $ \varepsilon_x = 10^{-6} $, $ \varepsilon_f = 10^{-12} $, $ \delta = 10^{-4} $, $ \rho = 0.2 $. 最大迭代步设为$ K = 20000. $ ADM算法中初始矩阵$ X_0 $和$ \Lambda_0 $, 和SOC算法中的初始矩阵$ P_0 $和$ \Lambda_0 $均取为合适维数的零矩阵. 罚参数$ \alpha $取为$ \alpha = 0.1 $, 且随系统维数自适应调整. Major算法中初始迭代矩阵取为$ X_0 = {\rm orth(rand}(n, p)) $, 其中rand$ (n, p) $表示$ n\times p $随机生成矩阵, orth$ (\cdot ) $表示Matlab中的列正交指令. OptStiefel和Mixed算法均取为所参考文献中给定的参数标准.

例4.1  按Matlab记号的形式给定系数矩阵$ A_1, A_2, B_1, B_2, C $如下

$ {A_1} = {\rm rand}({l, n}), \; {B_1} = {\rm rand}({p, s}), \; {A_2} = {\rm rand}({l, n}), \; {B_2} = {\rm rand}({p, s}), \; C = {\rm ones}({l, s}), $

其中ones$ ({l, s}) $表示$ l\times s $元素全为1的矩阵. 初始迭代矩阵$ X_0 = {\rm orth(rand}(n, p)) $.

对于不同维数参数产生的随机矩阵$ A_1, B_1, A_2, B_2 $和$ C $, 表 1给出了算法1求解问题(1.1) 的数值比较结果, 其中"CPU" 表示迭代时间, 单位(s), "IT" 表示迭代步, $ { X_{\rm feasi}} = \|X_k^TX_k-I_p\| $表示当前迭代点$ X_k $的可行度. 图 1给出了不同系数维数下, 各算法的黎曼梯度范数$ \| {\rm{grad}} f(X_k)\| $随迭代步$ k $的收敛曲线图. 结合表 1的数据结果和图 1的收敛曲线图, 我们分析得如下结论.

(1) 对于求解问题(1.1), 与其它几类可行算法相比, 算法1在迭代步和迭代时间上均有较大的优势. 与两类不可行算法相比, 算法1的迭代优势则更明显. 虽然算法1的总体迭代步较多, 但每步迭代成本低, 故达到相同精度下的总体迭代时间要少很多. 特别地, 文献[27]中提出的基于BB步长非单调搜索的梯度下降法(OptStiefel)是近年来求解流形优化问题的经典算法, 而对于问题模型(1.1), 本文算法1有与OptStiefel算法相当的迭代效率.

(2) Mixed算法因在$ X_k $的迭代搜索方向取为$ X_k $的切平面上两个方向的线性组合, 但最优组合系数无法确定. 本实验中均取为固定组合系数, 并未自适应调整, 故该算法总体迭代效率较差. ADM和SOC算法虽然迭代步较少, 但相对迭代时间长, 因为两类算法都需内迭代, 即每次迭代都需求解(4.2) 式中$ V $子问题和(4.5) 式中$ X $子问题所示的无约束矩阵优化问题. 另外, ADM和SOC算法中指标"$ { X_{\rm feasi}} $" 给出的是$ X $子问题序列所对应的可行度. 因为ADM和SOC算法都是不可行算法, 故$ X $子问题迭代更新矩阵不一定满足列正交约束条件. Majorization算法对于求解问题(1.1) 虽可行, 且能得到使目标函数单调下降的可行点序列, 但收敛速度相对较慢. 并且因为该算法需要用到Kronecker乘积后矩阵的最大特征值, 故不适用于大规模问题.

本文从Stiefel流形的角度设计黎曼MPRP共轭梯度法求解列正交约束下广义Sylvester方程极小化问题. 该算法可视为欧氏空间上MPRP共轭梯度法在Stiefel流形上的推广. 仿照欧氏空间上MPRP共轭梯度法的全局收敛性分析, 基于Stiefel流形的基本性质, 本文给出了该黎曼MPRP共轭梯度法的完整收敛性分析. 同时数值实验部分给出数值比较结果, 说明本文算法在迭代效率上与几类不可行算法和几类基于Stiefel流形的梯度下降法有较大的优势.