给定复矩阵$ A\in {\mathbb C}^{m\times n} $ ($ m\geq n $), 其奇异值分解(SVD) 为$ A=U\Sigma V^H=\sum\limits_{i=1}^{n}\sigma_iu_iv_i^H, $其中$ U=(u_1, u_2, \cdots , u_m)\in {\mathbb C}^{m\times m} $, $ V=(v_1, v_2, \cdots , v_n)\in {\mathbb C}^{n\times n} $均为酉阵, 即满足$ U^HU=UU^H=I_m $, $ V^HV=VV^H=I_n $, 其中$ I_n $为$ n $阶单位矩阵, $ u_i $和$ v_i $分别称为左、右奇异向量. $ \Sigma=(\Sigma_1^T, O)^T $, $ \Sigma_1=\mbox{diag}(\sigma_1, \sigma_2, \cdots , \sigma_n) $, $ \sigma_i $ ($ 1\leq i\leq n $) 为矩阵$ A $的奇异值, 且$ \sigma_1\geq \sigma_2\geq \cdots \geq \sigma_n $, $ O $表示适当维数的零矩阵. 从直观上讲, $ U $和$ V^H $可视为旋转操作, $ \Sigma $可视为缩放操作. 因此奇异值分解的含义可表述为, 若将矩阵看做一个变换, 那么任何这样的变换可以看做是两个旋转和一个缩放变换的复合. 它也可以通过对$ \Sigma $的调整实现降维, 通过$ U $和$ V^H $在高维和低维间转换. 截断奇异值分解(TSVD) 是SVD的变形, 只计算用户指定的最大的$ p $个奇异值. TSVD与SVD不同的是它可以产生一个指定维度的分解矩阵. 如对$ n\times n $阶矩阵, 通过SVD分解后仍然是一个$ n\times n $阶矩阵, 而TSVD可以生成指定维度的矩阵, 这样就可实现降维. 假设$ \sigma_p $是要保留的最小的非零奇异值, 则对于任意$ p<\mbox{rank}(A) $, 取$ A_p=\sum\limits_{i=1}^{p}\sigma_iu_iv_i^H $, 则矩阵$ A_p $是秩为$ p $的矩阵中最接近原矩阵$ A $的一个较低维矩阵.

TSVD经常被用在病态问题^{[1, 2]}的解决上. 在病态模型问题中, 模型的病态性主要体现在系数矩阵的小奇异值对参数及其方差的放大. 截断奇异值法的基本思想是将这些小的奇异值截掉并重构系数阵以削弱模型的病态性. 如在生物发光断层成像(BLT)^[3]中, 由于BLT重建问题中系数矩阵$ A $的病态性, 随着其奇异值逐渐趋近于零, 数据中的微小噪声将被放大, 所求得的解将会远远偏离问题的真实解. 考虑到在实际测量中, 噪声是不可避免的, 为了得到重建问题的稳定解, 可将TSVD方法作用于模型中的系数矩阵$ A $, 通过改造系数矩阵$ A $, 将容易造成解不稳定的较小的奇异值直接去除. 更多TSVD方法在微波成像, 核磁共振波谱分析, 压缩光谱图像分析和信号处理等领域中的应用见文献[4–5].

对于复矩阵$ A \in {\mathbb C}^{m\times n} $$ (m\geq n) $, 其截断奇异值分解就是寻找$ A $的$ p $个最大奇异值和与之相关的奇异向量. 与之密切相关的是求解如下矩阵优化问题^[6–9]

(1.1) $ \begin{equation} \begin{array}{rl} \mathop{{\rm maximize}}\limits_{U, V}& {{{\rm{Tr}}}}(U^HAV\Theta), \\ \mbox{subject to}& U\in {\mathbb C}^{m\times p}, V\in {\mathbb C}^{n\times p}, U^HU=V^HV=I_p, \end{array} \end{equation} $

其中$ \Theta=\mbox{diag}(\mu_1, \cdots, \mu_p) $, $ \mu_1> \cdots > \mu_p>0 $, $ 1\leq p\leq \mbox{rank}(A) $, $ {{{\rm{Tr}}}}(\cdot) $表示复数域内的内积, 即$ \forall B, C\in {\mathbb C}^{m\times n} $, 其定义为

$ {{{\rm{Tr}}}}(B^HC)={\rm Re}({{{\rm{tr}}}}(B^HC))={{{\rm{tr}}}}\left({\rm Re}(B)^T{\rm Re}(C)\right) + {{{\rm{tr}}}}\left({\rm Im}(B)^T{\rm Im}(C)\right), $

其中$ {{{\rm{tr}}}}(\cdot) $表示括号内矩阵的迹, $ {{{\rm{Re}}}}(\cdot) $和$ {{{\rm{Im}}}}(\cdot) $分别表示括号内复数取实部和虚部. 由此内积诱导出的范数定义为$ \|B\|=\sqrt{{{{\rm{Tr}}}}(B^HB)}=({{{\rm{Re}}}}({{{\rm{tr}}}}(B^HB)))^{{1}/{2}} $. 针对问题(1.1), Sato等^[7]证明若$ (U_*, V_*) $是问题(1.1)的最优解, 则$ U_* $的第$ j $列和$ V_* $的第$ j $列分别为$ A $的第$ j $个最大奇异值对应的左、右奇异向量. $ A $的前$ p $个最大奇异值$ \sigma_1\geq \cdots \geq \sigma_p $可通过$ U_*^HAV_*=\mbox{diag}(\sigma_1, \cdots , \sigma_p) $求得. 注意到约束矩阵集合$ \{V\in {\mathbb C}^{n\times p}: V^HV=I_p, n\geq p\} $通常称为复Stiefel流形, 记为$ {{{\rm{St}}}}(n, p, {\mathbb C}) $, 故问题(1.1)的可行集可视为复乘积流形$ {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $. 因此复矩阵$ A\in {\mathbb C}^{m\times n} $的截断奇异值分解问题可等价转化为如下复乘积流形下的黎曼优化问题:

(1.2) $ \begin{equation} \begin{array}{rl} \mbox{minimize}& F(U, V):=-{{{\rm{Tr}}}}(U^HAV\Theta), \\ \mbox{subject to}& (U, V)\in {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}). \end{array} \end{equation} $

其中$ \Theta=\mbox{diag}(\mu_1, \cdots, \mu_p) $, $ \mu_1> \cdots > \mu_p>0 $, $ 1\leq p\leq \mbox{rank}(A) $.

针对问题(1.2), Sato等^{[8, 9]}从黎曼优化角度设计了若干优化算法求解. 基于$ n\times n $阶复矩阵和$ 2n\times 2n $阶实分块矩阵的如下同构关系

$ D=E+{\rm i} F\in {\mathbb C}^{n\times n}, E, F\in {{\mathbb R}}^{n\times n} \simeq \widetilde{D} =\left[\begin{array}{ccccc}{E}{\quad} & {F} \\ {-F}{\quad} & {E}\end{array}\right]\in {{\mathbb R}}^{2n\times 2n}, $

以及任意$ 2n\times 2n $阶实矩阵$ \widehat{D} $具有类似于$ \widetilde{D} $的分块结构当且仅当

$ \widehat{D}J_n=J_n\widehat{D}, J_n:=\left[\begin{array}{ccccc}{O}{\quad} & {I_n} \\ {-I_n}{\quad} & {O}\end{array}\right]. $

对正整数$ p, q $定义拟辛集

$ {\mathcal SP}(p, q)=\big\{\widetilde{D}\in {{\mathbb R}} ^{2p\times 2q}|\widetilde{D}J_p=J_q\widetilde{D}\big\}\subset {{\mathbb R}} ^{2p\times 2q}, $

和同构映射

$ \phi^{(n, p)}: {\mathbb C}^{n\times p}\rightarrow {\mathcal SP}(n, p); \widetilde{D}=\phi^{(n, p)}(D)=\phi^{(n, p)}(E+{\rm i} F)= \left[\begin{array}{ccccc}{E} {\quad} & {F} \\ {-F} {\quad} & {E}\end{array}\right], $

则复Stiefel流形$ {{{\rm{St}}}}(n, p, {\mathbb C}) $同构于如下$ 2n\times 2p $阶实Stiefel流形和拟辛集的交集, 即

$ {{{\rm{St}}}}(n, p, {\mathbb C})\simeq \mbox{Stp}(n, p):={{{\rm{St}}}}(2n, 2p, {{\mathbb R}})\cap {\mathcal SP}(n, p), $

进而问题(1.2)可转化为如下实数域上乘积流形约束黎曼优化问题

(1.3) $ \begin{equation} \begin{array}{rl} \mbox{minimize}& -\mbox{tr}(\widetilde{U}^T\widetilde{A}\widetilde{V}\widetilde{\Theta}), \\ \mbox{subject to}& (\tilde{U}, \tilde{V})\in \mbox{Stp}(m, p)\times \mbox{Stp}(n, p), \end{array} \end{equation} $

其中$ \tilde{\Theta}=\left[\begin{array}{ccccc}{\Theta} {\quad} & {O} \\ {O}{\quad} & {\Theta}\end{array}\right]. $针对问题(1.3), Sato^[8]设计了基于尺度向量转移算子的黎曼共轭梯度法求解. 值得说明的是, 由于在$ \mbox{Stp}(n, p) $上计算量是$ {{{\rm{St}}}}(n, p, {\mathbb C}) $上的两倍, 且因为黎曼共轭梯度法只用到了目标函数的一阶导数信息, 故针对问题(1.3)设计的黎曼共轭梯度法其收敛速度和迭代解的精度都不太理想. Sato等^[9]针对问题(1.3) 进一步设计了黎曼牛顿法求解, 但针对导出的黎曼牛顿方程只给出了$ p=1 $时的求解方案. 若问题(1.1)退化到实数域, 即问题(1.2)中可行集为实数域上的Stiefel乘积流形$ {{{\rm{St}}}}(m, p, {{\mathbb R}}) $$ \times {{{\rm{St}}}}(n, p, {{\mathbb R}}) $, Sato等^[7]设计了基于乘积流形的梯度下降法、共轭梯度法及牛顿法, 但同样针对黎曼牛顿方程只给出了$ p=1 $时的求解方案. 针对实数域上问题(1.2) 黎曼牛顿法导出的复杂黎曼牛顿方程, 从降低系统维数和简化计算入手, Aihara和Sato^[10]进一步设计了有效求解算法, 即基于Kroncker积和反对称矩阵拉直算子将牛顿方程转化为易于求解的标准对称线性方程组, 并给出无需系数矩阵显示形式的迭代方案, 进而利用Krylov子空间方法求解线性方程组无需系数矩阵显示形式的特点设计相应的Krylov子空间方法求解. 值得注意的是, 因为没有给出系数矩阵的显示形式, 故无法区分由牛顿法收敛得到的驻点是局部最小值、局部最大值或鞍点. 另外, 因为牛顿法的局部收敛特性, 其对初值选取比较敏感, 当初始矩阵选取不好时, 算法可能不收敛. 特别地, Xu^[11]研究了如下一类矩阵优化模型

(1.4) $ \begin{equation} \min\limits_{U_1, \cdots, U_m\in {\mathbb U}_n}\Big| \text{Tr}\Big(cI_n+\prod\limits_{j=1}^mA_jU_j\Big)\Big|, \max\limits_{U_1, \cdots, U_m\in {\mathbb U}_n}\Big| \text{Tr}\Big(cI_n+\prod\limits_{j=1}^mA_jU_j\Big)\Big|, \end{equation} $

其中$ c\in {{\mathbb R}} $, $ A_1, \cdots , A_m $为$ n\times n $复矩阵, $ {\mathbb U}_n $为全体$ n\times n $正交阵构成的集合, 即$ {{{\rm{St}}}}(n, n, {\mathbb C}) $. 文献[11]基于Jocobi公式、矩阵指数和矩阵函数等给出了问题(1.4)解的具体解析表达式, 并给出了问题模型及所提算法在机械系统测试信号与航空发动机故障诊断方面的具体应用. 若令$ c=0, m=2 $和$ A_2=\Theta $, 则问题(1.4)可退化为问题(1.2).

针对以上分析, 本文在文献[7–9]的基础上进一步研究针对复数域上问题(1.2) 的基于目标函数二阶导数信息的黎曼牛顿求解算法. 对于黎曼牛顿法, 有效求解导出的牛顿方程和设计有效初值选取方案是关键. 针对牛顿方程, 本文基于乘积流形上切空间的等价形式, 写出黎曼海塞在切空间上的变换矩阵, 进而基于克罗内克积和复数域上的若干矩阵拉直算子, 将复数域上的黎曼牛顿方程转化为实数域上标准对称线性方程组, 进而利用相关的Krylov子空间方法求解; 针对初值选取, 本文利用黎曼优化中具有单步迭代成本低且存储量少等优点的经典一阶梯度算法OptStiefelGBB^[12]来选取迭代初值, 进而得到一类黎曼混合牛顿法. 充分的数值实验和数值比较说明, 本文针对问题(1.2) 设计的基于黎曼牛顿法的混合算法是切实可行的, 其中数值比较包括黎曼优化工具箱Manopt^[19]中的原始信赖域算法和黎曼BFGS算法, 以及适用于问题(1.2) 的具有不同形式的若干黎曼梯度类算法.

本文余下内容组织如下: 第二节给出复乘积流形$ {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $的若干几何性质, 问题(1.2)中目标函数的对应的黎曼梯度和黎曼海塞的具体形式以及相关拉直算子的定义及性质; 第三节给出黎曼混合牛顿法的迭代框架; 第四节给出黎曼混合信赖域算法的迭代框架; 第五节给出数值实验和数值比较; 最后第六节给出文章总结.

黎曼优化是指黎曼流形$ {\mathcal M} $上的优化, 其一般迭代格式为^[13]

(2.1) $ \begin{equation} x_{k+1}={\mathcal R}_{x_k}(t_k\eta_k), \end{equation} $

其中$ t_k $表示迭代步长, $ \eta_k $为位于点$ {x_k} $对应切平面$ {\bf{T}}_{x_k}{\mathcal M} $上的搜索方向, $ {\mathcal R} $表示$ {\bf{T}}_{x_k}{\mathcal M} $到流形$ {\mathcal M} $的收缩算子, 其定义见文献[13]. 对于复Stiefel流形, 常用的收缩算子有基于指数映射的收缩算子, 基于QR分解的收缩算子, 基于极分解的收缩算子和基于Cayley变换的收缩算子^[12–14].

下面简要给出复乘积流形$ {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $上的若干几何性质. 关于实Stiefel流形$ {{{\rm{St}}}}(m, p, {{\mathbb R}}) $的相关性质, 读者可参考文献[13]. 复Stiefel流形$ {{{\rm{St}}}}(m, p, {\mathbb C}) $和$ {{{\rm{St}}}}(n, p, {\mathbb C}) $的维数分别为$ 2pm-p^2 $和$ 2pn-p^2 $, 故

(2.2) $ \begin{eqnarray} {\rm dim} \big({\rm St}(m, p, {\mathbb C})\times {\rm St}(n, p, {\mathbb C})\big):&=&{\rm dim}\ {\rm St}(m, p, {\mathbb C})+{\rm dim}\ {\rm St}(n, p, {\mathbb C}){}\\ &=&2p(m+n-p). \end{eqnarray} $

复乘积流形$ {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $在点$ (U, V) $处的切空间可表示为^[13]

$ \begin{eqnarray*} && {\bf{T}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C})\simeq {\bf{T}}_{U}{{{\rm{St}}}}(m, p, {\mathbb C})\times {\bf{T}}_{V}{{{\rm{St}}}}(n, p, {\mathbb C})\\ &=&\big\{(\xi, \eta)\in {\mathbb C}^{m\times p}\times {\mathbb C}^{n\times p}|\xi^HU+U^H\xi=\eta^HV+V^H\eta=0\big\}. \end{eqnarray*} $

因复乘积流形$ {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $是欧式空间$ {\mathbb C}^{m\times p}\times {\mathbb C}^{n\times p} $上的嵌入子流形, 故其上的黎曼度量可定义为

(2.3) $ \begin{equation} \langle (\xi_1, \eta_1), (\xi_2, \eta_2)\rangle_{(U, V)}={{{\rm{Tr}}}}(\xi_1^H\xi_2)+{{{\rm{Tr}}}}(\eta_1^H\eta_2), \end{equation} $

其中$ (U, V)\in {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $, $ (\xi_1, \eta_1), (\xi_2, \eta_2)\in {\bf{T}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $, 其诱导范数为$ \|(\xi_1, \eta_1)\|_{(U, V)}={{{\rm{Tr}}}}(\xi_1^H\xi_1)+{{{\rm{Tr}}}}(\eta_1^H\eta_1) $. 为讨论方便, 下文用$ \langle\cdot, \cdot\rangle $和$ \|\cdot\| $分别表示复乘积流形$ {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $上的黎曼度量及其诱导范数. 在黎曼度量(2.3)下, $ \forall (B, C)\in {\mathbb C}^{m\times p}\times{\mathbb C}^{n\times p} $在$ (U, V) $对应的切空间$ {\bf{T}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $上的正交投影可表示为

(2.4) $ \begin{equation} {\bf P}_{( {U}, {V})}( {B}, {C})=\left({\bf P}_{ {U}}( {B}), \ {\bf P}_{ {V}}( {C}) \right), \end{equation} $

其中

(2.5) $ \begin{equation} \begin{array}{cc} {\bf P}_{ {U}}( {B})={B}- {U}\ {\rm her}( {U}^H {B}), & {\bf P}_{ {V}}( {C})={C}- {V}\ {\rm her}( {V}^H {C}), \end{array} \end{equation} $

分别表示$ B\in {\mathbb C}^{m\times p} $在切空间$ { {\bf{T}}}_{ {U}} {\rm St}(m, p, {\mathbb C}) $上的正交投影, 和$ C\in {\mathbb C}^{n\times p} $在切空间$ { {\bf{T}}}_{ {V}} {\rm St}(n, $$ p, {\mathbb C}) $上的正交投影. 另$ \forall W\in {\mathbb C}^{p\times p} $, 记$ {{{\rm{her}}}}(W):=\frac{1}{2}(W+W^H) $.

下面讨论问题(1.2)中目标函数$ {F}(U, V) $所对应的黎曼梯度和黎曼海塞. 记$ \bar{F} $为定义在$ {\mathbb C}^{m\times p}\times {\mathbb C}^{n\times p} $上与$ {F}(U, V) $具有相同形式的函数, 即

(2.6) $ \begin{equation} \bar{F}(U, V):=-{{{\rm{Tr}}}}(U^HAV\Theta), (U, V)\in {\mathbb C}^{m\times p}\times {\mathbb C}^{n\times p}. \end{equation} $

引理2.1   目标函数$ {F}(U, V) $在$ (U, V)\in {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $的黎曼梯度$ {{{\rm{grad }}}} F $可表示为

(2.7) $ \begin{equation} {{{\rm{grad }}}} F(U, V)=\left({\bf P}_U(-AV\Theta), {\bf P}_V(-A^HU\Theta)\right). \end{equation} $

证   因为复乘积流形$ {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $是欧式空间$ {\mathbb C}^{m\times p}\times {\mathbb C}^{n\times p} $上的黎曼子流形, 故目标函数$ {F}(U, V) $的黎曼梯度$ {{{\rm{grad }}}} F $等于函数$ \bar{F}(U, V) $的欧式梯度$ \nabla \bar{F} $在切空间$ {\bf{T}}_{(U, V)} $$ {{{\rm{St}}}}(m, p, {\mathbb C}) $$ \times {{{\rm{St}}}}(n, p, {\mathbb C}) $上的正交投影^[13]. 而通过简单计算可知, 函数$ \bar{F}(U, V) $在$ (U, V)\in{\mathbb C}^{m\times p}\times {\mathbb C}^{n\times p} $的欧式梯度$ \nabla \bar{F} $可表示为

(2.8) $ \begin{equation} \nabla\overline{F}(U, V)=\left(\nabla_U\overline{F}(U, V), \ \nabla_V\overline{F}(U, V)\right), \end{equation} $

其中

(2.9) $ \begin{equation} \nabla_U \bar{F}(U, V)=-AV\Theta, \nabla_V \bar{F}(U, V)=-A^HU\Theta. \end{equation} $

进而由(2.4)式和(2.5)式可得

(2.10) $ \begin{equation} {{{\rm{grad }}}} F(U, V)={ {\bf{P}}}_{(U, V)}\left(\nabla\overline{F}(U, V)\right)=\left({ {\bf{P}}}_U(-AV\Theta), { {\bf{P}}}_V(-A^HU\Theta)\right). \end{equation} $

引理2.1得证.

引理2.2   令$ (\xi, \eta)\in {\bf{T}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $, 目标函数$ F(U, V) $在$ (\xi, \eta) $处的黎曼海塞是切空间$ {\bf{T}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $上的线性变换, 其具体形式为

(2.11) $ \begin{equation} {{{\rm{Hess }}}} F(U, V)[(\xi, \eta)]=\left( {\bf{P}}_{U}\left(\xi{{{\rm{her}}}}(U^HAV\Theta)-A\eta\Theta\right), {\bf{P}}_{V}\left(\eta{{{\rm{her}}}}(V^HA^HU\Theta)-A^H\xi\Theta\right)\right). \end{equation} $

证   $ \forall (\xi, \eta)\in {\mathbb C}^{m\times p}\times {\mathbb C}^{n\times p} $, 经简单计算, 函数$ \bar{F}(U, V) $在$ (\xi, \eta) $的欧式海塞$ \nabla^2 \bar{F} $可表示

(2.12) $ \begin{equation} \nabla^2 \bar{F}(U, V)[(\xi, \eta)]=\big(-A\eta\Theta, -A^H\xi \Theta\big). \end{equation} $

利用$ \nabla^2 \overline{F} $和正交投影(2.4)–(2.5), $ F(U, V) $在切空间$ {\bf{T}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $的黎曼海塞可表示为^[15]

(2.13) $ \begin{equation} {{{\rm{Hess }}}} F(U, V)[(\xi, \eta)]={\bf P}_{(U, V)}\big({\mathcal D}({{{\rm{grad }}}} F)(U, V)[(\xi, \eta)]\big). \end{equation} $

注意到$ {{{\rm{grad }}}} {F}(U, V)={\bf P}_{(U, V)}\left(\nabla \bar{F}(U, V)\right) $. 我们把这个关系式的右边看成$ U $和$ V $的两个矩阵函数的乘积, 而不是一个映射和一个函数的复合. 参考文献[15], (2.13) 式可进一步表示为

(2.14) $ \begin{eqnarray} {{{\rm{Hess }}}} F(U, V)[(\xi, \eta)]&=&{\bf P}_{(U, V)}\Big( {\mathcal D}({{{\rm{grad }}}} F)(U, V)[(\xi, \eta)]\Big){}\\ & =&{\bf P}_{(U, V)}\Big( {\mathcal D}({\bf P}_{(U, V)}(\nabla \bar{F}))(U, V)[(\xi, \eta)]\Big){}\\ & =&{\bf P}_{(U, V)}\Big( {\mathcal D}{\bf P}_{(U, V)}[(\xi, \eta)](\nabla \bar{F}(U, V)) +{\bf P}_{(U, V)}\left({\mathcal D}(\nabla \bar{F})(U, V)[(\xi, \eta)]\right)\Big){}\\ & =&{\bf P}_{(U, V)}\Big( {\mathcal D}{\bf P}_{(U, V)}[(\xi, \eta)](\nabla \bar{F}(U, V)) +\nabla^2\bar{F}(U, V)[(\xi, \eta)]\Big), \end{eqnarray} $

在(2.14)式的最后一个等式中用到了$ {\bf P}^2_{(U, V)}={\bf P}_{(U, V)} $. 又因为

(2.15) $ \begin{eqnarray} {\mathcal D}{\bf P}_{(U, V)}[(\xi, \eta)](\nabla \bar{F}(U, V)) &=&D{\bf P}_{(U, V)}[(\xi, \eta)](\nabla_U \bar{F}(U, V), \nabla_V \bar{F}(U, V)){}\\ &=&{\mathcal D}\Big(\nabla_U \bar{F}(U, V)-U{\rm her}(U^H\nabla_U \bar{F}(U, V)), {}\\ & &\nabla_V \bar{F}(U, V)-V{\rm her} (V^H\nabla_V \bar{F}(U, V))\Big)[(\xi, \eta)]{}\\ &=&\Big(-U\ {\rm her}(\xi^H\nabla_U \bar{F}(U, V))-\xi\ {\rm her}(U^H\nabla_U \bar{F}(U, V)), {}\\ &&-V\ {\rm her}(\eta^H\nabla_V \bar{F}(U, V))-\eta\ {\rm her}(V^H\nabla_V \bar{F}(U, V))\Big). \end{eqnarray} $

记$ S_1 $和$ S_2 $均为$ p\times p $阶Herimite矩阵, 则有如下等式成立

(2.16) $ \begin{equation} {\bf P}_{(U, V)}(US_1, VS_2)=\big(US_1-U{{{\rm{her}}}}(U^HUS_1), VS_2-V{{{\rm{her}}}}(V^HVS_2)\big)=O. \end{equation} $

利用(2.16)式的结论可以简化(2.14)式最后一个等式中$ {\bf P}_{(U, V)}\left({\mathcal D}{\bf P}_{(U, V)}[(\xi, \eta)](\nabla \bar{F}(U, V))\right) $的计算. 事实上, 结合(2.15)式和(2.16)式可得

(2.17) $ \begin{equation} {\bf{P}}_{(U, V)}\Big({\mathcal D} {\bf{P}}_{(U, V)}[(\xi, \eta)](\nabla \overline{f}(U, V))\Big)= {\bf{P}}_{(U, V)}\Big(\xi {{{\rm{her}}}}(U^HAV\Theta), \ \eta {{{\rm{her}}}}(V^HA^HU\Theta)\Big). \end{equation} $

根据(2.14)式的最后一个等式, $ F(U, V) $在$ (\xi, \eta) $处的黎曼海塞$ {{{\rm{Hess }}}} F $可表示为

$ \begin{eqnarray*} \begin{array}{rl} {{{\rm{Hess }}}} F(U, V)[(\xi, \eta)]= {\bf{P}}_{(U, V)}\Big( {\mathcal D} {\bf{P}}_{(U, V)}[(\xi, \eta)](\nabla \bar{F}(U, V)) +\nabla^2\bar{F}(U, V)[(\xi, \eta)]\Big). \end{array} \end{eqnarray*} $

将(2.8)式和(2.17)代入上式即得(2.11). 引理2.2得证.

本节最后给出下文需用到的若干矩阵拉直算子的定义及相关性质. $ \forall B=(b_{ij})\in {\mathbb C}^{m\times n} $, 定义拉直算子

$ {{{\rm{vec}}}}(B)=(b_{11}, \cdots , b_{m1}, b_{12}, \cdots , b_{m2}, b_{1n}, \cdots , b_{mn})^T. $

易知$ \forall B, B'\in {\mathbb C}^{n\times n} $有$ {{{\rm{Tr}}}}(B^HB')={{{\rm{vec}}}}(B)^H{{{\rm{vec}}}}(B') $. 关于拉直算子$ {{{\rm{vec}}}} $和Kronecker积$ \otimes $, 有如下性质

1. $ \forall B\in {\mathbb C}^{m\times n} $, $ X\in {\mathbb C}^{n\times s} $和$ C\in {\mathbb C}^{s\times t} $, 有$ {{{\rm{vec}}}}(BXC)=(C^T\otimes B){{{\rm{vec}}}}(X) $.

2. $ \forall B\in {{\mathbb R}}^{m\times n} $, 存在$ mn\times mn $阶置换矩阵$ T_{(m, n)}=\sum\limits_{i=1}^n\sum\limits_{j=1}^m E_{ij}\otimes E_{ij}^T $ (参见文献[16]), 使得

(2.18) $ \begin{equation} {{{\rm{vec}}}}(B)=T_{(m, n)}{{{\rm{vec}}}}(B^T), \end{equation} $

其中$ E_{ij} $是指标位置$ (i, j) $处为$ 1 $其他位置为零的$ n\times m $阶矩阵. $ T_{(m, n)} $是置换矩阵故正交, 结合$ T_{(m, n)}T_{(n, m)}=I_{mn} $ (参见文献[16]), 可得$ T_{(m, n)}^{-1}=T_{(m, n)}^{T}=T_{(n, m)}, $且

(2.19) $ \begin{equation} {{{\rm{vec}}}}(B^T)=T_{(n, m)}{{{\rm{vec}}}}(B). \end{equation} $

3. $ \forall B\in {{\mathbb R}}^{n\times n} $, 记$ {{{\rm{sym}}}}(B)=\frac{1}{2}(B+B^T) $和$ {{{\rm{skew}}}}(B)=\frac{1}{2}(B-B^T) $, 有

(2.20) $ \begin{equation} {{{\rm{vec}}}}({{{\rm{sym}}}}(B))=\frac{1}{2}(I_{n^2}+T_{(n, n)}){{{\rm{vec}}}}(B), \ \ \ {{{\rm{vec}}}}({{{\rm{skew}}}}(B))=\frac{1}{2}(I_{n^2}-T_{(n, n)}){{{\rm{vec}}}}(B). \end{equation} $

记$ {\mathbb S}^{n\times n} $和$ {\mathbb AS}^{n\times n} $分别表示所有$ n\times n $阶对称矩阵和反对称矩阵构成的集合. 下面给出对称矩阵和反对称矩阵拉直算子的定义.

定义2.1   $ \forall B=(b_{ij})\in {\mathbb S}^{n\times n} $, 定义对称矩阵拉直算子$ {{{\rm{vec}}}}_{S}(B)\in {{\mathbb R}}^{{n(n+1)}/{2}} $如下

(2.21) $ \begin{equation} {{{\rm{vec}}}}_{S}(B)=(b_{11}, \sqrt{2}b_{21}, \cdots , \sqrt{2}b_{n1}, b_{22}, \sqrt{2}b_{32}, \cdots , \sqrt{2}b_{n2}, \cdots , b_{(n-1)(n-1)}, \sqrt{2}b_{n(n-1)}, b_{nn} )^T. \end{equation} $

$ \forall C=(c_{ij})\in {\mathbb AS}^{n\times n} $, 定义反对称矩阵拉直算子$ {{{\rm{vec}}}}_{K}(C)\in {{\mathbb R}}^{{n(n-1)}/{2}} $如下

(2.22) $ \begin{equation} {{{\rm{vec}}}}_{K}(C)=\sqrt{2}(c_{21}, c_{31}, \cdots , c_{n1}, c_{32}, c_{42}, \cdots , c_{n2}, \cdots , c_{(n-1)(n-2)}, c_{n(n-2)}, c_{n(n-1)} )^T. \end{equation} $

由定义易得如下等式

(2.23) $ \begin{eqnarray} && \text{tr}(BB')={{{\rm{vec}}}}_{S}(B)^T{{{\rm{vec}}}}_{S}(B'), \forall B, B'\in {\mathbb S}^{n\times n}, \end{eqnarray} $

(2.24) $ \begin{eqnarray} && \text{tr}(C^TC')={{{\rm{vec}}}}_{K}(C)^T{{{\rm{vec}}}}_{K}(C'), \forall C, C'\in {\mathbb AS}^{n\times n}. \end{eqnarray} $

引理2.3^[17]   $ \forall B\in {\mathbb S}^{n\times n} $和$ \forall C\in {\mathbb AS}^{n\times n} $, 存在$ S_n\in {{\mathbb R}}^{n^2\times {n(n+1)}/{2}} $和$ K_n\in {{\mathbb R}}^{n^2\times {n(n-1)}/{2}} $使得

$ {{{\rm{vec}}}}(B)=S_n{{{\rm{vec}}}}_S(B), {{{\rm{vec}}}}(C)=K_n{{{\rm{vec}}}}_K(C). $

矩阵$ S_n $和$ K_n $的具体表达式可见文献[17], 且$ S_n $和$ K_n $均为列正交矩阵, 即满足$ S_n^TS_n=I_{{n(n+1)}/{2}} $和$ K_n^TK_n=I_{{n(n-1)}/{2}} $.

为建立复数和实向量, 复矩阵和实矩阵对之间联系, 定义同构映射$ \cong $, 即

$ c_1+{\rm i} c_2=q \cong \ \tilde{q}=[c_1 c_2]. $

$ \forall B={{{\rm{Re}}}}(B)+{\rm i} {{{\rm{Im}}}}(B)\in {\mathbb C}^{m\times n} $, 有$ B\cong \widetilde{B}=[{{{\rm{Re}}}}(B) {{{\rm{Im}}}}(B)] $, 且

$ {{{\rm{vec}}}}({{{\rm{Re}}}}(B))+{\rm i} {{{\rm{vec}}}}({{{\rm{Im}}}}(B))={{{\rm{vec}}}}(B)\cong {{{\rm{vec}}}}(\widetilde{B})=\left[\begin{array}{ccccc}{{{\rm{vec}}}}({{{\rm{Re}}}}(B))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(B))\end{array}\right]. $

$ \forall B={{{\rm{Re}}}}(B)+{\rm i} {{{\rm{Im}}}}(B) $和$ C={{{\rm{Re}}}}(C)+{\rm i} {{{\rm{Im}}}}(C) $, 可得

$ ({{{\rm{Re}}}}(B)+{{{\rm{Re}}}}(C))+{\rm i} ({{{\rm{Im}}}}(B)+{{{\rm{Im}}}}(C))=(B+C)\cong \widetilde{B+C}=[{{{\rm{Re}}}}(B)+{{{\rm{Re}}}}(C) {{{\rm{Im}}}}(B)+{{{\rm{Im}}}}(C)]. $

下面给出拉直算子$ {{{\rm{vec}}}} $, $ {{{\rm{vec}}}}_S $, $ {{{\rm{vec}}}}_K $以及同构映射$ \cong $的若干结论, 其证明较简单, 故略.

1. $ \forall B\in {\mathbb C}^{m\times n} $, 可得

(2.25) $ \begin{equation} \begin{array}{l} {{{\rm{vec}}}}({B^H})\cong {{{\rm{vec}}}}({\widetilde{B^H}})=\left[\begin{array}{ccccc}T_{(n, m)}&O\\ O& -T_{(n, m)}\end{array}\right]\left[\begin{array}{ccccc}{{{\rm{vec}}}}({{{\rm{Re}}}}(B))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(B))\end{array}\right]. \end{array} \end{equation} $

$ \forall B\in {\mathbb C}^{n\times n} $, 记$ {{{\rm{skewH}}}}(B):=\frac{1}{2}(B-B^H) $, 可得

(2.26) $ \begin{equation} {{{\rm{vec}}}}({{{\rm{skewH}}}}(B))\cong {{{\rm{vec}}}}(\widetilde{{{{\rm{skewH}}}}(B)})=\left[\begin{array}{ccccc} { } \frac{1}{2}(I_{n^2}-T_{(n, n)})&O\\ O& { } \frac{1}{2}(I_{n^2}+T_{(n, n)})\end{array}\right]\left[\begin{array}{ccccc}{{{\rm{vec}}}}({{{\rm{Re}}}}(B))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(B))\end{array}\right]. \end{equation} $

2. 记$ {\mathbb AH}^{n\times n} $表示所有$ n\times n $阶反Hermite矩阵构成的集合. $ \forall B\in {\mathbb AH}^{n\times n} $, 可得

(2.27) $ \begin{equation} {{{\rm{vec}}}}(B)\cong {{{\rm{vec}}}}(\widetilde{B})\left[\begin{array}{ccccc}K_n\ &O\\ O\ & S_n\end{array}\right] \left[\begin{array}{ccccc}{{{\rm{vec}}}}_K({{{\rm{Re}}}}(B))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(B))\end{array}\right]. \end{equation} $

3. $ \forall B\in {\mathbb C}^{m\times n} $, $ X\in {\mathbb C}^{n\times s} $和$ C\in {\mathbb C}^{s\times t} $, 可得

(2.28) $ \begin{equation} {{{\rm{vec}}}}(BXC)\cong{{{\rm{vec}}}}(\widetilde{BXC})=\left[\begin{array}{ccccc}T_1\\ T_2\end{array}\right]\left[\begin{array}{ccccc}{{{\rm{vec}}}}({{{\rm{Re}}}}(X))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(X))\end{array}\right], \end{equation} $

其中$ T=\Big[ C^T\otimes B \ i(C^T\otimes B)\Big] $, $ T_1={{{\rm{Re}}}}(T) $, $ T_2={{{\rm{Im}}}}(T). $

4. $ \forall B\in {\mathbb C}^{m\times n} $, $ X\in {\mathbb AH}^{n\times n} $和$ C\in {\mathbb C}^{n\times t} $, 可得

(2.29) $ \begin{equation} {{{\rm{vec}}}}(BXC)\cong{{{\rm{vec}}}}(\widetilde{BXC})=\left[\begin{array}{ccccc}W_1\\ W_2\end{array}\right]\left[\begin{array}{ccccc}{{{\rm{vec}}}}_K({{{\rm{Re}}}}(X))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(X))\end{array}\right], \end{equation} $

其中$ W=\Big[ (C^T\otimes B)K_n \ i(C^T\otimes B)S_n\Big] $, $ W_1={{{\rm{Re}}}}(W) $, $ W_2={{{\rm{Im}}}}(W). $

对于迭代框架(2.1), 通过求解如下黎曼牛顿方程得到黎曼牛顿法在当前迭代步$ (U^{(k)}, $$ V^{(k)}) \in {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $的搜索方向$ (\xi^{(k)}, \eta^{(k)}) \in { {\bf{T}}}_{(U^{(k)}, V^{(k)})}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $

(3.1) $ \begin{equation} {{{\rm{Hess }}}} F(U^{(k)}, V^{(k)})[(\xi, \eta)]=-{{{\rm{grad }}}} F(U^{(k)}, V^{(k)}). \end{equation} $

若已求得$ (\xi^{(k)}, \eta^{(k)}) $, 则更新迭代步由下式给出

(3.2) $ \begin{equation} \begin{array}{l} (U^{(k+1)}, V^{(k+1)})={\mathcal R}_{(U^{(k)}, V^{(k)})}(\xi^{(k)}, \eta^{(k)}), \end{array} \end{equation} $

其中牛顿法的迭代步长通常设置为1. 将(2.7), (2.11)式代入(3.1) 式, 得到黎曼牛顿方程的具体形式为(为方便讨论省略上标)

(3.3) $ \begin{equation} \left\{ \begin{array}{l} {\bf{P}}_{U}\left(\xi {{{\rm{her}}}}(U^HAV\Theta)-A\eta\Theta\right) = {\bf{P}}_{U}\left(AV\Theta \right), \\ {\bf{P}}_{V}\left( \eta{{{\rm{her}}}}(V^HA^HU\Theta)-A^H\xi\Theta\right)= {\bf{P}}_{V}\left( A^HU\Theta\right), \end{array}\right. \end{equation} $

其中搜索方向$ (\xi, \eta)\in {{\bf T}}_{( {U}, {V})}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $.

下面讨论如何求解牛顿方程(3.3). 注意到搜索方向$ (\xi, \eta)\in { {\bf{T}}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $, 即变量$ (\xi, \eta) $须同时满足$ \xi^HU+U^H\xi=O $和$ \eta^HV+V^H\eta=O $, 故该牛顿方程求解比较复杂. 注意到黎曼海塞是切空间$ { {\bf{T}}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $上的线性变换, 考虑从该线性变换的变换矩阵入手, 将牛顿方程(3.3) 转化为容易求解的线性方程组. 为此需引入复乘积流形切空间的另一种表达方式.

引理3.1^[12]   复乘积流形$ {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $在$ (U, V)\in {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $处切空间$ {\bf{T}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $的等价表达方式为

$ \begin{eqnarray*} && {\bf{T}}_{(U, V)}\big({{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C})\big)\\ &=&\Big\{(UE+U_{\bot}F, VM+V_{\bot}N)\ \big|\ E, M\in {\mathbb AH}^{p\times p}, F\in {\mathbb C}^{(m-p)\times p}, N\in {\mathbb C}^{(n-p)\times p}\Big\}. \end{eqnarray*} $

其中$ U_{\bot}\in {\mathbb C}^{m\times (m-p)} $, $ V_{\bot}\in {\mathbb C}^{n\times (n-p)} $分别表示矩阵$ U $和$ V $的正交补, 即满足$ U^HU_{\bot}=O $, $ UU^H+U_{\bot}U_{\bot}^H=I_{m} $和$ V^HV_{\bot}=O $, $ VV^H+V_{\bot}V_{\bot}^H=I_{n} $.

由引理3.1可知

$ \begin{eqnarray*} \mbox{dim}\left( {\bf{T}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C})\right) &=&\mbox{dim}(E)+\mbox{dim}(F)+\mbox{dim}(M)+\mbox{dim}(N)\\ &=&2p(m+n-p). \end{eqnarray*} $

令$ (\xi, \eta)\in {\bf{T}}_{(U, V)}\big({{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C})\big) $是在点$ (U, V) $处的搜索方向. 由引理3.1, $ (\xi, \eta) $可表示为

(3.4) $ \begin{equation} (\xi, \eta)=(UE+U_{\bot}F, VM+V_{\bot}N), \end{equation} $

其中$ E, M\in {\mathbb AH}^{p\times p}, F\in {\mathbb C}^{(m-p)\times p}, N\in {\mathbb C}^{(n-p)\times p} $. 因$ {{{\rm{Hess }}}} F(U, V)[(\xi, \eta)] $也在切空间上, 故也存在唯一的矩阵$ E_H, M_H\in {\mathbb AH}^{p\times p}, F_H\in {\mathbb C}^{(m-p)\times p}, N_H\in {\mathbb C}^{(n-p)\times p} $, 使得

(3.5) $ \begin{equation} {{{\rm{Hess }}}} F(U, V)[(\xi, \eta)]=(UE_H+U_{\bot}F_H, VM_H+V_{\bot}N_H). \end{equation} $

下面引理给出了由$ E, M, F, N $表示$ E_H, M_H, F_H, N_H $的具体表示方式.

引理3.2   $ \forall (U, V)\in {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $, 搜索方向$ (\xi, \eta)\in {\bf{T}}_{(U, V)}\big({{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C})\big) $表示为(3.4)式, 则对应黎曼海塞可表示为$ {{{\rm{Hess }}}} F(U, V)[(\xi, \eta)]=(UE_H+U_{\bot}F_H, $$ VM_H+V_{\bot}N_H), $其中

(3.6) $ \begin{eqnarray} && E_H={{{\rm{skewH}}}}\left(E{{{\rm{her}}}}(U^HAV\Theta)-U^HAVM\Theta-U^HAV_{\bot}N\Theta\right), \end{eqnarray} $

(3.7) $ \begin{eqnarray} && F_H=F{{{\rm{her}}}}(U^HAV\Theta)-U_{\bot}^HAVM\Theta-U_{\bot}^HAV_{\bot}N\Theta, \end{eqnarray} $

(3.8) $ \begin{eqnarray} && M_H={{{\rm{skewH}}}}\left(M{{{\rm{her}}}}(V^HA^HU\Theta)-V^HA^HUE\Theta-V^HA^HU_{\bot}F\Theta\right), \end{eqnarray} $

(3.9) $ \begin{eqnarray} && N_H=N{{{\rm{her}}}}(V^HA^HU\Theta)-V_{\bot}^HA^HUE\Theta-V_{\bot}^HA^HU_{\bot}F\Theta. \end{eqnarray} $

证   结合引理2.2中黎曼海塞的表达式(2.11)和(3.5), 可知

(3.10) $ \begin{eqnarray} &&UE_H+U_{\bot}F_H= {\bf{P}}_{U}\left(\xi{{{\rm{her}}}}(U^HAV\Theta)-A\eta\Theta\right), \end{eqnarray} $

(3.11) $ \begin{eqnarray} &&VM_H+V_{\bot}N_H= {\bf{P}}_{V}\left(\eta{{{\rm{her}}}}(V^HA^HU\Theta)-A^H\xi\Theta\right). \end{eqnarray} $

注意到$ \forall W_1 \in {\mathbb C}^{m \times p} $有$ U^{H} {\bf{P}}_{U}(W_1)={{{{\rm{skewH}}}}}\left(U^{H} W_1\right) $和$ U_{\perp}^{H} {\bf{P}}_{U}(W_1)=U_{\perp}^{H} W_1 $. (3.10)式两边同乘$ U^H $, 结合$ U^{H} U=I_{p} $和$ U^{H} U_{\perp}=O $可得

(3.12) $ \begin{equation} E_{H}=U^{H} {\bf{P}}_{U}\left(\xi{{{\rm{her}}}}(U^HAV\Theta)-A\eta\Theta\right)={{{{\rm{skewH}}}}}\left(U^{H}\xi{{{\rm{her}}}}(U^HAV\Theta)-U^{H}A\eta\Theta\right). \end{equation} $

进一步由$ U^H\xi=U^H(UE+U_{\bot}F)=E $, 将$ \eta=VM+V_{\bot}N $代入到(3.12)式可得(3.6)式. 等式(3.10)两边同乘$ U_{\bot}^H $, 可得

$ \begin{eqnarray*} F_{H}=U_\bot^{H} {\bf{P}}_{U}\left(\xi{{{\rm{her}}}}(U^HAV\Theta)-A\eta\Theta\right)=U_\bot^{H}\xi{{{\rm{her}}}}(U^HAV\Theta)-U_\bot^{H}A\eta\Theta. \end{eqnarray*} $

由$ U_{\perp}^{H} \xi=U_{\perp}^{H} (U E+U_{\perp} F)=F $, 上式代入$ \eta=VM+V_{\bot}N $可得(3.7)式. 同理, (3.11)式两边分别同乘$ V^H $和$ V_{\bot}^H $, 并结合$ \forall W_2 \in {\mathbb C}^{n \times p} $有$ V^{H} {\bf{P}}_{V}(W_2)={{{{\rm{skewH}}}}}\left(V^{H} W_2\right) $和$ V_{\perp}^{H} {\bf{P}}_{V}(W_2)=V_{\perp}^{H} W_2 $可得

$ \begin{eqnarray*} &&M_H={{{\rm{skewH}}}}\left(V^H\eta{{{\rm{her}}}}(V^HA^HU\Theta)-V^HA^H\xi\Theta \right), \\ &&N_H=V_\bot^H\eta{{{\rm{her}}}}(V^HA^HU\Theta)-V_\bot^HA^H\xi\Theta. \end{eqnarray*} $

进而由$ V^H\eta=M $和$ V_{\bot}^H\eta=N $, 并代入$ \xi=UE+U_{\bot}F $, 分别可得(3.8)式和(3.9)式. 引理3.2得证.

记$ \kappa=2p(m+n-p) $. 因黎曼海塞是切空间$ {\bf{T}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $上的线性变换, 注意到切空间$ {\bf{T}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $和向量空间$ {{\mathbb R}}^{\kappa} $在黎曼海塞这个线性映射下是同构的, 故黎曼海塞也可视为$ {{\mathbb R}}^{\kappa} $上的线性变换, 即把一个$ \kappa $维向量变换为另一个$ \kappa $维向量. 下面引理给出该线性变换的变换矩阵.

引理3.3   令$ {\mathcal H} $为$ {{\mathbb R}} ^{\kappa} $上作用于$ \big({{{\rm{vec}}}}_K({{{\rm{Re}}}}(E))^T, $$ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(E))^T, $$ {{{\rm{vec}}}}({{{\rm{Re}}}}(F))^T, $$ {{{\rm{vec}}}}({{{\rm{Im}}}}(F))^T, $$ {{{\rm{vec}}}}_K({{{\rm{Re}}}}(M))^T, $$ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(M))^T, $$ {{{\rm{vec}}}}({{{\rm{Re}}}}(N))^T, $$ {{{\rm{vec}}}}({{{\rm{Im}}}}(N))^T\big)^T $的线性变换, $ E, $$ M\in {\mathbb AH}^{p\times p}, $$ F\in {\mathbb C}^{(m-p)\times p}, N\in {\mathbb C}^{(n-p)\times p} $, 且

(3.13) $ \begin{equation} {\mathcal H}\left[\begin{array}{ccccc} {{{\rm{vec}}}}_K({{{\rm{Re}}}}(E))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(E))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(F))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(F))\\ {{{\rm{vec}}}}_K({{{\rm{Re}}}}(M))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(M))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(N))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(N)) \end{array}\right]= \left[\begin{array}{ccccc}{{{\rm{vec}}}}_K({{{\rm{Re}}}}(E_H))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(E_H))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(F_H))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(F_H))\\ {{{\rm{vec}}}}_K({{{\rm{Re}}}}(M_H))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(M_H))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(N_H))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(N_H)) \end{array}\right], \end{equation} $

其中$ E_H $, $ F_H $, $ M_H $, $ N_H $在(3.6)–(3.9)式给出, 线性变换$ {\mathcal H} $的变换矩阵$ H_{A} $由下式给出

(3.14) $ \begin{equation} H_A=\left[\begin{array}{cccccccc} Q_{E, 1}&\ Q_{F, 1}\ &Q_{M, 1}&\ Q_{N, 1}\\ Q_{E, 2}&Q_{F, 2}&Q_{M, 2}&\ Q_{N, 2}\\ R_{E, 1}&R_{F, 1}&R_{M, 1}&\ R_{N, 1}\\ R_{E, 2}&R_{F, 2}&R_{M, 2}&\ R_{N, 2}\\ S_{E, 1}&S_{F, 1}&S_{M, 1}&\ S_{N, 1}\\ S_{E, 2}&S_{F, 2}&S_{M, 2}&\ S_{N, 2}\\ T_{E, 1}&T_{F, 1}&T_{M, 1}&\ T_{N, 1}\\ T_{E, 2}&T_{F, 2}&T_{M, 2}&\ T_{N, 2} \end{array}\right]\in {{\mathbb R}}^{\kappa\times \kappa}, \end{equation} $

其中

(3.15) $ \begin{equation} \begin{array}{ll} { } Q_{E, 1}=\frac{1}{2} K_{p}^{T}\left(I_{p^{2}}-T_{(p, p)}\right) {\rm Re}\left(Q_{E}\right), \quad Q_{E, 2}=\frac{1}{2} S_{p}^{T}\left(I_{p^{2}}+T_{(p, p)}\right) {\rm Im}\left(Q_{E}\right). \\ Q_{F, 1}={O}, \quad Q_{F, 2}={O}. \\ { } Q_{M, 1}=-\frac{1}{2} K_{p}^{T}\left(I_{p^{2}}-T_{(p, p)}\right) {\rm Re}\left(Q_{M}\right), \quad Q_{M, 2}=-\frac{1}{2} S_{p}^{T}\left(I_{p^{2}}+T_{(p, p)}\right) {\rm Im}\left(Q_{M}\right). \\ { } Q_{N, 1}=-\frac{1}{2} K_{p}^{T}\left(I_{p^{2}}-T_{(p, p)}\right) {\rm Re}\left(Q_{N}\right), \quad Q_{N, 2}=-\frac{1}{2} S_{p}^{T}\left(I_{p^{2}}+T_{(p, p)}\right) {\rm Im}\left(Q_{N}\right). \end{array} \end{equation} $

(3.16) $ \begin{equation} \begin{array}{ll} R_{E, 1} ={O}, R_{E, 2}={O}, R_{F, 1} ={\rm Re}\left(R_{F}\right), R_{F, 2}={\rm Im}\left(R_{F}\right). \\ R_{M, 1} =-{\rm Re}\left(R_{M}\right), R_{M, 2}=-{\rm Im}\left(R_{M}\right), R_{N, 1} =-{\rm Re}\left(R_{N}\right), R_{N, 2} =-{\rm Im}\left(R_{N}\right). \end{array} \end{equation} $

(3.17) $ \begin{equation} \begin{array}{ll} { } S_{E, 1}=-\frac{1}{2} K_{p}^{T}\left(I_{p^{2}}-T_{(p, p)}\right) {\rm Re}\left(S_{E}\right), \quad S_{E, 2}=-\frac{1}{2} S_{p}^{T}\left(I_{p^{2}}+T_{(p, p)}\right) {\rm Im}\left(S_{E}\right). \\ { } S_{F, 1}=-\frac{1}{2} K_{p}^{T}\left(I_{p^{2}}-T_{(p, p)}\right) {\rm Re}\left(S_{F}\right), \quad S_{F, 2}=-\frac{1}{2} S_{p}^{T}\left(I_{p^{2}}+T_{(p, p)}\right) {\rm Im}\left(S_{F}\right). \\ [3mm] { } S_{M, 1}=\frac{1}{2} K_{p}^{T}\left(I_{p^{2}}-T_{(p, p)}\right) {\rm Re}\left(S_{M}\right), \quad S_{M, 2}=\frac{1}{2} S_{p}^{T}\left(I_{p^{2}}+T_{(p, p)}\right) {\rm Im}\left(S_{M}\right). \\ S_{N, 1}={O}, \quad S_{N, 2}={O}. \end{array} \end{equation} $

(3.18) $ \begin{equation} \begin{array}{llll} T_{E, 1}=-{\rm Re}\left(T_{E}\right), T_{E, 2}=-{\rm Im}\left(T_{E}\right), T_{F, 1} =-{\rm Re}\left(T_{F}\right), T_{F, 2}=-{\rm Im}\left(T_{F}\right), \\ T_{M, 1} =O, T_{M, 2}=O, T_{N, 1} ={\rm Re}\left(T_{N}\right), T_{N, 2} ={\rm Im}\left(T_{N}\right), \end{array} \end{equation} $

其中

(3.19) $ \begin{equation} \begin{array}{ll} Q_{E} =\Big[\big(({{{\rm{her}}}}(U^HAV\Theta))^T\otimes I_p\big)K_{p} i\big(({{{\rm{her}}}}(U^HAV\Theta))^T\otimes I_p\big)S_{p}\Big], \\ Q_{M} =\Big[(\Theta\otimes(U^HAV))K_{p} i(\Theta\otimes(U^HAV))S_{p}\Big], \\ Q_{N} =\Big[\Theta\otimes(U^HAV_{\bot}) i(\Theta\otimes(U^HAV_{\bot}))\Big]; \end{array} \end{equation} $

(3.20) $ \begin{equation} \begin{array}{ll} R_{F}=\Big[({{{\rm{her}}}}(U^HAV\Theta))^T\otimes I_{m-p} i\big(({{{\rm{her}}}}(U^HAV\Theta))^T\otimes I_{m-p}\big)\Big], \\ R_{M} =\Big[(\Theta\otimes(U_{\bot}^HAV))K_{p} i(\Theta\otimes(U_{\bot}^HAV))S_{p}\Big], \\ R_{N}=\Big[\Theta\otimes(U_{\bot}^HAV_{\bot}) i(\Theta\otimes(U_{\bot}^HAV_{\bot}))\Big]; \end{array} \end{equation} $

(3.21) $ \begin{equation} \begin{array}{ll} S_{E} =\Big[(\Theta\otimes(V^HA^HU))K_{p} i(\Theta\otimes(V^HA^HU))S_{p}\Big], \\ S_{F}=\Big[\Theta\otimes(V^HA^HU_{\bot}) i(\Theta\otimes(V^HA^HU_{\bot}))\Big], \\ S_{M} =\Big[\big(({{{\rm{her}}}}(V^HA^HU\Theta))^T\otimes I_p\big)K_{p} i\big(({{{\rm{her}}}}(V^HA^HU\Theta))^T\otimes I_p\big)S_{p}\Big]; \end{array} \end{equation} $

(3.22) $ \begin{equation} \begin{array}{ll} T_{E}=\Big[(\Theta\otimes(V_{\bot}^HA^HU))K_{p} i(\Theta\otimes(V_{\bot}^HA^HU))S_{p}\Big], \\ T_{F} =\Big[\Theta\otimes(V_{\bot}^HA^HU_{\bot}) i(\Theta\otimes(V_{\bot}^HA^HU_{\bot}))\Big], \\ T_{N}=\Big[({{{\rm{her}}}}(V^HA^HU\Theta))^T\otimes I_{n-p} i\big(({{{\rm{her}}}}(V^HA^HU\Theta))^T\otimes I_{n-p}\big)\Big]. \end{array} \end{equation} $

证   因为$ E_H\in{\mathbb AH}^{p\times p} $, 由(3.6)式, 结合(3.19)式中$ Q_E, Q_M $和$ Q_N $的定义和(2.25)–(2.29)式, $ {{{\rm{vec}}}}(E_H)\cong {{{\rm{vec}}}}(\widetilde{E_H}) $可计算如下

$ \begin{eqnarray*} && \left[\begin{array}{ccccc} K_p\ & O \\ O\ & S_p \end{array}\right]\left[\begin{array}{ccccc} {\rm vec_K}({\rm Re}(E_H)) \\ {\rm vec_S}({\rm Im}(E_H)) \end{array}\right]\\ &=&\left[\begin{array}{ccccc} { } \frac{1}{2}\left(I_{p^{2}}-T_{(p, p)}\right) & O \\ O & { } \frac{1}{2}\left(I_{p^{2}}+T_{(p, p)}\right) \end{array}\right]\Bigg\{\left[\begin{array}{ccccc} {\rm Re}(Q_E) \\ {\rm Im}(Q_E) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec_K}({\rm Re}(E)) \\ {\rm vec_S}({\rm Im}(E))\end{array}\right]\\ &&-\left[\begin{array}{ccccc} {\rm Re}(Q_M) \\ {\rm Im}(Q_M) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec_K}({\rm Re}(M)) \\ {\rm vec_S}({\rm Im}(M))\end{array}\right] -\left[\begin{array}{ccccc} {\rm Re}(Q_N) \\ {\rm Im}(Q_N) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(N)) \\ {\rm vec}({\rm Im}(N))\end{array}\right]\Bigg\}. \end{eqnarray*} $

注意到$ K_p^TK_p=I_{{n(n-1)}/{2}} $, $ S_p^TS_p=I_{{n(n+1)}/{2}} $, 结合(3.15)式中$ Q_{E, 1} $, $ Q_{E, 2} $, $ Q_{F, 1} $, $ Q_{F, 2} $, $ Q_{M, 1} $, $ Q_{M, 2} $, $ Q_{N, 1} $和$ Q_{N, 2} $的定义可得

$ \begin{eqnarray*} &&\left[\begin{array}{ccccc} {\rm vec}_{K}({\rm Re}(E_H)) \\ {\rm vec}_{S}({\rm Im}(E_H)) \end{array}\right]\\ & =&\left[\begin{array}{ccccc} Q_{E, 1} \\ Q_{E, 2} \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}_{K}({\rm Re}(E)) \\ {\rm vec}_{S}({\rm Im}(E)) \end{array}\right] +\left[\begin{array}{ccccc} Q_{M, 1} \\ Q_{M, 2} \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}_{K}({\rm Re}(M)) \\ {\rm vec}_{S}({\rm Im}(M)) \end{array}\right]+\left[\begin{array}{ccccc} Q_{N, 1} \\ Q_{N, 2} \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(N)) \\ {\rm vec}({\rm Im}(N)) \end{array}\right]. \end{eqnarray*} $

由(3.7)式, 结合(2.25)–(2.29)式以及(3.16)式和(3.20)式, $ {\rm vec}(F_H) \cong {\rm vec}(\widetilde{F_H}) $可计算如下

$ \begin{eqnarray*} &&\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(F_H)) \\ {\rm vec}({\rm Im}(F_H)) \end{array}\right] \\ &=& \left[\begin{array}{ccccc} {\rm Re}(R_F) \\ {\rm Im}(R_F) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(F)) \\ {\rm vec}({\rm Im}(F)) \end{array}\right] -\left[\begin{array}{ccccc} {{{\rm{Re}}}}(R_M) \\ {{{\rm{Im}}}}(R_M) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}_{K}({\rm Re}(M)) \\ {\rm vec}_{S}({\rm Im}(M)) \end{array}\right]-\left[\begin{array}{ccccc} {{{\rm{Re}}}}(R_N) \\ {{{\rm{Im}}}}(R_N) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(N)) \\ {\rm vec}({\rm Im}(N)) \end{array}\right]\\ &=&\left[\begin{array}{ccccc} R_{F, 1} \\ R_{F, 1} \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(F)) \\ {\rm vec}({\rm Im}(F)) \end{array}\right] +\left[\begin{array}{ccccc} R_{M, 1} \\ R_{M, 1} \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}_{K}({\rm Re}(M)) \\ {\rm vec}_{S}({\rm Im}(M)) \end{array}\right]+\left[\begin{array}{ccccc} R_{N, 1} \\ R_{N, 1} \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(N)) \\ {\rm vec}({\rm Im}(N)) \end{array}\right]. \end{eqnarray*} $

类似于$ {{{\rm{vec}}}}(E_H) $和$ {{{\rm{vec}}}}(F_H) $的推导, 注意到$ M_H\in{\mathbb AH}^{p\times p} $, 结合(3.8), (3.9), (2.25)–(2.29), (3.17), (3.18), (3.21)式和(3.22)式可得

$ \begin{eqnarray*} &&\left[\begin{array}{ccccc} {\rm vec_K}({\rm Re}(M_H)) \\ {\rm vec_S}({\rm Im}(M_H)) \end{array}\right] \\ &=&\left[\begin{array}{ccccc} { } \frac{1}{2}K_p^T\left(I_{p^{2}}-T_{(p, p)}\right) & O \\ O &{ } \frac{1}{2}S_p^T\left(I_{p^{2}}+T_{(p, p)}\right) \end{array}\right]\Bigg\{-\left[\begin{array}{ccccc} {\rm Re}(S_E) \\ {\rm Im}(S_E) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec_K}({\rm Re}(E)) \\ {\rm vec_S}({\rm Im}(E))\end{array}\right]\\ &&-\left[\begin{array}{ccccc} {\rm Re}(S_F) \\ {\rm Im}(S_F) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(F)) \\ {\rm vec}({\rm Im}(F))\end{array}\right]+\left[\begin{array}{ccccc} {\rm Re}(S_M) \\ {\rm Im}(S_M) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec_K}({\rm Re}(M)) \\ {\rm vec_S}({\rm Im}(M))\end{array}\right]\Bigg\}\\ &=& \left[\begin{array}{ccccc} S_{E, 1} \\ S_{E, 2} \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}_{K}({\rm Re}(E)) \\ {\rm vec}_{S}({\rm Im}(E)) \end{array}\right] +\left[\begin{array}{ccccc} S_{F, 1} \\ S_{F, 2} \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(F)) \\ {\rm vec}({\rm Im}(F)) \end{array}\right]+\left[\begin{array}{ccccc} S_{M, 1} \\ S_{M, 2} \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}_K({\rm Re}(M)) \\ {\rm vec}_S({\rm Im}(M)) \end{array}\right]. \end{eqnarray*} $

$ \begin{eqnarray*} &&\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(N_H)) \\ {\rm vec}({\rm Im}(N_H)) \end{array}\right] \\ &= &-\left[\begin{array}{ccccc} {{{\rm{Re}}}}(T_E) \\ {{{\rm{Im}}}}(T_E) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}_K({\rm Re}(E)) \\ {\rm vec}_S({\rm Im}(E)) \end{array}\right] -\left[\begin{array}{ccccc} {{{\rm{Re}}}}(T_F) \\ {{{\rm{Im}}}}(T_F) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(F)) \\ {\rm vec}({\rm Im}(F)) \end{array}\right]+\left[\begin{array}{ccccc} {{{\rm{Re}}}}(T_N) \\ {{{\rm{Im}}}}(T_N) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(N)) \\ {\rm vec}({\rm Im}(N)) \end{array}\right]\\ &=&\left[\begin{array}{ccccc} T_{E, 1} \\ T_{E, 1} \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}_K({\rm Re}(E)) \\ {\rm vec}_S({\rm Im}(E)) \end{array}\right] +\left[\begin{array}{ccccc} T_{F, 1} \\ T_{F, 1} \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}_{K}({\rm Re}(F)) \\ {\rm vec}_{S}({\rm Im}(F)) \end{array}\right]+\left[\begin{array}{ccccc} T_{N, 1} \\ T_{N, 1} \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}({\rm Re}(N)) \\ {\rm vec}({\rm Im}(N)) \end{array}\right]. \end{eqnarray*} $

引理3.3得证.

推论3.1   在黎曼度量(2.3)下, 矩阵$ H_A $是对称的.

证   $ \forall (\xi, \eta), (\xi', \eta')\in {\bf{T}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $, 由(3.4)式得, $ (\xi, \eta), (\xi', \eta') $可分别表示为

$ (\xi, \eta)=(UE+U_{\bot}F, VM+V_{\bot}N), \ \ \ \ (\xi', \eta')=(UE'+U_{\bot}F', VM'+V_{\bot}N'), $

其中$ E, M, E', M'\in {\mathbb AH}^{p\times p}, F, F'\in {\mathbb C}^{(m-p)\times p}, N, N'\in {\mathbb C}^{(n-p)\times p} $. 结合(2.3), (2.23)式和(2.24)式, 可得

(3.23) $ \begin{eqnarray} \langle (\xi, \eta), \ (\varepsilon', \eta')\rangle &=&{{{\rm{Tr}}}}(\xi^H\xi')+{{{\rm{Tr}}}} (\eta^H\eta'){}\\ &=&{{{\rm{Tr}}}}(E^HE')+{{{\rm{Tr}}}}(F^HF')+{{{\rm{Tr}}}}(M^HM')+{{{\rm{Tr}}}}(N^HN'){}\\ &=&{{{\rm{tr}}}}({{{\rm{Re}}}}(E)^T{{{\rm{Re}}}}(E'))+{{{\rm{tr}}}}({{{\rm{Im}}}}(E)^T{{{\rm{Im}}}}(E'))+{{{\rm{tr}}}}({{{\rm{Re}}}}(F)^T{{{\rm{Re}}}}(F')){}\\ &&+{{{\rm{tr}}}}({{{\rm{Im}}}}(F)^T{{{\rm{Im}}}}(F'))+{{{\rm{tr}}}}({{{\rm{Re}}}}(M)^T{{{\rm{Re}}}}(M'))+{{{\rm{tr}}}}({{{\rm{Im}}}}(M)^T{{{\rm{Im}}}}(M')){}\\ &&+{{{\rm{tr}}}}({{{\rm{Re}}}}(N)^T{{{\rm{Re}}}}(N'))+{{{\rm{tr}}}}({{{\rm{Im}}}}(N)^T{{{\rm{Im}}}}(N')){}\\ &=&{{{\rm{vec}}}}_K({{{\rm{Re}}}}(E))^T{{{\rm{vec}}}}_K({{{\rm{Re}}}}(E'))+{{{\rm{vec}}}}_S({{{\rm{Im}}}}(E))^T{{{\rm{vec}}}}_S({{{\rm{Im}}}}(E')){}\\ &&+{{{\rm{vec}}}}({{{\rm{Re}}}}(F))^T{{{\rm{vec}}}}({{{\rm{Re}}}}(F'))+{{{\rm{vec}}}}({{{\rm{Im}}}}(F))^T{{{\rm{vec}}}}({{{\rm{Im}}}}(F')){}\\ &&+ {{{\rm{vec}}}}_K({{{\rm{Re}}}}(M))^T{{{\rm{vec}}}}_K({{{\rm{Re}}}}(M'))+{{{\rm{vec}}}}_S({{{\rm{Im}}}}(N))^T{{{\rm{vec}}}}_S({{{\rm{Im}}}}(N')){}\\ &&+{{{\rm{vec}}}}({{{\rm{Re}}}}(N))^T{{{\rm{vec}}}}({{{\rm{Re}}}}(N'))+{{{\rm{vec}}}}({{{\rm{Im}}}}(N))^T{{{\rm{vec}}}}({{{\rm{Im}}}}(N')). \end{eqnarray} $

因为黎曼海塞是对称变换^[13], 即$ {{{\rm{Hess }}}} F(U, V) $在黎曼度量(2.3)下是对称的, 即

$ \langle(\xi, \eta), {{{\rm{Hess }}}} F(U, V)\left[\left(\xi^{\prime}, \eta^{\prime}\right)\right]\rangle=\langle{\rm Hess} F(U, V)[(\xi, \eta)], \left(\xi^{\prime}, \eta^{\prime}\right)\rangle. $

结合引理3.3, 上式等价于

$ \left[\begin{array}{ccccc} {{{\rm{vec}}}}_K({{{\rm{Re}}}}(E))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(E))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(F))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(F))\\ {{{\rm{vec}}}}_K({{{\rm{Re}}}}(M))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(M))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(N))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(N)) \end{array}\right]^TH_A\left[\begin{array}{ccccc} {{{\rm{vec}}}}_K({{{\rm{Re}}}}(E'))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(E'))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(F'))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(F'))\\ {{{\rm{vec}}}}_K({{{\rm{Re}}}}(M'))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(M'))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(N'))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(N')) \end{array}\right]=\left(H_A\left[\begin{array}{ccccc} {{{\rm{vec}}}}_K({{{\rm{Re}}}}(E))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(E))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(F))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(F))\\ {{{\rm{vec}}}}_K({{{\rm{Re}}}}(M))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(M))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(N))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(N)) \end{array}\right] \right)^T\left[\begin{array}{ccccc} {{{\rm{vec}}}}_K({{{\rm{Re}}}}(E'))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(E'))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(F'))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(F'))\\ {{{\rm{vec}}}}_K({{{\rm{Re}}}}(M'))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(M'))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(N'))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(N')) \end{array}\right]. $

进而由$ (\xi, \eta) $和$ (\xi', \eta') $的任意性可得$ H_A=H_A^T. $推论3.1得证.

接下来求解牛顿方程(3.3). 由(3.5)式, 牛顿方程可以改写成

(3.24) $ \begin{eqnarray} &&UE_H+U_{\bot}F_H={\bf P}_U(AV\Theta), \end{eqnarray} $

(3.25) $ \begin{eqnarray} &&VM_H+V_{\bot}N_H={\bf P}_V(A^HU\Theta). \end{eqnarray} $

在(3.24)式两边分别左乘$ U^H $和$ U_{\bot}^H $, 在(3.25)式两边分别左乘$ V^H $和$ V_{\bot}^H $, 得到

(3.26) $ \begin{equation} \begin{array}{rl} &E_H=U^H{\bf P}_U(AV\Theta), F_H=U_{\bot}^H{\bf P}_U(AV\Theta), \\ &M_H=V^H{\bf P}_V(A^HUN), N_H=V_{\bot}^H{\bf P}_V(A^HUN). \end{array} \end{equation} $

再$ \forall W_1\in {\mathbb C}^{m\times p} $, $ U^H {\bf{P}}_{U}(W_1)={{{\rm{skewH}}}}(U^HW_1) $和$ U_{\bot}^H {\bf{P}}_U(W_1)=U_{\bot}^HW_1 $; $ \forall W_2\in {\mathbb C}^{n\times p} $, $ V^H {\bf{P}}_{V}(W_2)={{{\rm{skewH}}}}(V^HW_2) $和$ V_{\bot}^H {\bf{P}}_V(W_2)=V_{\bot}^HW_2 $. 应用拉直算子$ {\rm vec} $作用于(3.26)式, 结合(2.25)–(2.29)式可得

$ \begin{eqnarray*} &&\left[\begin{array}{ccccc} {\rm vec}_{K}\left({\rm Re}\left(E_{H}\right)\right) \\ {\rm vec}_{S}\left({\rm Im}\left(E_{H}\right)\right) \end{array}\right]=\left[\begin{array}{ccccc} { } \frac{1}{2} K_{p}^{T}\left(I_{p^{2}}-T_{(p, p)}\right) & O \\ O& { } \frac{1}{2} S_{p}^{T}\left(I_{p^{2}}+T_{(p, p)}\right) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}\left({\rm Re}\left(U^HAV\Theta\right)\right. \\ {\rm vec}\left({\rm Im}\left(U^HAV\Theta\right)\right. \end{array}\right], \\ &&\left[\begin{array}{ccccc} {\rm vec}\left({\rm Re}\left(F_{H}\right)\right) \\ {\rm vec}\left({\rm Im}\left(F_{H}\right)\right) \end{array}\right]=\left[\begin{array}{ccccc} {\rm vec}\left({\rm Re}\left(U_{\bot}^HAV\Theta\right)\right. \\ {\rm vec}\left({\rm Im}\left(U_{\bot}^HAV\Theta\right)\right. \end{array}\right], \\ &&\left[\begin{array}{ccccc} {\rm vec}_{K}\left({\rm Re}\left(M_{H}\right)\right) \\ {\rm vec}_{S}\left({\rm Im}\left(M_{H}\right)\right) \end{array}\right]=\left[\begin{array}{ccccc} { } \frac{1}{2} K_{p}^{T}\left(I_{p^{2}}-T_{(p, p)}\right) & O \\ O &{ } \frac{1}{2} S_{p}^{T}\left(I_{p^{2}}+T_{(p, p)}\right) \end{array}\right]\left[\begin{array}{ccccc} {\rm vec}\left({\rm Re}\left(V^HA^HU\Theta\right)\right. \\ {\rm vec}\left({\rm Im}\left(V^HA^HU\Theta\right)\right. \end{array}\right], \\ &&\left[\begin{array}{ccccc} {\rm vec}\left({\rm Re}\left(N_{H}\right)\right) \\ {\rm vec}\left({\rm Im}\left(N_{H}\right)\right) \end{array}\right]=\left[\begin{array}{ccccc} {\rm vec}\left({\rm Re}\left(V_{\bot}^HA^HU\Theta\right)\right. \\ {\rm vec}\left({\rm Im}\left(V_{\bot}^HA^HU\Theta\right)\right. \end{array}\right]. \end{eqnarray*} $

最后由引理3.3, 牛顿方程(3.3)可转化为如下线性方程组形式

(3.27) $ \begin{equation} H_Ax=g, \end{equation} $

其中$ H_A $由(3.14)式给出,

(3.28) $ \begin{equation} x=\left[\begin{array}{ccccc} {\rm vec}_{K}\left({\rm Re}\left(E\right)\right) \\ {\rm vec}_{S}\left({\rm Im}\left(E\right)\right) \\ {\rm vec}\left({\rm Re}\left(F\right)\right) \\ {\rm vec}\left({\rm Im}\left(F\right)\right) \\ {\rm vec}_{K}\left({\rm Re}\left(M\right)\right) \\ {\rm vec}_{S}\left({\rm Im}\left(M\right)\right) \\ {\rm vec}\left({\rm Re}\left(N\right)\right) \\ {\rm vec}\left({\rm Im}\left(N\right)\right) \end{array}\right], \ \ \ \ \ g=\left[\begin{array}{ccccc} { } \frac{1}{2} K_{p}^{T}\left(I_{p^{2}}-T_{(p, p)}\right) {\rm vec}\left({\rm Re}\left(U^HAV\Theta\right)\right)\\ { } \frac{1}{2} S_{p}^{T}\left(I_{p^{2}}+T_{(p, p)}\right) {\rm vec}\left({\rm Im}\left(U^HAV\Theta\right)\right)\\ {\rm vec}\left({\rm Re}\left(U_{\bot}^HAV\Theta\right)\right) \\ {\rm vec}\left({\rm Im}\left(U_{\bot}^HAV\Theta\right)\right)\\ { } \frac{1}{2} K_{p}^{T}\left(I_{p^{2}}-T_{(p, p)}\right) {\rm vec}\left({\rm Re}\left(V^HA^HU\Theta\right)\right)\\ { } \frac{1}{2} S_{p}^{T}\left(I_{p^{2}}+T_{(p, p)}\right) {\rm vec}\left({\rm Im}\left(V^HA^HU\Theta\right)\right)\\ {\rm vec}\left({\rm Re}\left(V_{\bot}^HA^HU\Theta\right)\right)\\ {\rm vec}\left({\rm Im}\left(V_{\bot}^HA^HU\Theta\right)\right)\\ \end{array}\right]. \end{equation} $

因变换矩阵$ H_A $对称, 故(3.27)式为实对称线性方程组. 求解该实对称线性方程组, 可得$ {\rm vec}_{K}\left({\rm Re}\left(E\right)\right) $, $ {\rm vec}_{S}\left({\rm Im}\left(E\right)\right) $, $ {\rm vec}\left({\rm Re}\left(F\right)\right) $, $ {\rm vec}\left({\rm Im}\left(F\right)\right) $, $ {\rm vec}_{K}\left({\rm Re}\left(M\right)\right) $, $ {\rm vec}_{S}\left({\rm Im}\left(M\right)\right) $, $ {\rm vec}\left({\rm Re}\left(N\right)\right) $和$ {\rm vec}\left({\rm Re}\left(N\right)\right) $. 进一步通过拉直算子逆运算可得$ E={\rm Re}\left(E\right)+{\rm i} {\rm Im}\left(E\right) \in {\mathbb AH}^{p \times p} $, $ F={\rm Re}\left(F\right)+{\rm i} {\rm Im}\left(F\right) \in {\mathbb C}^{(m-p) \times p} $, $ M={\rm Re}\left(M\right)+{\rm i} {\rm Im}\left(M\right) \in {\mathbb AH}^{p \times p} $和$ N={\rm Re}\left(N\right)+{\rm i} {\rm Im}\left(N\right) \in {\mathbb C}^{(n-p) \times p}. $从而得到牛顿方程(3.3)的解, 即搜索方向$ (\xi, \eta)=\left(U E+U_{\bot} F, V M+V_{\bot} N\right) $.

基于以上讨论, 求解问题(1.2)的黎曼牛顿法的完整迭代格式如下:

算法1   求解问题(1.2)的黎曼牛顿法(Newton)

1. 给定初始点$ (U^0, V^0)\in {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $.

2. 对$ k=0, 1, 2, \cdots, $ 执行下面的循环

3. 计算$ U^{(k)}_{\bot} $和$ V^{(k)}_{\bot} $满足$ {U^{(k)}}^{H} U_{\bot}^{(k)}=O $, $ {U_{\bot}^{(k)}}^{H}U_{\bot}^{(k)}=I_{m-p} $, $ {V^{(k)}}^{H} V_{\bot}^{(k)}=O $和$ {V_{\perp}^{(k)}}^{H}V_{\perp}^{(k)}=I_{n-p}. $

4. 由引理3.3和(3.28)分别计算$ H_A^{(k)} $和$ g^{(k)} $. 计算$ x^{(k)}\in{{\mathbb R}} ^{\kappa} $满足如下实对称线性方程组

(3.29) $ \begin{equation} H_A^{(k)}x^{(k)}=g^{(k)}. \end{equation} $

5. 将$ x^{(k)} $分块为$ {x^{(k)}}^T=\left[{x_1^{(k)}}^T {x_2^{(k)}}^T {x_3^{(k)}}^T {x_4^{(k)}}^T {x_5^{(k)}}^T {x_6^{(k)}}^T {x_7^{(k)}}^T {x_8^{(k)}}^T\right] $, 其中$ x_{1}^{(k)}\in{{\mathbb R}} ^{{p(p-1)}/{2}}, $$ x_{2}^{(k)}\in{{\mathbb R}} ^{{p(p+1)}/{2}}, $$ x_{3}^{(k)}\in{{\mathbb R}} ^{(m-p)\times p}, $$ x_{4}^{(k)}\in{{\mathbb R}} ^{(m-p)\times p}, $$ x_{5}^{(k)}\in{{\mathbb R}} ^{{p(p-1)}/{2}}, $$ x_{6}^{(k)}\in{{\mathbb R}} ^{{p(p+1)}/{2}}, $$ x_{7}^{(k)}\in{{\mathbb R}} ^{(n-p)\times p}, $$ x_{8}^{(k)}\in{{\mathbb R}} ^{(n-p)\times p}. $

6. 重构$ E, M, \in {\mathbb AH}^{p\times p}, F\in {\mathbb C}^{(m-p)\times p}, N\in {\mathbb C}^{(n-p)\times p} $使其满足

$ \begin{array}{lll} {{{\rm{vec}}}}_K\big({{{\rm{Re}}}}(E^{(k)})\big)=x_1^{(k)}, \ & {{{\rm{vec}}}}_S\big({{{\rm{Im}}}}(E^{(k)})\big)=x_2^{(k)}, \ & E^{(k)}={{{\rm{Re}}}}(E^{(k)})+{\rm i}{{{\rm{Im}}}}(E^{(k)});\\ {{{\rm{vec}}}}\big({{{\rm{Re}}}}(F^{(k)})\big)=x_3^{(k)}, \ & {{{\rm{vec}}}}\big({{{\rm{Im}}}}(F^{(k)})\big)=x_4^{(k)}, \ & F^{(k)}={{{\rm{Re}}}}(F^{(k)})+{\rm i}{{{\rm{Im}}}}(F^{(k)});\\ {{{\rm{vec}}}}_K\big({{{\rm{Re}}}}(M^{(k)})\big)=x_5^{(k)}, \ & {{{\rm{vec}}}}_S\big({{{\rm{Im}}}}(M^{(k)})\big)=x_6^{(k)}, \ & M^{(k)}={{{\rm{Re}}}}(M^{(k)})+{\rm i}{{{\rm{Im}}}}(M^{(k)});\\ {{{\rm{vec}}}}\big({{{\rm{Re}}}}(N^{(k)})\big)=x_7^{(k)}, \ & {{{\rm{vec}}}}\big({{{\rm{Im}}}}(N^{(k)})\big)=x_8^{(k)}, \ & N^{(k)}={{{\rm{Re}}}}(N^{(k)})+{\rm i}{{{\rm{Im}}}}(N^{(k)}). \end{array} $

7. 计算$ \xi^{(k)}=U^{(k)}E^{(k)}+U_{\bot}^{(k)}F^{(k)}, $$ \eta^{(k)}=V^{(k)}M^{(k)}+V_{\bot}^{(k)}N^{(k)} $.

8. 计算更新迭代步

$ (U^{(k+1)}, V^{(k+1)}) ={\mathcal R}_{(U^{(k)}, V^{(k)})}\left(\xi^{(k)}, \eta^{(k)}\right). $

9. 循环结束.

注3.1   算法Newton中第3步的$ U_{\perp}^{(k)} $和$ V_{\perp}^{(k)} $可通过QR分解来实现. 实际计算中, 利用MATLAB中的qr函数作用于$ U^{(k)} \in {{{\rm{St}}}}(m, p, {\mathbb C}) $可得$ U^{(k)}=Q \left[\begin{array}{ccccc}{R_{1}} \\ {O} \end{array}\right] $, 其中$ Q \in {\mathbb C}^{m \times m} $满足$ Q^HQ=I_m $, $ R_{1} $为$ p \times p $为上三角阵. 令$ Q=\left[\begin{array}{ccccc}Q_{1}\ Q_{2}\end{array}\right] $, $ Q_{1} \in {\rm St}(m, p, {\mathbb C}) $, $ Q_{2} \in {\rm St}(m, m-p, {\mathbb C}) $, 则有$ Q_{2}^{T} U^{(k)}=\left[\begin{array}{ccccc}O\ & I_{m-p}\end{array}\right] \left[\begin{array}{ccccc}{R_{1}} \\ {O}\end{array}\right]=O $. 故可取$ Q_{2}=U_{\perp}^{(k)} $. 类似地可计算$ V_{\perp}^{(k)} $.

注3.2  对于算法Newton中的实对称线性方程组$ H_{A}^{(k)} x^{(k)}=g^{(k)} $, 其中$ H_{A}^{(k)}\in {\mathbb S}^{\kappa\times \kappa} $, 若$ H_{A}^{(k)} $是可逆的则可直接求证. 当问题维度较大时, 直接求逆比较困难, 此时可通过Krylov子空间方法中的若干方法求解, 如共轭残量法(Conjugate Residual method, CR)^[18]. 因对于求解线性方程组$ Hx=g $, 其中$ H $是对称可逆矩阵且$ g $非零, 共轭残量法只要求系数矩阵对称, 不要求对称正定. 若$ H_{A}^{(k)} $进一步对称正定, 则可采用共轭梯度法(CG)求解. 另外, 因对未知向量$ x $没有附加的结构要求, 若将(3.29)式视为一般线性方程组, 则可采用迭代效率更高的广义极小残差(GMRES) 或稳定双共轭梯度(BiCGSTAB)等方法求解.

注3.3   一般来说, 由(3.1)式得到的搜索方向$ (\xi^{(k)}, \eta^{(k)}) $并不一定是$ F(U, V) $的下降方向. 实际上

$ {\mathcal D} F(U, V)[(\xi, \eta)]=\langle {{{\rm{grad }}}} F(U, V), (\xi, \eta)\rangle=-\langle {{{\rm{grad }}}} F(U, V), {{{\rm{Hess }}}} F(U, V)^{-1}{{{\rm{grad }}}} F(U, V)\rangle. $

由上式可知, 若要使得$ (\xi, \eta) $为下降方向, 一个可行的充分条件是黎曼海塞$ {{{\rm{Hess }}}} F(U, V)[(\xi, \eta)] $对称正定, 即$ \forall (\xi, \eta)\neq 0 $, $ \langle (\xi, \eta), {{{\rm{Hess }}}} F(U, V)[(\xi, \eta)]\rangle>0 $. 注意到

$ \forall (\xi, \eta)\in {\bf{T}}_{(U, V)}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}), (\xi, \eta)=(UE+U_{\bot}F, VM+V_{\bot}N), $

其中$ E, M\in {\mathbb AH}^{p\times p}, F\in {\mathbb C}^{(m-p)\times p}, N\in {\mathbb C}^{(n-p)\times p} $, 有

$ \langle{{{\rm{Hess }}}} F(U, V)[(\xi, \eta)], \ (\xi, \eta)\rangle= \left[\begin{array}{ccccc} {{{\rm{vec}}}}_K({{{\rm{Re}}}}(E))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(E))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(F))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(F))\\ {{{\rm{vec}}}}_K({{{\rm{Re}}}}(M))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(M))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(N))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(N)) \end{array}\right]^TH_A\left[\begin{array}{ccccc} {{{\rm{vec}}}}_K({{{\rm{Re}}}}(E))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(E))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(F))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(F))\\ {{{\rm{vec}}}}_K({{{\rm{Re}}}}(M))\\ {{{\rm{vec}}}}_S({{{\rm{Im}}}}(M))\\ {{{\rm{vec}}}}({{{\rm{Re}}}}(N))\\ {{{\rm{vec}}}}({{{\rm{Im}}}}(N)) \end{array}\right], $

故可采用变换矩阵$ H_A $的正定性来判定黎曼海塞$ {{{\rm{Hess }}}} F(U, V)[(\xi, \eta)] $的正定性, 进而可判定由算法Newton得到的驻点是否为问题(1.2) 的局部极小点.

注3.4   类似于无约束优化, $ \left\|{{{{\rm{grad }}}}}F(U^{(k)}, V^{(k)})\right\| \leq \varepsilon $ ($ \varepsilon>0 $) 可以作为算法Newton的迭代终止条件. 事实上, 问题(1.2) 对应的拉格朗日函数为

$ {\mathcal L}:={\mathcal L}( {U}, {V}, {\Lambda}, {\Omega})=\bar{F}( {U}, {V})-\frac{1}{2}{{{\rm{Tr}}}}( {\Lambda}( {U}^H {U}- {I}_p))-\frac{1}{2}{{{\rm{Tr}}}}( {\Omega}( {V}^H {V}- {I}_p)), $

其中$ {\Lambda}, {\Omega} $分别为约束条件$ U^HU=I_p $和$ V^HV=I_p $所对应的Hermite拉格朗日乘子. 由此可知问题(1.2) 的一阶最优性条件为

$ {U}^H {U}= {I}_p, \ \ \ \ \frac{\partial}{\partial {U}}{\mathcal L}=\nabla_U\bar{F}- {U}{\nabla_U\bar{F}}^H {U}=O, $

$ {V}^H {V}= {I}_p, \ \ \ \ \frac{\partial}{\partial {V}}{\mathcal L}=\nabla_V\bar{F}- {V}{\nabla_V\bar{F}}^H {V}=O, $

其中$ {\Lambda}=\nabla_U\bar{F}^H {U} $和$ {\Omega}=\nabla_V\bar{F}^H {V} $. 由$ {U}^H {U}= {I}_p $和$ {V}^H {V}= {I}_p $, $ \frac{\partial}{\partial {U}}{\mathcal L}=O $和$ \frac{\partial}{\partial {V}}{\mathcal L}=O $等价于$ {{{\rm{grad }}}} f( {U}, {V})=O $, 由(2.10), (2.4)式和(2.5)式有

$ {\bf P}_{ {U}}\Big(\nabla_U\bar{F}\Big)= ( {I}_p-\frac{1}{2} {U} {U}^H)\Big( \nabla_U\bar{F}- {U}{\nabla_U\bar{F}}^H {U}\Big)=( {I}_p-\frac{1}{2} {U} {U}^H)\frac{\partial}{\partial {U}}{\mathcal L}, $

$ {\bf P}_{ {V}}\Big(\nabla_V\bar{F}\Big)= ( {I}_p-\frac{1}{2} {V} {V}^H)\Big( \nabla_V\bar{F}- {V}{\nabla_V\bar{F}}^H {V}\Big)=( {I}_p-\frac{1}{2} {V} {V}^H)\frac{\partial}{\partial {V}}{\mathcal L}. $

一般来说牛顿法不能保证全局收敛, 且对迭代初值的选取比较敏感. 这里我们应用黎曼优化中的经典一阶梯度类算法OptStiefelGBB^[12]来产生一个合适的迭代初值. 该一阶算法的优势在于单步迭代成本低且存储量少. 对于迭代框架(2.1), 算法OptStiefelGBB利用非单调Armijo原则来得到步长$ t_k $, 即给定$ \rho, \varrho, \delta \in(0, 1) $, $ t_{k}=\gamma_{k} \delta^{h} $, 其中$ \gamma_{k} $为初始步长, $ h $为满足以下不等式的最小正整数

(3.30) $ \begin{equation} F\left({\mathcal R}_{(U^{(k)}, V^{(k)})}\left(\xi^{(k)}, \eta^{(k)}\right)\right) \leq C_{k}+\rho t_{k}\langle{{{{\rm{grad }}}}} F(U^{(k)}, V^{(k)}), (\xi^{(k)}, \eta^{(k)})\rangle. \end{equation} $

为加速收敛, 初始步长$ \gamma_{k} $通过Barzilai-Borwein (BB)方法来确定, 即

(3.31) $ \begin{eqnarray} &&\gamma_{k}^{(1)}=\frac{{{{\rm{Tr}}}}(S_{U, k-1}^{H} S_{U, k-1})+{{{\rm{Tr}}}}(S_{V, k-1}^{H} S_{V, k-1})}{ \left| {{{{\rm{Tr}}}}}(Z_{U, k-1}^{H} S_{U, k-1})+ {{{{\rm{Tr}}}}}(Z_{V, k-1}^{H} S_{V, k-1})\right|} {}\\ \mbox{ 或 }&& \gamma_{k}^{(2)}=\frac{\left|{{{{\rm{Tr}}}}}(S_{U, k-1}^{H} Z_{U, k-1})+{{{{\rm{Tr}}}}}(S_{V, k-1}^{H} Z_{V, k-1})\right|}{{{{{\rm{Tr}}}}}(Z_{U, k-1}^{H} Z_{U, k-1})+{{{{\rm{Tr}}}}}(Z_{V, k-1}^{H} Z_{V, k-1})}, \end{eqnarray} $

其中$ (S_{U, k-1}, S_{V, k-1})=(U^{(k)}, V^{(k)})-(U^{(k-1)}, V^{(k-1)}) $, $ (Z_{V, k-1}, Z_{U, k-1})={{{{\rm{grad }}}}} F(U^{(k)}, V^{(k)})-{{{{\rm{grad }}}}} F(U^{(k-1)}, V^{(k-1)}). $为保证全局收敛性, (3.30)式中参考函数值$ C_{k} $由以下方式确定

$ C_{k}=\frac{1}{\theta_{k}}\left(\varrho \theta_{k-1} C_{k-1}+F(U^{(k)}, V^{(k)})\right), $

其中$ C_{0}=F(U^{(0)}, V^{(0)}), \theta_{k}=\varrho \theta_{k-1}+1 $, $ \theta_{0}=1 $. 算法OptStiefelGBB求解问题(1.2)的完整迭代格式如下

算法2   求解问题(1.2)的黎曼梯度类算法(OptStiefelGBB)

1. 任取初始值$ (U^{(0)}, V^{(0)})\in {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $. 令$ k=0, \gamma_{\min } \in[0, 1], \gamma_{\max } \geq 1 $. $ C_{0}=F(U^{(0)}, V^{(0)}), \theta_{0}=1 $.

2. 对$ k=0, 1, 2, \cdots $执行下面的循环

3. 计算搜索方向$ (\xi^{(k)}, \eta^{(k)})\in {\bf{T}}_{(U^{(k)}, V^{(k)})}{{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $,

$ (\xi^{(k)}, \eta^{(k)}) =\left( {\bf{P}}_{U^{(k)}}(-AV^{(k)}\Theta), {\bf{P}}_{V^{(k)}}(-A^HU^{(k)}\Theta)\right). $

4. 根据(3.31)式计算$ \gamma_{k} $, 取$ \gamma_{k}=\max \left\{\gamma_{\min }, \min \left\{\gamma_{k}, \gamma_{\max }\right\}\right\}. $计算$ C_{k}, \theta_{k} $并且找到满足(3.30) 式的步长$ t_{k} $.

5. 计算迭代更新步

$ (U^{(k+1)}, V^{(k+1)}) ={\mathcal R}_{(U^{(k)}, V^{(k)})}\left(t_k(\xi^{(k)}, \eta^{(k)})\right). $

6. 循环结束.

为加速收敛, 我们提出一类基于黎曼牛顿法的混合算法, 其迭代框架如下

算法3   求解问题(1.2)的黎曼牛顿混合算法(Hybrid-Newton)

1. 任取初始值$ (U^{(0)}, V^{(0)})\in {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $, 参数$ \varepsilon_1>\varepsilon>0 $. 令$ k=0 $.

2. 对$ \left\|{{{{\rm{grad }}}}} F(U^{(k)}, V^{(k)})\right\| \leq \varepsilon_1 $ 执行下面的循环

3. 执行算法OptStiefelGBB中的第3–5步.

4. $ k:=k+1 $.

5. 循环结束.

6. 令$ (U^{(0)}, V^{(0)}):=(U^{(k)}, V^{(k)}) $和$ k:=0 $.

7. 执行算法Newton的第2–9步.

算法OptStiefelGBB具有全局收敛性^[12], Stiefel流形上黎曼牛顿法具有局部二阶收敛速度, 见文献[13, 定理6.3.2], 故算法Hybrid-Newton集合了两种算法的优点, 即有如下结论

定理3.1   任意给定初始值$ (U^{(0)}, V^{(0)})\in {{{\rm{St}}}}(m, p, {\mathbb C})\times {{{\rm{St}}}}(n, p, {\mathbb C}) $, 若$ \varepsilon_{1}>0 $足够小, 算法Hybrid-Newton得到的序列$ \{(U^{(k)}, V^{(k)})\} $全局收敛且局部二次收敛于问题(1.2)的局部极小点.

本节给出算法Hybrid-Newton求解问题(1.2)的数值结果. 所有的数据均通过MATALB (R2019), Intel$ ^{\circledR} $ Core i5处理器, 3.20GHz的PC中获得. 算法Newton采用基于Cayley变换的收缩算子, 即给定当前迭代步$ (U^{(k)}, V^{(k)}) $对应的搜索方向$ (\xi^{(k)}, \eta^{(k)}) $, 迭代更新步的具体计算公式为

(4.1) $ \begin{eqnarray} (U^{(k+1)}, V^{(k+1)})&=&{\mathcal R}_{(U^{(k)}, V^{(k)})}\left(\xi^{(k)}, \eta^{(k)}\right){}\\ &=&\bigg(U^{(k)}+P_{\xi^{(k)}}\left(I_{2p}-\frac{1}{2} Q_{\xi^{(k)}}^H P_{\xi^{(k)}}\right)^{-1} Q_{\xi^{(k)}}^H U^{(k)}, {}\\ && V^{(k)}+P_{\eta^{(k)}}\left(I_{2p}-\frac{1}{2} Q_{\eta^{(k)}}^H P_{\eta^{(k)}}\right)^{-1} Q_{\eta^{(k)}}^H V^{(k)}\bigg), \end{eqnarray} $

其中

$ P_{\xi^{(k)}}=\left[\Pi_{U^{(k)}} \xi^{(k)} U^{(k)}\right], Q_{\xi^{(k)}}=\left[U -\Pi_{U^{(k)}} \xi^{(k)}\right], \Pi_{U^{(k)}}=I_n-\frac{1}{2} U^{(k)} {U^{(k)}}^H; $

$ P_{\eta^{(k)}}=\left[\Pi_{V^{(k)}} \eta^{(k)} V^{(k)}\right], Q_{\eta^{(k)}}=\left[V -\Pi_{V^{(k)}} \eta^{(k)}\right], \Pi_{V^{(k)}}=I_n-\frac{1}{2} V^{(k)} {V^{(k)}}^H. $

在实际计算中, 算法Newton中实对称线性方程组(3.29)均采用CR方法求解, 其迭代终止标准取为

$ \frac{\|r^{(k)}_L\|_2}{\|g^{(k)}\|_2}\leq 10^{-6}\ \ \ \ \mbox{或}\ \ \ \ \frac{\|x_L^{(k)}-x_{L-1}^{(k)}\|_2}{\|x_L^{(k)}\|_2}\leq 10^{-6}, $

其中$ \|\cdot\|_2 $表示向量2范数, $ r^{(k)}_L $表示$ x^{(k)}_L $对应的残差$ g^{(k)}-H_A^{(k)}x^{(k)}_L $, $ x^{(k)}_L $为算法Newton中第$ k $次迭代对应的第$ L $次近似值. 算法CR的最大迭代步设为1000. 算法Hybrid-Newton中算法转换参数$ \varepsilon_1=0.01 $. 由注3.4, 终止标准均取为

$ \|{{{\rm{grad }}}} f(U^{(k)}, V^{(k)})\|<\varepsilon, $

其中$ \varepsilon=10^{-6} $. 本节所有表格中, CT., IT., Grad. 和Obj. 分别表示达到停机标准时的迭代时间(s), 迭代步, 黎曼梯度范数$ \|{{{\rm{grad }}}} F(U^{(k)}, V^{(k)})\| $和目标函数值$ F(U^{(k)}, V^{(k)}) $. 特别地, 算法Hybrid-Newton中的IT.由两部分组成, 第一部分表示运行算法OptStifelGBB直到生成合适迭代初值所需的迭代步, 第二部分表示转到算法Newton后继续迭代直到满足停机标准所需的迭代步.

为数值比较, 本节将算法Hybrid-Newton与黎曼优化工具箱Manopt^[19]中的黎曼信赖域方法(记为Existing RTR) 和黎曼BFGS方法(记为Existing RBFGS), 以及黎曼梯度类方法进行比较. 值得说明的是文献[20]将欧式空间中的Dai非线性共轭梯度法^[21]推广到实Stiefel流形, 文献[22]将其推广到实乘积流形$ \mbox{Stp}(m, p)\times \mbox{Stp}(n, p) $, 本文将其进一步推广到复乘积流形$ {{{\rm{St}}}}(p, m, {\mathbb C})\times {{{\rm{St}}}}(p, n, {\mathbb C}) $, 并作适当修改应用于求解问题(1.2). 为此先简要给出该黎曼共轭梯度法(记为RCG-Cayley) 的迭代框架. 算法RCG-Cayley的迭代格式为

$ (U^{(k+1)}, V^{(k+1)})={\mathcal R}_{(U^{(k)}, V^{(k)})}\big(t_{k} (\xi^{(k)}, \eta^{(k)})\big), $

其中$ {\mathcal R} $采用基于Cayley变换的收缩计算公式(4.1), 迭代步长$ t_k $采用如下非单调线搜索方式确定

(4.2) $ \begin{eqnarray} F\left({\mathcal R}_{(U^{(k)}, V^{(k)})}(t_{k} (\xi^{(k)}, \eta^{(k)}))\right) &\leq& \max\{F(U^{(k)}, V^{(k)}), \cdots, F({U}^{(k-h_k)}, {V}^{(k-h_k)})\}{}\\ & &+\delta t_k\langle{{{{\rm{grad }}}}} F(U^{(k)}, V^{(k)}), (\xi^{(k)}, \eta^{(k)})\rangle, \end{eqnarray} $

其中$ \delta\in (0, 1) $, $ h_k=\min\{h-1, k\} $, $ h $为正整数. 搜索方向更新方案为

$ (\xi^{(k+1)}, \eta^{(k+1)})=-{{{\rm{grad }}}} F(U^{(k)}, V^{(k)})+\beta_k{\mathcal T}_{t_k(\xi^{(k)}, \eta^{(k)})}(\xi^{(k)}, \eta^{(k)}), $

其中$ {\mathcal T} $为基于微分收缩的向量转移算子(Vector transport), 其定义其具体计算方式见文献[13, 20]. 参数$ \beta_{k+1}\in [0, \beta_{k+1}^D] $, 其中

$ \beta_{k+1}^D=\frac{\|{{{\rm{grad }}}} F(U^{(k+1)}, V^{(k+1)})\|^2} {\max\left\{M, \ \ -\langle {{{{\rm{grad }}}}} F(U^{(k)}, V^{(k)}), \ \ ({\xi}^{k}, {\eta}^{k})\rangle\right\}}, $

这里$ M=\langle {{{\rm{grad }}}} F(U^{(k+1)}, V^{(k+1)}), \ \ {\mathcal T}_{t_k({\xi}^{k}, {\eta}^{k})}({\xi}^{k}, {\eta}^{k})\rangle -\langle {{{\rm{grad }}}} F(U^{(k)}, V^{(k)}), \ ({\xi}^{k}, {\eta}^{k})\rangle $. 在实际计算中, 算法RCG-Cayley 的具体参数选取参考文献[20], 其停机标准取为$ \|{{{\rm{grad }}}} F(U^{(k)}, V^{(k)})\| $$ \leq 10^{-6} $.

例4.1   本例说明算法Hybrid-Newton的全局收敛性和算法Newton中CR方法对于求解实对称线性组(3.29) 的有效性. 取定$ m=200, n=20, p=5 $. 参考文献[7], 按如下方式生成问题(1.2)中的矩阵$ A\in {\mathbb C}^{m\times n} $和$ \Theta\in {{\mathbb R}}^{p\times p} $

$ A = U_r \mbox{diag}(\sigma_1, \cdots, \sigma_n)V_r^H, \Theta=\mbox{diag}(p, \cdots, 1), $

其中$ U_r\in {\mathbb C}^{m\times n} $和$ V_r\in {\mathbb C}^{n\times n} $为任意生成的(列)酉阵, 即满足$ U^HU=I_n $和$ V^HV=I_n $, $ \sigma_1\geq \cdots\geq \sigma_n $为区间[0, 10]随机生成的$ A $的$ n $个奇异值. 令$ Z_u^{0}={\rm rand}(m, p)+{\rm rand}(m, p) $i, $ Z_v^{0}={\rm rand}(n, p)+{\rm rand}(n, p) $i. 针对算法Hybrid-Newton, 按下列方式生成初始迭代矩阵($ qf(\cdot) $表示对括号内矩阵进行QR分解(缩减版)所对应的Q因子)

(4.3) $ \begin{equation} \begin{array}{l} (a). (U^{(0)}, V^{(0)})=\left(qf(Z_u^{0}), qf(Z_v^{0})\right);\\ (b). (U^{(0)}, V^{(0)})=\left(qf(U_*+10*Z_u^{0}), qf(V_*+10*Z_v^{0})\right);\\ (c). (U^{(0)}, V^{(0)})=\left(qf(U_*+10^2*Z_u^{0}), qf(V_*+10^2*Z_v^{0})\right);\\ (d). (U^{(0)}, V^{(0)})=\left(qf(U_*+10^3*Z_u^{0}), qf(V_*+10^3*Z_v^{0})\right), \end{array} \end{equation} $

其中$ (U_*, V_*) $为算法OptStiefelGBB生成的问题(1.2)的解, 迭代初值取为$ \left(qf(Z_u^{0}), qf(Z_v^{0})\right) $, 终止标准取为$ \|{{{\rm{grad }}}} F(U^{(k)}, V^{(k)})\|\leq 10^{-6} $.

为说明算法的全局收敛性, 算法Hybrid-Newton取不同的迭代初值$ (U^{(0)}, V^{(0)}) $, 即(4.3)式中的(a)(b)(c)(d), 表 1给出了两算法的数值结果, 图 1给出了经算法OptStiefelGBB 生成满足条件的迭代初值的Newton法的黎曼梯度范数$ \|{{{\rm{grad }}}} f(U^{(k)}, V^{(k)})\| $的收敛曲线图. 从表 1和图 1中可知, 算法Hybrid-Newton 对初值选取依赖性不大, 这说明算法具有全局收敛性. 另外若由算法OptStiefelGBB得到合适的初值, 算法Newton 只需几步就能得到满足精度要求的迭代近似值.

为说明CR方法求解实对称线性方程组(3.29) 的有效性, 算法Hybrid-Newton取不同的迭代初值$ (U^{(0)}, V^{(0)}) $, 表 2给出了算法Newton每次循环中CR方法求解方程(3.29)的数值结果. 从表 2中可知, 应用CR方法求解由黎曼牛顿方程导出的实对称线性方程组是可行的.

例4.2   本例给出算法Hybrid-Newton与黎曼优化工具箱Manopt中的两类经典二阶算法, 即Exsiting-RTR和Exsiting-RBFGS, 以及黎曼梯度类算法OptStiefelGBB和RCG-Cayley的数值比较结果. 按下列方式生成问题(1.2) 中的矩阵$ A\in {\mathbb C}^{m\times n} $和$ \Theta\in {{\mathbb R}}^{p\times p} $:

$ A={\rm randn}(m, n)+{\rm randn}(m, n){\rm i}, {\quad} \Theta=\mbox{diag}(1, \cdots, p). $

为数值比较, 各算法的迭代初值均取为

$ (U^{(0)}, V^{(0)})=\left(qf(Z_u^{0}), qf(Z_v^{0})\right), $

其中$ Z_u^{0}={\rm randn}(m, p)+{\rm randn}(m, p){\rm i}, $$ Z_v^{0}={\rm randn}(n, p)+{\rm randn}(n, p) $i. 取不同的系数维数参数$ (m, n, p) $, 表 3给出了不同算法的数值比较结果, 图 2给出了不同算法黎曼梯度范数$ \|{{{\rm{grad }}}} f(U^{(k)}, V^{(k)})\| $随迭代步的收敛曲线图. 结合表 3的数据结果和图 2的收敛曲线图, 分析可知对任意迭代初值$ (U^{(0)} $, $ V^{(0)}) $, 算法OptStiefelGBB和算法RCG-Cayley均能收敛. 但当接近最优点时, 收敛速度会变得较慢.若以此为初值执行黎曼牛顿法, 因其具有局部二次收敛特性, 故一般只需额外几步就能得到满足精度要求的解. 黎曼牛顿法虽然每次迭代需要求解实对称线性方程组, 但因为迭代次数少, 总体迭代时间可以接受.

本文应用黎曼优化中的黎曼牛顿法和黎曼信赖域算法求解复矩阵截断奇异值分解问题(1.2). 对于黎曼牛顿法, 关键问题是如何有效求解对应的黎曼牛顿方程. 本文利用若干矩阵拉直算子和克罗内克积, 将其转化为容易求解的实对称线性方程组, 达到降维和简化计算的目的. 因为牛顿法对迭代初值比较敏感, 本文利用具有全局收敛性, 且迭代成本低, 存储量少的黎曼一阶经典算法OptStiefelGBB产生初值. 进而设计的混合算法同时具有全局收敛和局部二次收敛特性. 数值实验和数值比较验证, 对于中小规模的复矩阵截断奇异值分解问题, 本文设计的算法是切实有效的.

在数值实验中我们发现, 对于大规模的复矩阵截断奇异值问题, 黎曼牛顿法中求解实对称线性方程组需要花费大量时间. 故寻求更有效地方式求解黎曼牛顿方程是下一阶段的研究工作. 另外, 以迭代法求解黎曼牛顿方程导出的实对称线性方程组实质上是非精确地方式处理黎曼牛顿方程, 故寻求更有效方式非精确求解黎曼牛顿方程, 并细致分析对应非精确黎曼牛顿法, 如Guass-Newton, CG-Newton等的收敛性是下一步工作重点. 另外, 在目标函数中加入正则化项, 进而采用正则化牛顿法或邻近牛顿法求解, 也是下一步的研究工作.