嵌套子空间的增广拉格朗日法求解稀疏广义特征值问题

Nested Subspace Augmented Lagrangian Method for Solving Sparse Generalized Eigenvalue Problems

摘要: 稀疏广义特征值问题作为计算数学领域的重要研究课题，其核心在于确定满足广义特征方程$Ax=\lambda Bx$ 的特征对$\left(\lambda ,x\right)$ ，其中$A\in {ℝ}^{n\times n},B\in {ℝ}^{n\times n}$ 为给定的系数矩阵。值得注意的是，在工程结构分析、计算金融学、量子物理模拟及航空航天设计等实际应用场景中，所涉及的矩阵往往具有大规模、稀疏性以及非对称性等复杂特性，这对传统数值求解方法提出了严峻挑战。针对对称广义特征值问题的特殊数学结构，本文研究创新性地提出一种基于子空间构造与优化算法协同的数值计算方法。该方法通过建立具有良好逼近特性的子空间，结合增广拉格朗日半光滑牛顿法(ALMSsn)实现特征系统的高效求解，为解决大规模稀疏矩阵的广义特征值问题提供了新的技术途径。

Abstract: The sparse generalized eigenvalue problem, as a pivotal research subject in computational mathematics, focuses on determining eigenpairs $\left(\lambda ,x\right)$ that satisfy the generalized eigenvalue equation $Ax=\lambda Bx$ , where $A\in {ℝ}^{n\times n},B\in {ℝ}^{n\times n}$ are given coefficient matrices. It is noteworthy that in practical applications such as structural engineering analysis, computational finance, quantum physics simulations, and aerospace design, the involved matrices often exhibit complex characteristics including large-scale dimensions, sparsity, and asymmetry, posing significant challenges to conventional numerical solvers. Targeting the specific mathematical structure of symmetric generalized eigenvalue problems, this study innovatively proposes a numerical computational methodology that integrates subspace construction with optimization algorithms. By establishing subspaces with superior approximation properties and combining them with an augmented Lagrangian semi-smooth Newton method, the proposed approach enables efficient resolution of eigen-systems, thereby providing a novel technical pathway for addressing large-scale sparse generalized eigenvalue problems.

1. 引言

广义特征值问题作为科学计算领域的核心课题之一，在结构动力学分析、量子化学计算以及电磁场模拟等科学与工程领域具有重要应用价值。广义特征值问题目的是寻找一个矩阵的特征值和对应的特征向量，这些特征值和特征向量满足广义特征值方程：

$Ax=\lambda Bx$ ，$A\in {ℝ}^{n\times n}$ ，$B\in {ℝ}^{n\times n}$ ，

其中，A和B分别是已知的$n\times n$ 实数或复数矩阵，$\lambda$ 是特征值，$x$ 是n维列向量，是对应的特征向量。

当涉及对称广义特征值问题时(即A、B均为对称矩阵且B正定)，其典型应用场景如大型结构振动模态分析或电子能带计算中，系数矩阵往往呈现高维度、强稀疏性等特性。在此情形下，QR分解、QZ算法等传统数值方法往往面临计算效率瓶颈：一方面，算法复杂度随矩阵维度呈三次方增长；另一方面，矩阵稀疏性在正交变换过程中会遭到破坏，导致存储需求显著增加。

近三十年来，针对对称广义特征值问题的数值解法研究取得了显著进展。当前主流方法主要采用投影技术框架，其数学本质在于构建适当的低维子空间$S\subset {ℝ}^n$ ，将原始问题投影为$\left(S^{\text{T}}AS,S^{\text{T}}BS\right)$ 的降维广义特征值问题，进而通过求解子空间特征系统获得原问题的近似解。投影方法的计算效能与子空间构造策略密切相关，在众多子空间选择方案中，基于Krylov子空间的方法因其优异的收敛特性备受关注。其中经典的Lanczos算法通过将原问题投影至Krylov子空间：

${\mathcal{K}}_m\left(H,x\right):=\text{span}\left\{x,Hx,\cdots ,H^{m-1}x\right\}$ ，

其中，$H:=B^{-1}A$ ，x是某个初始向量，Lanczos算法构造了${\mathcal{K}}_m\left(H,x\right)$ 的基$Z$ ，然后计算投影问题$\left(ZAZ,ZBZ\right)$ 的特征对$\left(\theta ,u\right)$ ，特征值$\theta$ 称为Ritz值，并被视为特征方程$\left(A,B\right)$ 的近似特征值，其对应的近似特征向量$Zu$ ，称为Ritz向量。

Lanczos算法中，在构造$Z$ 时，并不需要显式形成$H$ 。它的收敛性得到了广泛的研究，该方法可以很好地逼近分离良好的极端特征对。然而，当极端特征值没有很好地分离时，Ritz值的收敛速度通常相当差。为了提高收敛速度，需要使用预条件来对原问题进行优化使其特征值更加聚拢[1]，可以将Lanczos算法应用在$B{\left(A-\mu B\right)}^{-1}BX=\widehat{\lambda }Bx$ 上，其中$\mu$ 是特定选取的，然后将$\widehat{\lambda }$ 转换为特征方程$\left(A,B\right)$ 的特征值，这种预条件策略有很快的收敛速度，但是这样做需要在每次迭代中求解与原问题大小相同的线性系统。由于问题规模的大小，若通过直接方法如LU分解法求解预条件的线性系统是不可行的。在这种情况下，可以使用迭代方法来计算近似解。然而，由于Lanczos算法对扰动非常敏感，尤其是在迭代的早期阶段，因此这种嵌套了线性系统的求解方法对每步线性系统的解有非常高的精度要求，这限制了此方法的实用性。

JDQZ方法[2]，也利用了预条件策略，它对预条件线性系统的近似解的精确度要求低很多。这种方法对所需的特征对进行初始近似，并依次求解校正方程来扩展搜索子空间。它有助于校正方程选择“右”向量进行扩展，从而以最佳方式细化特征对的近似值。在这里，预条件策略适用于逼近校正方程的解，但是这种方法仍需确定不精确解的程度，从而获得更精确的特征对。另外使用GMRES或BiCG等Krylov子空间方法的预处理策略已被广泛研究[3] [4]。

另一方面，可以使用梯度类方法，例如最速下降法来找到最小的特征值。这种方法不需要任何形式的求$B^{-1}$ 的运算，但是，当极端特征值没有很好地分离时，它们收敛性通常很差。一些已经研究了的梯度方法包括LOBPCG方法[5]、DACG方法[6]和逆迭代的变体，称为预条件逆迭代[7]。在这里，不是应用之前的预条件策略，而是使用一种不同类型的预处理策略来加速迭代。通常，选择适用于用矩阵$A$ 求解线性系统的前提条件并将其应用于原特征方程解的迭代中。尽管在许多情况下这些方法结果很好，但预处理策略在加速特征方程中的作用尚不清楚。

Golub和Ye提出并分析了一种不需要逆运算的Krylov子空间方法来计算对称广义特征方程$\left(A,B\right)$ 的极端特征值[8]。与梯度类型方法一样，该算法不需要随时计算$B^{-1}$ ，但与许多梯度方法不同，该算法利用了Krylov子空间提供更好的近似解，此外，还有用一种预处理策略来加速算法收敛到极值特征值的方法[8]。然而，由于向量迭代的性质，逆自由Krylov子空间算法在存在多个紧密聚类的特征值的情况下收敛较差。Quillen在Golub和Ye的工作基础上开发上一方法的块泛化，该方法具有相似的收敛特性[9]。

Alimisis和Saad在2023年提出了一种用于计算近似不变子空间的子空间迭代算法的变体，他们重新审视标准子空间迭代方法，并开发了利用梯度类型算法结合了Grassmann流形的新变体[10]。该算法有时优于标准的基于切比雪夫的子空间迭代，而且不需要估计最优参数。

Zhang和Shen在2018年将经典的特征值问题Lanczos方法扩展到了信赖域子问题[11]。由于子空间的维数对于良态的信赖域子问题影响通常较小，但对于病态问题影响则可能很大。会导致计算成本、内存需求和数值稳定性方面的数值困难。他们提出了一种高效的嵌套重启策略的Lanczos方法，解决了这一数值难题，并进行了收敛性分析。

本文提出了一种基于嵌套Lanczos算法与残差驱动子空间扩展的稀疏广义特征值问题求解方法。通过嵌套迭代动态生成多维子空间，在初始化阶段以初始向量构建基向量，后续通过Lanczos过程逐步扩展子空间，并结合优化算法更新当前解。每次迭代中，计算残差并对其进行正交化处理以保持子空间多样性，将正交化后的残差向量作为新增基向量，形成高维搜索空间以增强全局优化能力。在子空间内部通过增广拉格朗日半光滑牛顿算法求解投影后的问题，该算法显著提升了大规模稀疏问题的收敛效率与解的精度。

2. 预备知识

2.1. 广义特征值问题

设$A=\left(a_{ij}\right)\in {ℝ}^{n\times n}$ 是n阶实对称矩阵，$B=\left(b_{ij}\right)\in {ℝ}^{n\times n}$ 是n阶实对称正定矩阵，使得

$Ax=\lambda Bx$ ，

有非零解向量$x\in {ℝ}^n$ ，则称$\lambda$ 是矩阵A相对于矩阵B的特征值，且x是属于$\lambda$ 的特征向量，当$B\ne I$ 时，该问题为广义特征值问题。

定义1 如果一个函数$R\left(A,B,x\right)$ 满足$R\left(A,B,x\right)=\frac{x^{H}Ax}{x^{H}Bx}$ ，其中x为非零向量，A，B是n阶Hermitan矩阵，且B是正定矩阵，则称这样的函数为广义Rayleigh商。

2.2. Krylov子空间方法

Krylov子空间方法是特殊的Rayleigh-Ritz方法，其主要思想是在一个维数较小的子空间$\mathcal{K}\in {ℝ}^n$ 中寻找近似解。这类方法也被看作是一种投影方法，即寻找真解在某个子空间中的投影，可以是正交投影，也可以是斜投影。近似解的计算是通过迭代递归地找残差向量，第n个残差向量$r_n$ 可通过A的某个多项式与初始残差向量$r_0$ 相乘得到，即：$r_n=P\left(A\right)r_0$ 。

Krylov子空间方法有两个特征：1) 极小残差性，来保证收敛速度。2) 每一步迭代的计算量和储存量较小，保证计算的高效。

定义2 令$v\in {ℝ}^n,v\ne 0$ ，由$v,Av,\cdots ,A^{m-1}v$ 生成的子空间，称为由A和v生成的m维Krylov子空间，记作${\mathcal{K}}_m\left(A,v\right)$ ，即：${\mathcal{K}}_m\left(A,v\right):=\text{span}\left\{v,Av,\cdots ,A^{m-1}v\right\}$ 。

性质1 ${\mathcal{K}}_m\left(A,v\right)$ 是一个Krylov子空间，则：

1) Krylov子空间是嵌套的，即：${\mathcal{K}}_1\subseteq {\mathcal{K}}_2\subseteq \cdots \subseteq {\mathcal{K}}_m\subseteq \cdots$ ；

2) ${\mathcal{K}}_m\left(A,v\right)$ 的维数不超过m；

3) ${\mathcal{K}}_m\left(A,v\right)=\left\{x=P\left(A\right)v:P�次�不超�m-1的多�式\right\}$ 。

定义3 对于方程$Ax=b$ ，若${\mathcal{K}}_m\left(A,v\right)$ 是一个Krylov子空间，要求残量$r=b-Ax$ 满足m个正交性条件，即：$r=b-A\tilde{x}\perp ℒ$ 。其中，$\tilde{x}$ 是近似解，$ℒ$ 是另一个m维子空间，不同的$ℒ$ 对应不同的投影方法当$ℒ=\mathcal{K}$ 时，我们称为正交投影法，否则称为斜投影法。将$ℒ$ 称为约束空间，$\mathcal{K}$ 称为搜索空间。

下面给出投影方法的数学描述，给定初值$x^{\left(0\right)}\in {ℝ}^n$ ，$r_0=b-x^{\left(0\right)}$ 为初始残量，

$\begin{array}{l}\text{find}\tilde{x}=x^{\left(0\right)}+Vy\\ \text{s}.\text{t}.W^{\text{T}}AVy=W^{\text{T}}{r}_0.\end{array}$

其中，$V=\left[v_1,v_2,\cdots ,v_m\right]$ ，$W=\left[w_1,w_2,\cdots ,w_m\right]$ 分别是${\mathcal{K}}_m$ 和${ℒ}_m$ 的一组基。近似解$\tilde{x}$ 存在的唯一条件是矩阵$W^{\text{T}}AV$ 非奇异。

定理2 若$A,\mathcal{K},ℒ$ 满足下面条件之一：

(1) A正定，且$ℒ=\mathcal{K}$ ；

(2) A非奇异，且$ℒ=A\mathcal{K}$ ，

则$W^{\text{T}}AV$ 非奇异。

若$W^{\text{T}}AV$ 非奇异，可解得$y={\left(W^{\text{T}}AV\right)}^{-1}{W}^{\text{T}}{r}_0$ ，因此

$\tilde{x}=x^{\left(0\right)}+V{\left(W^{\text{T}}A\right)}^{-1}{W}^{\text{T}}{r}_0.$

下面我们介绍如何通过Arnoldi过程和Lanczos过程找到搜索空间$\mathcal{K}$ 的一组正交基。

2.3. Arnoldi过程

Arnoldi过程是求解非对称矩阵的一种正交投影方法，其主要思想就是通过对一个非对称矩阵A，产生Krylov子空间的一组标准正交基的方法。

算法1：基于Gram-Schmidt正交化的Arnoldi过程

1 给定非零向量r，计算$v_1=r/{\Vert r\Vert }_2$

2 当$j=1,2,\cdots ,m-1$ 执行

3 $w_j=A{v}_j{h}_{ij}$ ；

4 当$i=1,2,\cdots ,j$ 执行

5 $h_{ij}=\left(v_j,w_j\right)$ ；

6 $w_j=w_j-{\displaystyle {\sum }_{i=1}^j{h}_{ij}{v}_j}$ ；

7 $h_{j+1,j}={\Vert w_j\Vert }_2$ ；

8 若$h_{j+1,j}=0$ ，则

9 停止；

10 $v_{j+1}=w_j/h_{j+1,j}$ 。

令$H_{m+1,m}=\left[\begin{array}{cccc}{h}_{11}& h_{12}& \cdots & h_{1m}\\ h_{21}& h_{22}& \cdots & h_{2m}\\ 0& h_{32}& \cdots & h_{3m}\\ ⋮& ⋮& ⋮& ⋮\\ 0& \cdots & h_{m-1,m}& h_{m,m}\\ 0& \cdots & 0& h_{m+1,m}\end{array}\right]$

设$V_m=\left[v_1,v_2,\cdots ,v_m\right]$ ，$H_m=H_{m+1,m}\left(1:m,1:m\right)\in {ℝ}_{m\times m}$ ，即$H_m$ 是由$H_{m+1,m}$ 的前m行组成的Hessenberg矩阵。则：

$A{V}_m=V_{m+1}{H}_{m+1,m}=V_m{H}_m+h_{m+1,m}{V}_{m+1}+e_m^{\text{T}}$ ，

$V_m^{\text{T}}A{V}_m=H_m$ 。

定理3 如果到第$k\left(k<m\right)$ 步时有$h_{k+1,k}=0$ ，则算法将提前终止。此时$A{v}_k$ 必定可由$v_1,v_2,\cdots ,v_k$ 线性表出。则有：$A{V}_k=V_k{H}_k$ ，即${\mathcal{K}}_k$ 是A的一个不变子空间。

定理4 如果$h_{k+1,k}\ne 0$ ，$\left(k=1,2,\cdots ,m-1\right)$ ，则向量$v_1,v_2,\cdots ,v_m$ 构成${\mathcal{K}}_m$ 的一组标准正交基，其中${\mathcal{K}}_m\left(A,r\right)=\text{span}\left\{r,Ar,\cdots ,A_{m-1}r\right\}$ 。

2.4. Lanczos过程

对于方程$Ax=b$ ，若A是对称的，由Arnoldi过程的性质$H_m=V_m^{\text{T}}A{V}_m$ 可知$H_m$ 也是对称的，而$H_m$ 又是上Hessenberg矩阵，故$H_m$ 是一个对称三对角矩阵，记为$T_m$ ，记${\alpha }_j=h_{j,j}$ ，${\beta }_j=h_{j,j+1}$ ，则：

$T_m=\left[\begin{array}{cccc}{\alpha }_1& {\beta }_1& & \\ {\beta }_1& \ddots & \ddots & \\ & \ddots & \ddots & {\beta }_{m-1}\\ & & {\beta }_{m-1}& {\alpha }_m\end{array}\right]$

与Arnoldi过程类似，Lanczos过程具有如下性质：

$A{V}_m=V_m{T}_m+{\beta }_m{v}_{m+1}+e^{\text{T}}$ (1)

$V_m^{\text{T}}A{V}_m=T_m$ 。

考察关系式(1)两边的第j列，有：

${\beta }_j{v}_{j+1}=A{v}_j-{\alpha }_j{v}_j-{\beta }_{j-1}{v}_{j-1}$ ，

其中，$j=1,2,\cdots$ ，$v_0=0$ ，${\beta }_0=0$ ，基于这个三项递推公式，Anoldi过程简化为Lanczos过程：

算法2：Lanczos过程

1 给定非零向量r，令$v_0=0$ ，${\beta }_0=0$ ；

2 令$v_1=r/{\Vert r\Vert }_2$

3 当$j=1,2,\cdots ,m$ 执行

4 $w=A{v}_j-{\beta }_{j-1}{v}_{j-1}$ ；

5 ${\alpha }_j=\left(v_j,w_j\right)$ ；

6 $w_j=w_j-{\alpha }_j{v}_j$ ；

7 ${\beta }_j={\Vert w_j\Vert }_2$ ；

8 若${\beta }_j=0$ ，则

9 停止；

10 $v_{j+1}=w_j/{\beta }_j$ 。

2.5. 邻近算子

定义4 (邻近算子[12])对于一个凸函数$h\left(x\right)$ 定义它的邻近算子为：

${\text{prox}}_h\left(x\right)=\underset{u\in \text{dom}h}{\arg \min }\left\{h\left(u\right)+\frac{1}{2}{\Vert u-x\Vert }^2\right\}$ 。

邻近算子的目的是求解一个据x不算太远的点，并使函数值$h\left(x\right)$ 也相对较小，上述问题的解是唯一的，也就是对于凸函数，它的邻近算子是唯一的，下面定理给出此问题的解的存在唯一性。

定理6 (邻近算子与次梯度的关系)若$h\left(x\right)$ 是适当闭凸函数，则

$u={\text{prox}}_h\left(x\right)\iff x-u\in \partial h\left(u\right)$ 。

Proof. ($\Rightarrow$ )若$u={\text{prox}}_h\left(x\right)$ ，则由最优性条件得$0\in \partial h\left(u\right)+\left(u-x\right)$ ，因此$x-u\in \partial h\left(u\right)$ 。

($\Leftarrow$ )若$x-u\in \partial h\left(u\right)$ ，由次梯度的定义可得到：

$h\left(v\right)\ge h\left(u\right)+{\left(x-u\right)}^{\text{T}}\left(v-u\right)$ ，$\forall v\in \text{dom}h$ ，

两边同时加上$\frac{1}{2}{\Vert v-x\Vert }^2$ ，即有：

$\begin{array}{c}h\left(v\right)+\frac{1}{2}{\Vert v-x\Vert }^2\ge h\left(u\right)+{\left(x-u\right)}^{\text{T}}\left(v-u\right)+\frac{1}{2}{\Vert v-x\Vert }^2\\ \ge h\left(u\right)+\frac{1}{2}{\Vert u-x\Vert }^2,\forall v\in \text{dom}h,\end{array}$

因此我们得到$u={\text{prox}}_h\left(x\right)$ 。◻

邻近算子的计算可以看成是次梯度算法的隐式格式，对于非光滑情景，由于次梯度不唯一，显示迭代格式就不唯一，而隐式格式能得到唯一解，在步长选择上面，隐式格式也优于显示格式。

2.6. 近似点梯度法(PGA)

考虑如下复合优化问题：

$\min \psi \left(x\right)=f\left(x\right)+h\left(x\right)$ ，

其中，函数f为可微函数，定义域$\text{dom}f={ℝ}^n$ ，函数h是凸函数，可以是非光滑的。

近似点梯度法的思想是对目标函数$\psi \left(x\right)$ 的两部分，对一部分$f\left(x\right)$ 做梯度下降，对非光滑部分使用邻近算子，则近似点梯度法[12] [13]的迭代公式为：

$x^{k+1}={\text{prox}}_{t_{k}h}\left(x^k-t_k\nabla f\left(x^k\right)\right)$ ，

其中，$t_k$ 为每次迭代的步长，它可以是常数或者由线搜索得出。注意到，当$h\left(x\right)=0$ 时，迭代公式变为梯度下降法：

$x^{k+1}=x^k-t_k\nabla f\left(x^k\right)$ ；

当$h\left(x\right)=I_C\left(x\right)$ 时，$I_C\left(x\right)$ 是指示函数，则迭代公式变为投影梯度法：

$x^{k+1}=P_C\left(x^k-t_k\nabla f\left(x^k\right)\right)$ 。

下面给出近似点梯度法的算法：

算法3 邻近点梯度法

1 输入：函数$f\left(x\right)$ ，$h\left(x\right)$ ，初始点$x^0$ ，初始化$k\leftarrow 1$ ；

2 当未达到收敛准则 执行

3 $x^{k+1}={\text{prox}}_{t_{k}h}\left(x^k-t_k\nabla f\left(x^k\right)\right)$ ；

4 $k\leftarrow k+1$ 。

当我们根据邻近算子的定义，把迭代公式展开：

$\begin{array}{c}{x}^{k+1}=\underset{u}{\arg \min }\left\{h\left(u\right)+\frac{1}{2t_k}{\Vert u-x^k+t_k\nabla f\left(x^k\right)\Vert }^2\right\}\\ =\underset{u}{\arg \min }\left\{h\left(u\right)+f\left(x_k\right)+\nabla f{\left(x^k\right)}^{\text{T}}\left(u-x^k\right)+\frac{1}{2t_k}{\Vert u-x^k\Vert }^2\right\},\end{array}$

可以发现，近似点梯度法的实质就是将问题的光滑部分线性展开再加上二次项并保留非光滑部分，然后求极小作为每一步的估计。

此外根据定理6，还可以将近似点梯度法写成

$x^{k+1}=x^k-t_k\nabla f\left(x^k\right)-t_k{g}^k$ ，$g^k\in \partial h\left(x^{k+1}\right)$ ，

其本质是对光滑部分做显式梯度下降，非光滑部分做隐式梯度下降。在固定步长$t_k=t\in \left(0,\frac{1}{L}\right]$ 的情况下，迭代点$x^k$ 处的函数值$\psi \left(x^k\right)$ 以$\mathcal{O}\left(\frac{1}{k}\right)$ 的速率收敛到${\psi }^{\text{*}}$ 。

对于近似点梯度法，它的优点是：1) 简单易实现：近似点梯度法相对于其他复杂的优化算法而言，实现起来相对简单；2) 适用性广泛：近似点梯度法可以应用于各种类型的优化问题，包括凸优化和非凸优化问题；3) 收敛性保证：在某些条件下，近似点梯度法可以保证收敛到优化问题的局部最优解。

3. 稀疏广义特征值问题的嵌套子空间方法

3.1. 求解稀疏广义特征值问题的子空间

对于广义特征值问题

$\begin{array}{l}Ax=\lambda Bx\hfill \end{array}$ 。

其中，A是对称矩阵，B为正定对称矩阵，为了计算该问题最小特征值，将问题转化为带非光滑正则项和等式约束的优化问题：

$\begin{array}{ll}\underset{x}{\text{min}}\hfill & x^{\text{T}}Ax+\phi \left(x\right)\hfill \\ \text{s}.\text{t}.\hfill & x^{\text{T}}Bx=1.\hfill \end{array}$ (2)

我们的想法是通过Lanczos算法嵌套的迭代生成一个子空间$\text{span}\left(U\right)$ ，在子空间中求解问题(2)，具体而言，子空间的生成过程如下：

初始化子空间：在嵌套迭代的初始阶段，子空间的生成开始，将初始向量作为子空间的第一个基向量，初始向量可以是一个随机生成向量，通常服从标准正态分布或均匀分布。在后续的迭代过程中，根据Lanczos算法，逐步生成子空间的附加基向量，在每次迭代中，通过优化算法得到当前解$x_k$ 。计算当前解$x_k$ 与上一次迭代的解$x$ 之间的差值，即残差$d_s=x_k-x$ 。

对残差$d_s$ 进行正交化操作，以保持子空间的多样性。在正交化过程中，将$d_s$ 与子空间中的已有基向量进行正交化，确保$d_s$ 与子空间中的基向量线性无关，保证新基向量不重复已有方向，避免子空间退化。将正交化后的$d_s$ 加入到之前生成的子空间中，得到新的子空间。通过不断地迭代优化过程，每次生成的残差$d_s$ 都会加入到子空间中，从而形成一个多维子空间，这个多维子空间能够提供更多的搜索方向，帮助算法更好地搜索最优解。同时，通过正交化操作，保持了子空间的多样性，避免了子空间退化为低维空间的情况。这样，算法可以更有效地在高维空间中进行搜索，提高了收敛速度和全局搜索能力。

在嵌套Lanczos算法中的外部迭代过程中，每次迭代通过PGA算法求解如下子问题，以找到更优的更新方向$d$ 。

$\begin{array}{ll}\underset{x}{\text{min}}\hfill & g^{\text{T}}\left(x_k+d\right)+\frac{\sigma }{2}{\Vert d\Vert }^2+\phi \left(x_k+d\right)\hfill \\ \text{s}.\text{t}.\hfill & {\left(B{x}_k\right)}^{\text{T}}d=0.\hfill \end{array}$

其中$g=2Ax$ ，令$d=x-x_k$ ，则问题转换为

$\begin{array}{ll}\underset{x}{\text{min}}\hfill & g^{\text{T}}x+\frac{\sigma }{2}{\Vert x-x_k\Vert }^2+\phi \left(x\right)\hfill \\ \text{s}.\text{t}.\hfill & {\left(B{x}_k\right)}^{\text{T}}\left(x-x_k\right)=0.\hfill \end{array}$

给出上述问题的增广拉格朗日函数：

$\begin{array}{c}L\left(x,\lambda \right)=g^{\text{T}}x+\frac{\sigma }{2}{\Vert x-x_k\Vert }^2+\phi \left(x\right)-\lambda {\left(B{x}_k\right)}^{\text{T}}\left(x-x_k\right)\\ ={\left(g-\lambda B{x}_k\right)}^{\text{T}}x+\frac{\sigma }{2}{\Vert x-x_k\Vert }^2+\phi \left(x\right)+\lambda B{x}_k^{\text{T}}{x}_k\\ =\frac{\sigma }{2}{\Vert x-x_k+\frac{1}{\sigma }\left(g-\lambda B{x}_k\right)\Vert }^2+\phi \left(x\right)+\lambda B{x}_k^{\text{T}}{x}_k-\frac{1}{2\sigma }{\Vert g-\lambda B{x}_k\Vert }^2.\end{array}$

则$x_{k+1}$ ，$d_{k+1}$ 更新公式如下：

$x_{k+1}={\text{prox}}_{\phi /\sigma }\left(x_k-\frac{1}{\sigma }\left(g-\lambda B{x}_k\right)\right)$ ，

$d_{k+1}={\text{prox}}_{\phi /\sigma }\left(x_k-\frac{1}{\sigma }\left(g-\lambda B{x}_k\right)\right)-x_k$ 。

于是有

${\left({\text{prox}}_{\phi /\sigma }\left(x_k-\frac{1}{\sigma }\left(g-\lambda B{x}_k\right)\right)-x_k\right)}^{\text{T}}B{x}_k=0$ 。

通过以上子空间生成过程，嵌套迭代能够生成一个逐步优化的子空间，用于求解稀疏最小广义特征值问题。子空间的生成是嵌套迭代的核心部分，为算法提供了一种高效的方式来近似求解问题的最优解。

3.2. 算法投影分析

在本节中我们分析投影到子空间后的小问题：

$\begin{array}{ll}\underset{x\in x_k+\text{span}\left(U\right)}{\text{argmin}}\hfill & x^{\text{T}}Ax+\phi \left(x\right)\hfill \\ \text{s}.\text{t}.\hfill & x^{\text{T}}Bx=1.\hfill \end{array}$

其中，$\text{span}\left(U\right)$ 是求解问题所在的子空间，$Q_k$ 是$\text{span}\left(U\right)$ 中的一组正交基，则下一次迭代值可以表示为

$\tilde{x}=x_k+Q_{k}h$ 。

进而优化函数可表示为：

$\begin{array}{c}F\left(\tilde{x}\right)={\left(x_k+Q_{k}h\right)}^{\text{T}}A\left(x_k+Q_{k}h\right)+\phi \left(x_k+Q_{k}h\right)\\ =x_k^{\text{T}}A{x}_k+h^{\text{T}}{Q}_k^{\text{T}}A{Q}_{k}h+2h^{\text{T}}{Q}_k^{\text{T}}A{x}_k+\phi \left(x_k+Q_{k}h\right)\end{array}$。

约束条件${\tilde{x}}^{\text{T}}B\tilde{x}=1$可以转化为：

$\begin{array}{l}1={\left(x_k+Q_{k}h\right)}^{\text{T}}B\left(x_k+Q_{k}h\right),\\ 1=x_k^{\text{T}}B{x}_k+h^{\text{T}}{Q}_k^{\text{T}}B{Q}_{k}h+2h^{\text{T}}{Q}_k^{\text{T}}B{x}_k,\\ {\Vert Q_k^{\text{T}}B{x}_k\Vert }_2^2={\Vert h\Vert }_2^2+2h^{\text{T}}{Q}_k^{\text{T}}B{x}_k+{\Vert Q_k^{\text{T}}B{x}_k\Vert }_2^2,\\ {\Vert Q_k^{\text{T}}B{x}_k\Vert }_2^2={\Vert h+Q_k^{\text{T}}B{x}_k\Vert }_2^2.\end{array}$

令$z=h+Q_k^{\text{T}}B{x}_k$ ，则上式约束条件可以转化为：

${\Vert z\Vert }_2^2={\Vert Q_k^{\text{T}}B{x}_k\Vert }_2^2\to \Delta$ 。

由于$h=z-Q_k^{\text{T}}B{x}_k$ ，优化函数可以进一步转变：

$\begin{array}{l}{h}^{\text{T}}{Q}_k^{\text{T}}A{Q}_{K}h+2h^{\text{T}}{Q}_k^{\text{T}}A{x}_k+\phi \left(x_k+Q_{k}h\right)\\ ={\left(z-Q_k^{\text{T}}B{x}_k\right)}^{\text{T}}{Q}_k^{\text{T}}A{Q}_K\left(z-Q_k^{\text{T}}B{x}_k\right)+2{\left(z-Q_k^{\text{T}}B{x}_k\right)}^{\text{T}}{Q}_k^{\text{T}}A{x}_k+\phi \left(x_k+Q_k\left(z-Q_k^{\text{T}}B{x}_k\right)\right)\\ =z^{\text{T}}{Q}_k^{\text{T}}A{Q}_{k}z+2z^{\text{T}}\left(Q_k^{\text{T}}A{x}_k-Q_k^{\text{T}}A{Q}_k^{\text{T}}{Q}_{k}B{x}_k\right)+{\left(Q_k^{\text{T}}B{x}_k\right)}^{\text{T}}{Q}_k^{\text{T}}A{Q}_k{Q}_k^{\text{T}}B{x}_k+\phi \left(x_k+Q_k\left(z-Q_k^{\text{T}}B{x}_k\right)\right)\\ =z^{\text{T}}{Q}_k^{\text{T}}A{Q}_{k}z+2z^{\text{T}}{Q}_k^{\text{T}}A\left(x_k-Q_k^{\text{T}}{Q}_{k}B{x}_k\right)+{\left(Q_k^{\text{T}}B{x}_k\right)}^{\text{T}}{Q}_k^{\text{T}}A{Q}_k{Q}_k^{\text{T}}B{x}_k+\phi \left(x_k+Q_{k}z-Q_k{Q}_k^{\text{T}}B{x}_k\right)\\ =z^{\text{T}}{Q}_k^{\text{T}}A{Q}_{k}z+2z^{\text{T}}{Q}_k^{\text{T}}A\left(I-Q_k^{\text{T}}{Q}_{k}B\right)x_k+{\left(Q_k^{\text{T}}B{x}_k\right)}^{\text{T}}{Q}_k^{\text{T}}A{Q}_k{Q}_k^{\text{T}}B{x}_k+\phi \left(Q_{k}z+\left(I-Q_k{Q}_k^{\text{T}}B\right)x_k\right)\end{array}$ 。

令$\left(I-Q_k^{\text{T}}{Q}_{k}B\right)x_k=c_k$ ，$Q_k^{\text{T}}A{Q}_k=V$ ，

通过上述投影变换，将原问题转换成如下优化问题：

$\begin{array}{ll}\underset{z}{\text{min}}\hfill & z^{\text{T}}Vz+2z^{\text{T}}{Q}_k^{\text{T}}A{c}_k+\phi \left(Q_{k}z+c_k\right)\hfill \\ \text{s}.\text{t}.\hfill & {\Vert z\Vert }_2^2=\Delta .\hfill \end{array}$ (3)

3.3. 增广拉格朗日法求解投影问题

对于问题(3)，令$y=Q_{k}z+c_k$ ，得到：

$\begin{array}{ll}\underset{z}{\text{min}}\hfill & z^{\text{T}}Vz+2z^{\text{T}}{Q}_k^{\text{T}}A{c}_k+\phi \left(Q_{k}z+c_k\right)\hfill \\ \text{s}\text{.t}\hfill & {\Vert z\Vert }_2^2=\Delta ,\hfill \\ \hfill & y=Q_{k}z+c_k.\hfill \end{array}$

其中，$\Delta ={\Vert Q_k^{\text{T}}B{x}_k\Vert }_2^2$ ，这是一个带约束的流形优化问题，可以采用流形上的增广拉格朗日法，将问题转换成如下形式：

$\begin{array}{ll}\underset{z\in ℳ}{\text{min}}\hfill & z^{\text{T}}Vz+2z^{\text{T}}{Q}_k^{\text{T}}A{c}_k+\phi \left(Q_{k}z+c_k\right)\hfill \\ \text{s}.\text{t}.\hfill & y=Q_{k}z+c_k.\hfill \end{array}$

其增广拉格朗日函数为：

$\begin{array}{c}{L}_{\sigma }\left(z,y,\tau \right)=z^{\text{T}}Vz+2z^{\text{T}}{Q}_k^{\text{T}}A{c}_k+\phi \left(y\right)+{\tau }^{\text{T}}\left(Q_{k}z+c_k-y\right)+\frac{\sigma }{2}{\Vert Q_{k}z+c_k-y\Vert }_2^2\\ =z^{\text{T}}Vz+2z^{\text{T}}{Q}_k^{\text{T}}A{c}_k+\phi \left(y\right)+\frac{\sigma }{2}{\Vert Q_{k}z+c_k-y+\frac{\tau }{\sigma }\Vert }_2^2-\frac{{\Vert \tau \Vert }_2^2}{2\sigma }.\end{array}$

其中，$\tau$ 是对应线性约束$y=Q_{k}z+c_k$ 的拉格朗日乘子，我们需要求解如下问题：

$\begin{array}{ll}\underset{z,y,\tau }{\text{min}}\hfill & L_{\sigma }\left(z,y,\tau \right)\hfill \\ \text{s}.\text{t}.\hfill & z\in ℳ.\hfill \end{array}$ (4)

下面给出增广拉格朗日算法解决上述问题：

算法4：增广拉格朗日法求解投影问题

1 输入：${\tau }_0\in {ℝ}^m$ ；

3 当$j=0,1,\cdots$ 执行

4 $\left(z^{j+1},y^{j+1}\right)=\underset{z\in ℳ,y\in \text{span}\left(U\right)}{\arg \min }{L}_{\sigma }\left(z,y,{\tau }^j\right)$ ；

5 ${\tau }^{j+1}={\tau }^j-\tau \left(Q_k{z}^{j+1}-c_k{y}^{j+1}\right)$ ；

8 若达到收敛准则，则

9 停止；

10 增加$\sigma$

下面分析增广拉格朗日子问题如何求解。

3.4. 半光滑牛顿法求解子问题(4)

在这一小节中，我们用半光滑牛顿法解决子问题(4)，对$z,y$ 同时最小化等价于：

$\underset{z\in ℳ}{\min }{z}^{\text{T}}Vz+2z^{\text{T}}{Q}_k^{\text{T}}A{c}_k+{\phi }^{\sigma }\left(Q_{k}z+c_k+\frac{\tau }{\sigma }\right)$ 。

其中，${\phi }^{\sigma }$ 是Moreau-Yosida正则化，$w=Q_{k}z+c_k+\frac{\tau }{\sigma }$ ，此函数是连续可微的：

${\phi }^{\sigma }\left(w\right):=\underset{y\in \text{span}\left(U\right)}{\min }\phi \left(y\right)+\frac{\sigma }{2}{\Vert w-y\Vert }_2^2.$ (5)

公式(5)与计算对应函数的近似点映射${\text{prox}}_{f/\sigma }\left(x\right):=\underset{y}{\arg \min }f\left(y\right)+\frac{\sigma }{2}{\Vert x-y\Vert }_2^2$ 是等价的，故：

$y^{j+1}={\text{prox}}_{\phi /\sigma }\left(Q_k{z}^{j+1}+c_k+\frac{{\tau }^j}{\sigma }\right).$ (6)

令

$\Phi \left(z\right):=z^{\text{T}}Vz+2z^{\text{T}}{Q}_k^{\text{T}}A{c}_k+{\phi }^{\sigma }\left(Q_{k}z+c_k+\frac{\tau }{\sigma }\right),$ (7)

其中$z\in ℳ$ ，则

$\begin{array}{c}\nabla \Phi \left(z\right)=2Vz+2Q_k^{\text{T}}A{c}_k+{\nabla }_z{\phi }^{\sigma }\left(Q_{k}z+c_k+\frac{\tau }{\sigma }\right)\\ =2Vz+2Q_k^{\text{T}}A{c}_k+\sigma Q_k^{\text{T}}\left[Q_{k}z+c_k+\frac{\tau }{\sigma }-{\text{prox}}_{\phi /\sigma }\left(Q_{k}z+c_k+\frac{{\tau }^j}{\sigma }\right)\right]\\ =2Vz+2Q_k^{\text{T}}A{c}_k+\sigma Q_k^{\text{T}}{Q}_{k}z+\sigma Q_k^{\text{T}}{c}_k+\tau Q_k^{\text{T}}-\sigma \left[{\text{prox}}_{\phi /\sigma }\left(Q_{k}z+c_k+\frac{{\tau }^j}{\sigma }\right)\right].\end{array}$

因此我们得到如下命题来描述问题(4)的一阶必要条件。

命题1 给定${\tau }^j$ ，对任意$y^{j+1}$ 满足公式(6)和$\Vert \nabla \Phi \left(z^j\right)\Vert =0$ ，则

$\left(z^{j+1},y^{j+1}\right)=\underset{z\in ℳ,y\in \text{span}\left(U\right)}{\arg \min }{L}_{\sigma }\left(z,y,\tau \right).$

除此之外对于适当的足够大的$\sigma$ ，存在很小的$\epsilon >0$ ，使得

$\Vert Q_k{z}^{j+1}+c_k-y^{j+1}\Vert =Q_{k}z+c_k-{\text{prox}}_{\phi /\sigma }\left(Q_k{z}^j{v}_{0}A{v}_0+c_k+\frac{{\tau }^j}{\sigma }\right)\le \epsilon .$

命题2 令${\partial }^2\Phi \left(z\right)$ 为$\Phi \left(z\right)$ 的Clark广义Jacbbi，则

$\begin{array}{c}{\partial }^2\Phi \left(z\right)=\partial \nabla \Phi \left(z\right)\\ =2V+\sigma Q_k^{\text{T}}{Q}_k-\sigma Q_k^{\text{T}}\partial {\text{prox}}_{\phi /\sigma }\left(Q_{k}z+c_k+\frac{\tau }{\sigma }\right)Q_k\\ =2V+\sigma Q_k^{\text{T}}{Q}_k-\sigma Q_k^{\text{T}}W{Q}_k.\end{array}$

其中，$W\in \partial {\text{prox}}_{\frac{\phi }{\sigma }}\left(Q_{k}z+c_k+\frac{\tau }{\sigma }\right)$ 。

下面给出半光滑牛顿法具体算法：

算法5：半光滑牛顿法求解子问题(4)

1 输入：初始点$z_0=z_j\in ℳ$ ，$t\in \left(0,1\right)$ ，$\delta \in \left(0,1/2\right)$ ，$k=0$ ，${\tau }^j$

续表

2 当达到收敛准则 执行

3 选择$H^k\in \partial \nabla \Phi \left(z^k\right)$ ；

4 解d通过求解广义牛顿方程：$H_{k}d=-\partial \Phi \left(z_k\right)$ ；

5 令$\alpha =1$ ；

6 若$\Phi \left(z^k+\alpha d^k\right)>\Phi \left(z^k\right)+\alpha \Phi {\left(z^k\right)}^{\text{T}}{d}^k$ 执行

7 $\alpha =t\alpha$ ；

8 更新${\alpha }^k=\alpha$ 且$z^{k+1}=z^k+{\alpha }^k{d}^k$ ，$k=k+1$

半光滑牛顿法子问题，可以直接调用黎曼流形上的信赖域方法求解，根据$z=h+Q_k^{\text{T}}B{x}_k$ 计算可以得到$x_k$ ，进而能求出对应的特征值。下面给出半光滑牛顿法的收敛性分析。

定义7 如果F在x点所B-微分的元素$J\in {\partial }_{B}F\left(x\right)$ 都是非奇异的，那么称F在x点是BD-正则的。

引理6 ([14])若$\Phi$ 是半光滑的和BD-正则的，则存在常数$c>0$ ，$\kappa >0$ 和一个小邻域$N\left(z^{*},{\epsilon }_0\right)$ 使得对于任意的$U\in N\left(z^{*},{\epsilon }_0\right)$ 和$J\in {\partial }_B\Phi \left(U\right)$ ，下面的结论成立：

1) $z^{*}$ 是一个孤立解；

2) J是非奇异的并且${\Vert J^{-1}\Vert }_F\le c$ ；

3) 局部误差界条件对于$\Phi \left(z\right)$ 在邻域$N\left(z^{*},{\epsilon }_0\right)$ 上成立，也就是说${\Vert U-z^{*}\Vert }_2\le \kappa \Vert \Phi \left(U\right)\Vert$ 。

定理9 若$\nabla \Phi \left(z\right)$ 是半光滑的且BD-正则的，并且z是优化问题(4)的最优解。那么z的迭代是良定义的，且存在一个小邻域$N\left(z^{*},\epsilon \right)$ ，使得对于任意的k有$X{z}^k\in N\left(z^{*},\epsilon \right)$ ，迭代是超线性收敛的。如果$\nabla \Phi \left(x\right)$ 是强半光滑的，迭代是二次收敛的。

Proof. 根据引理6，迭代是良定义的。我们可以简单地推出

$\begin{array}{c}{\Vert z^{k+1}-z^{*}\Vert }_2={\Vert z^k-J_k^{-1}\nabla \Phi \left(z^k\right)-z^{*}\Vert }_2\\ \le {\Vert J_k^{-1}\Vert }_F\cdot {\Vert \nabla \Phi \left(z^k\right)-\nabla \Phi \left(z^{*}\right)-J_k\left(z^k-z^{*}\right)\Vert }_2\\ =o\left({\Vert z^k-z^{*}\Vert }_2\right).\end{array}$

其中最后一个等式来源于$\nabla \Phi \left(z\right)$ 的半光滑性。如果$\nabla \Phi \left(z\right)$ 是强半光滑的，则其可以改为

$O\left({\Vert z^k-z^{*}\Vert }_2^2\right)$ 。

4. 数值实验

4.1. 实验

在本节中，我们评估了本文算法在解决对称稀疏广义特征值问题时的性能。在Matlab (R2017a)中编码，并在运行Windows 10 (64位)系统、配备Intel (R) Core (TM) i5-10210U CPU @1.60 GHz和8 GB内存的PC上进行评估。

对于函数$\phi \left(x\right)$ ，我们选用$l_1$ 范数进行测试，并将算法与邻近梯度法进行对比实验。

$l_1$ 范数

$l_1={\Vert X\Vert }_1$ ，

$l_1$ 函数的邻近算子为

$\begin{array}{c}{\text{prox}}_{{\phi }_L/\sigma }\left(x\right)=\underset{y}{\arg \min }{\phi }_L\left(y\right)+\frac{\sigma }{2}{\Vert x-y\Vert }_2^2\\ =\text{sign}\left(x\right)\text{max}\left\{\left|x\right|-1/\sigma ,0\right\},\end{array}$

$l_1$ 函数的广义雅可比函数为

${\left[J_L\right]}_{ii}=\{\begin{array}{ll}0,\hfill & \left|x_i\right|<1/\sigma \hfill \\ \left[0,1\right],\hfill & \left|x_i\right|=1/\sigma \hfill \\ 1,\hfill & \left|x_i\right|>1/\sigma \hfill \end{array}$

4.2. 实验结果

参数设置：Krylov子空间${\mathcal{K}}_m\left(A,r\right)$ 维度m = (10, 20)，Krylov子空间${\mathcal{K}}_m\left(A,Bx\right)$ 维度m = 1，外层最大子空间维度50，收敛容忍度1e-5，稀疏度1e-3。

Krylov子空间${\mathcal{K}}_m\left(A,r\right)$ 维度m = 10的nested ALM实验结果见表1，Krylov子空间${\mathcal{K}}_m\left(A,r\right)$ 维度m = 20的nested ALM实验结果见表2，PGA实验结果见表3。

Table 1. Nested ALM experimental results (m = 10)

表1. 嵌套子空间的增广拉格朗日法实验结果(m = 10)

n	CPU时间	残差	外层迭代步数	目标函数值
1000	3.1051	1.4287	295	−0.0879
1500	7.8445	1.4308	457	−0.0723
2000	26.5017	1.3940	906	−0.0627
2500	28.7841	1.4191	766	−0.0564
3000	67.3727	1.3800	1228	−0.0514

Table 2. Nested ALM experimental results (m = 20)

表2. 嵌套子空间的增广拉格朗日法实验结果(m = 20)

n	CPU时间	残差	外层迭代步数	目标函数值
1000	2.6639	1.3963	247	−0.0879
1500	11.3557	1.3356	603	−0.0723
2000	25.6663	1.3885	855	−0.0627
2500	33.2940	1.4138	872	−0.0564
3000	63.5617	1.4169	1126	−0.0514

Table 3. PGA experimental results

表3. 邻近梯度法实验结果

n	CPU时间	残差	迭代步数	目标函数值
1000	7.7420	1.4298	1369	−0.0879
1500	8.3456	1.4003	800	−0.0723
2000	37.8297	1.3911	1934	−0.0627
2500	39.7651	1.4258	1529	−0.0564
3000	79.2420	1.4543	2044	−0.0514

由实验结果可以看出，本文提出嵌套子空间的增广拉格朗日法，当子空间维数${\mathcal{K}}_m\left(A,r\right)$ 维度为20时的时间优于维度为10的时候，但残差略逊于子空间维度取10时，在实际应用中可以根据具体问题适当调整参数，可以根据残差下降速率自适应调整维度。例如，若残差下降缓慢，则增加维度以提升搜索能力。无论子空间维度取何值本文算法在CPU时间上都明显优于PGA，说明本算法对求解对称稀疏广义特征值问题具有较高的求解效率和稳定性。

5. 结论

本文实现了一种解决对称稀疏广义特征值问题的算法，该算法结合了嵌套结构和增广拉格朗日算法的思想，用于该问题的迭代求解。算法的核心是使用PGA算法来确定搜索方向并将其加入到子空间中，并使用子空间的信息进行进一步优化，整个算法动态地生成子空间，嵌套结构的引入可以进一步提高算法的性能，使其在复杂问题中寻找较好的解。并且实验结果表明，本算法对求解对称稀疏广义特征值问题具有较高的求解效率和稳定性，为了取得更好的性能，在实际应用中可以根据具体问题和数据适当调整算法参数。