大规模矩阵对角化方法：Lanczos

1. A的本征值即其Rayleigh商的极值

Rayleigh商定义：

对于nxn实对称阵A，任意非零n维矢量\(\vec{x}\neq 0\)，定义Rayleigh商为
\(r_A( \vec{x} ) = \frac{ \vec{x}^\top A \vec{x} }{ \vec{x}^\top \vec{x} }.\)
显然，Rayleigh商对于\(\vec{x}\)的模并不依赖。

Rayleigh商的范围

假设A的本征值是\(\lambda_i, i=1,...,n\),从小到大排列：\(\lambda_1 \leq \lambda_2 ... \leq \lambda_n\)。
本征矢是\(\vec{v}_i, i=1,...,n\)，有\(A\vec{v}_i = \lambda_i \vec{v}_i\)。
将非零矢量\(\vec{x}\)用\(\vec{v}_i\)表达，有\(\vec{x} = \sum_i c_i \vec{v}_i\)，则Rayleigh商为
\(r_A ( \vec{x} ) = \frac{ \sum_i c_i^2 \lambda_i^2 }{ \sum_i c_i^2 } \in [ \lambda_1， \lambda_n ]\)。

Rayleigh商的极值

计算Rayleigh商对\(\vec{x}\)的梯度
\(\nabla r_A(\vec{x}) = \frac{2}{\vec{x}^\top \vec{x}}( A \vec{x} - r_A(\vec{x})\vec{x} ).\)
取\(\nabla r_A(\vec{x}) = 0 \Leftrightarrow A \vec{x} = r_A(\vec{x})\vec{x}\)。
所以，A的Rayleigh商\(r_A(\vec{x})\)在其本征矢\(\vec{v}\)处取得极值，极值即A的本征值\(\lambda_i\)。

2. 子空间中的Rayleigh商

若有\(k\)个（\(k\leq n\)）正交归一的n维矢量\(\{ \vec{q}_1, \vec{q}_1, \cdots \vec{q}_k \}\)，构成子空间，记nxk的矩阵\(Q_k \equiv \begin{bmatrix} \vec{q}_1, \vec{q}_1, \cdots, \vec{q}_{k} \end{bmatrix}\)，
这个子空间中任意非零矢量为\(\vec{x} = \sum^k_{i=1} y_i \vec{q}_i = Q_k \vec{y}\)，\(\vec{y}\)是一个k维非零矢量，则A在这个子空间中的Rayleigh商为
\(r_A(Q_k \vec{y}) = \frac{ (Q_k \vec{y})^\top A (Q_k \vec{y}) }{ (Q_k \vec{y})^\top (Q_k \vec{y}) } = \frac{ \vec{y}^\top (Q^\top_k A Q_k) \vec{y} }{ \vec{y}^\top \vec{y} } = r_{Q^\top_k A Q_k}(\vec{y})\)。
转化为\(Q^\top_k A Q_k\)的Rayleigh商。

3. Krylov子空间：快速逼近极端本征值

根据上面的分析，要想构造很小的子空间，使得该子空间内存在最大/最小本征值对应的本征矢\(\vec{v}_1, \vec{v}_n\)的近似矢量，则需要构造\(Q_k\)，使得\(r_{Q^\top_k A Q_k}\)的最小值尽量小，最大值尽量大。
我们可以从选定的\(\vec{q}_1\)开始，逐渐增加\(\vec{q}_2, \vec{q}_3, ...\)，逐渐逼近目标。
一个有用的性质是：\(\nabla r_A(\vec{q}_1) = \frac{2}{\vec{q}^\top_1 \vec{q}_1}( A \vec{q}_1 - r_A(\vec{q}_1)\vec{q}_1 )\)，这个梯度是\(A\vec{q}_1, \vec{q}_1\)的线性叠加。
于是，美妙的想法产生了：何不如此构造子空间呢：\(\kappa(A, \vec{q}_1, k) = \{ \vec{q}_1, A\vec{q}_1, A^2\vec{q}_1, \cdots, A^{k-1}\vec{q}_1 \}\)，这个子空间叫做Krylov子空间。这样构造的话，对于\(Q_k\)，A的Rayleigh商在\(Q_k\)构造的子空间中任意一处的梯度矢量都在\(Q_{k+1}\)的子空间中。
所以，每次扩大子空间时：\(Q_1, Q_2, \cdots\)，Rayleigh商的最大值都会尽可能地增大，Rayleigh商的最小值都会尽可能地减小。
在这个意义上，Lanczos方法可以看做是Rayleigh商\(r_A\)的一个梯度下降/上升过程。
这和共轭梯度法很相似，共轭梯度法相当于：每次在\(Q_k\)内找一个下降方向，并保证这些下降方向线性无关，所以最多只会有\(n\)个下降方向，因此\(n\)步就能到达极值。当然，在大规模对角化问题中，如果只需要极端本征值，\(n\)步是奢侈的，Lanczos中迭代次数远小于\(n\)。

4. Lanczos方法

如上所述的策略很美妙，但还有一个小小的问题：我们如何判断子空间已经足够大了呢？换一句话说，我们如何判断子空间中的极端本征值已经足够大，或足够小？
在\(k\)逐渐增大时，我们需要在\(\kappa(A, \vec{q}_1, k)\)中不断地对角化，得到各个子空间中的极端本征值，看它们是否收敛。如果随着\(k\)的增大，极端本征值的变化极其微小，我们则经验地认为已经收敛。
要在\(\kappa(A, \vec{q}_1, k)\)中对角化A，就需要正交归一化\(\vec{q}_1, A\vec{q}_1, \cdots\)。
Lanczos方法即如下步骤，构造\(\vec{q}_1, \vec{q}_2, \cdots\)，既保证每次增加矢量，得到的子空间都是Krylov子空间，又保证\(\vec{q}_1, \cdots, \vec{q}_2, \cdots\)是正交归一矢量。

正交化方案：

给定单位矢量\(\vec{q}_1\)，定义\(\vec{r}_2 = A \vec{q}_1 - ( \vec{q}_1^\top A \vec{q}_1 ) \vec{q}_1, \vec{q}_2 = \vec{r}_2 / | \vec{r}_2 |\)。即从\(A\vec{q}_1\)中剔除\(\vec{q}_1\)的分量，然后归一化。
对于\(k\geq 2\)：\(\vec{r}_{k+1} = A \vec{q}_k - (\vec{q}^\top_k A \vec{q}_k) \vec{q}_k - (\vec{q}^\top_{k-1} A \vec{q}_k) \vec{q}_{k-1}, \vec{q}_{k+1} = \vec{r}_{k+1} / | \vec{r}_{k+1} |\)。即从\(A\vec{q}_k\)中剔除\(\vec{q}_{k-1}, \vec{q}_k\)的分量，然后归一化。

正交化吗？归纳法证明

显然，如上构造的子空间\(\{\vec{q}_1, \vec{q}_2\} = \{\vec{q}_1, A\vec{q}_1 \}, \{\vec{q}_1, \vec{q}_2, \vec{q}_3 \} = \{ \vec{q}_1, A\vec{q}_1, A^2\vec{q}_1 \}\)。并且\(\vec{q}_1, \vec{q}_2, \vec{q}_3\)正交。

以此为起点，开始归纳法证明：如果\(i=1,2,\cdots,k\)，都有\(\{\vec{q}_1, \cdots, \vec{q}_i\} = \{ \vec{q}_1, A\vec{q}_1, \cdots, A^{i-1}\vec{q}_1 \}\)并且\(\vec{q}_1, \cdots, \vec{q}_i\)两两垂直,那么，如上构造\(\vec{q}_{k+1}\)以后，也会有\(\{\vec{q}_1, \cdots, \vec{q}_i\} = \{ \vec{q}_1, A\vec{q}_1, \cdots, A^{i-1}\vec{q}_{k+1} \}\)并且\(\vec{q}_1, \cdots, \vec{q}_{k+1}\)两两垂直。

\(\{\vec{q}_1, \cdots, \vec{q}_k, \vec{q}_{k+1}\} = \{ \vec{q}_1, A\vec{q}_1, \cdots, A^{k}\vec{q}_1 \}\)

由于\(A\vec{q}_1, A\vec{q}_2, \cdots, A\vec{q}_{k-2}\)都是子空间\(\{ \vec{q}_1, \vec{q}_2, \cdots, \vec{q}_{k-1} \}\)中的矢量，所以都垂直于\(\vec{q}_k\)，即
\(\vec{q}_k^\top A \vec{q}_{k-2} = \vec{q}^\top_k A \vec{q}_{k-3} = \cdots \vec{q}^\top_k A \vec{q}_1 = 0\)。
构造\(\vec{r}_{k+1} = A \vec{q}_k - (\vec{q}^\top_k A \vec{q}_k) \vec{q}_k - (\vec{q}^\top_{k-1} A \vec{q}_k) \vec{q}_{k-1}, \vec{q}_{k+1} = \vec{r}_{k+1} / | \vec{r}_{k+1} |\)，
则显然有\(\{\vec{q}_1, \cdots, \vec{q}_k, \vec{q}_{k+1}\} = \{ \vec{q}_1, A\vec{q}_1, \cdots, A^{k}\vec{q}_1 \}\)，剩下的任务是证明\(\vec{q}_{k+1}\)与\(\vec{q}_1, \cdots, \vec{q}_k\)垂直。

\(\vec{q}_{k+1}\)与\(\vec{q}_k, \vec{q}_{k-1}\)都正交

由于\(\vec{q}_{k+1}\)是\(A\vec{q}_k\)剔除\(\vec{q}_k, \vec{q}_{k-1}\)的成分构造的，自然与\(\vec{q}_k, \vec{q}_{k-1}\)正交。

\(\vec{q}_{k+1}\)与\(\vec{q}_1, \vec{q}_2, \cdots \vec{q}_{k-2}\)都正交

因为\(\vec{q}_k\)与子空间\(\{ \vec{q}_1, \cdots \vec{q}_{k-1} \}\)正交，所以\(\vec{q}_k\)与\(A\vec{q}_1, A\vec{q}_2, \cdots, A\vec{q}_{k-2}\)都正交，即
\(\vec{q}^\top_k A \vec{q}_1 = \vec{q}^\top_k A \vec{q}_2 = \cdots = \vec{q}^\top_k A \vec{q}_{k-2} = 0\)。
利用A的对称性，上式即
\(\vec{q}^\top_1 A \vec{q}_k = \vec{q}^\top_2 A \vec{q}_k = \cdots = \vec{q}^\top_{k-2} A \vec{q}_{k} = 0\)。

因此，\(\vec{q}_{k+1}\)与\(\{ \vec{q}_1, \cdots, \vec{q}_k \}\)正交。
至此，命题得证。Lanczos方案确实可以得到正交归一化的Krylov子空间。

三对角化矩阵

由于上文中说明的\(\vec{q}^\top_k A \vec{q}_{k-2}=0, k \geq 2\)，所以，\(T_k = Q^\top_k A Q_k\)是一个三对角矩阵。\(A\)在\(Q_k\)子空间中的Rayleigh商即\(T_k = Q^\top_k A Q_k\)的Rayleigh商。随着\(k\)的增加，不断对角化\(T_k\)，直到\(T_k\)的本征值收敛。

截断误差

因为截断误差的存在，执行上面的策略以后，也不能保证\(\vec{q}_1, \vec{q}_2, \cdots\)的正交性，所以每次计算出\(A\vec{q}_{k+1}\)，都要依次与\(\vec{q}_i, i=1,\cdots,k-1\)正交化：
\(\vec{q}_{k+1} = \vec{q}_{k+1} - (\vec{q}^\top_i \vec{q}_{k+1}) \vec{q}_i\)
总结下来，Lanczos的做法即：构造\(A\vec{q}_k\)，然后依次剔除其中的\(\vec{q}_k, \vec{q}_{k-1}, \cdots, \vec{q}_1\)的组分，注意这个剔除顺序是逆序的，逆序优于顺序，因为顺序可能出现大数减小数的情况，扣不掉截断误差。

5. Lanczos 代码

下面是我写的Lanczos代码，考虑了截断误差，不变子空间等问题。但是可读性比较差，因为没有时间所有也没有仔细优化。如果在实际问题中优化，主要是针对 矩阵x矢量 的部分优化。
/*
* Lanczos(...) Lanczos method to diagonalize symmetric matrix A[n][n], the output is eigenvalues in z[] and eigenvectors in eigvec[]
* int n dimmension of the real symmetric matrix to be diagonalized
* int nkeep number of states to be kept, Lanczos(...) stops iterating if the lowest nkeep states converge
* double precision precision control parameter, it is used in
* error criterion = precision * max|A_{ij}|
* error = sqrt( 1/nkeep * sum^{nkeep}{i=0} (z^{k} [i] - z^{k-1} [i] )^2 )
* invariant subspace criterion: | q | < precision * max|A_{ij}|
*/
void Lanczos(int n, double **A, int nkeep, double *z, double *eigvec, double precision){

    if(nkeep >n) nkeep = n;

    for(int i=0;i<n;i++) z[i] = 0;

    double maxele=0;
    for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                    if( maxele < fabs( A[i][j] ) )
                            maxele = fabs( A[i][j] );
    int i, j, k, l, maxiter=nkeep * 500;
    int upperk = n>maxiter ? maxiter : n;
    double y, error = precision * maxele * 1E6, criterion = precision * maxele;
    double * eigket = new double [ upperk * upperk ];
    double *d = new double [upperk];
    double *e = new double [upperk];
    double *dd = new double [upperk];
    double *ee = new double [upperk];
    bool finish = false;
    vec omega, v; omega.init(n); v.init(n);

    vec * q = new vec [upperk];
    for(int i=0;i<upperk;i++) q[i].init(n);

    k=0;
    while( k < upperk && error> criterion ){
            // get new Lanczos pivot qk
            omega.SetZero();
            while( sqrt( omega.InnerProduct( omega ) ) < precision ){
                    omega.RandomNormal(); omega.ProjectOut(k, q);
                    if( sqrt( omega.InnerProduct( omega ) ) < precision )
                            continue;
                    omega.Normalize();
            }

            while( error > criterion && k < upperk ){
                    // check orthogonalization
                    v.copy(1,omega); v.ProjectOut(k,q);
                    //cout<<"After projection, ||qk||^2 = "<<v.InnerProduct(v)<<endl;
                    // When it comes here, omega is the fresh lanczos vector
                    q[k].copy(1, omega); k++;
                    // diagonalize tridiagonal matrix
                    v.Av(A, omega);
                    d[k-1] = v.InnerProduct( omega ); // tridiagonal element
                    if(k==1) e[k-1] = 0;
                    else e[k-1] = v.InnerProduct( q[k-2] );
                    for(int i=0;i<k;i++){
                            dd[i] = d[i]; ee[i] = e[i];
                    }
                    //cout<<"dd: "; for(int i=0;i<k;i++)cout<<dd[i]<<"\t";cout<<endl;
                    //cout<<"ee: "; for(int i=0;i<k;i++)cout<<ee[i]<<"\t";cout<<endl;
                    tqls(k, dd, ee, eigket);// O(k^3)
                    sort_eig(k, eigket, dd);
                    // calculate error
                    if( k > nkeep ){// error = sum of squared deviations
                            error = 0;
                            for(i=0;i<nkeep;i++)
                                    error += (dd[i] - z[i])*(dd[i] - z[i]);
                            error = sqrt(error/nkeep);
                    }

                    if( k % 20 ==0 ){
                            cout<<"Lanczos: after "<<k<<" iteractions, error = "<< error << (error > criterion ? " > ": " < ")<<criterion<<endl;
                            cout<<"\t lowest 10 lanczos eigen values:"<<endl;
                            for(int i=0;i<(k>10?10:k);i++)
                                    cout<<"\t"<<dd[i]<<endl;
                    }

                    if( error < criterion ){
                            finish = true;
                            break;
                    }
                    for(i=0;i<k;i++) z[i] = dd[i];
                    //cout<<"z:"; for(int i=0;i<k;i++)cout<<z[i]<<"\t"; cout<<endl;


                    // construct next Lanczos vector
                    v.ProjectOut(k, q); omega.copy(1, v);
                    if( sqrt(omega.InnerProduct(omega)) < precision * maxele )//invariant subspace
                            break;
                    omega.Normalize();
            }
            if(finish) break;
    }
    if(k==upperk && upperk < n ){
            cout<<"Lanczos: failed to converge after "<<upperk<<"iterations"<<endl;
            exit(1);
    }
    if(n==1){
            z[0] = A[0][0];
            eigket[0] = 1;
    }
    for(i=0;i<nkeep;i++){
            for(j=0;j<n;j++){
                    y=0;
                    for(l=0;l<k;l++){
                            y += eigket[i*k+l] * (q[l].v)[j];
                    }
                    eigvec[i*n+j] = y;
            }
    }
    delete [] d; delete [] e; delete [] dd; delete [] ee; delete [] eigket;
    omega.dele(); v.dele();
    for(int i=0;i<upperk;i++) q[i].dele(); delete [] q;
}

6 性能测试

    #include<iostream>
    using namespace std;
    #include<fstream>
    #include<cmath>
    #include<ctime>

    #include"class.h"
    #include"matrix.h"

    extern "C" void dsyev_(char *a,char *b,int *n,double *A,int *nn,double *e,double *work,int *lwork,int *info);

    int main(){

            int n=100;
            //cout<<"\t n="; cin>>n;

            clock_t tstart = clock();

            double **A = new double * [n];
            for(int i=0;i<n;i++) A[i] = new double [n];

            double *w = new double [n];
            double *v = new double [n];
            double *z = new double [n];
            double *eigvec = new double [n*n];
            clock_t tend = clock();
            cout<<" dynamic memory allocation takes "<< (double)(tend-tstart)/CLOCKS_PER_SEC <<" s"<<endl;

            ifstream fin("OveEven.txt");
            string line;
            getline(fin, line);
            tstart = clock();
            for(int i=0;i<n;i++){
                    for(int j=0;j<n;j++){
                            //fin >> A[i][j];
                            //if(i==j) A[i][j] = 1;
                            //else A[i][j] = 0;
                            //A[i][j] = i*i + j*j;
                            A[i][j] = (double)rand()/RAND_MAX;
                            A[j][i] = A[i][j];
                    }
            }
            tend = clock();
            cout<<" Value assignment of A[][] takes "<< (double)(tend-tstart)/CLOCKS_PER_SEC <<" s"<<endl;
            cout<<"A assigned"<<endl;
            for(int i=0;i<n;i++) v[i] = (double)rand()/RAND_MAX;
            cout<<"assignment of values finished"<<endl;
            double y;
            tstart = clock();
            for(int i=0;i<n;i++){
                    y = 0;
                    for(int j=0;j<n;j++)
                            y += A[i][j] * v[j];
                    w[i] = y;
            }
            tend = clock();
            cout<<" matrix multiplication takes "<< (double)(tend-tstart)/CLOCKS_PER_SEC <<" s"<<endl;

            Lanczos( n, A, 1, z, eigvec, 1E-3 );
            //SimpleLanczos(n, A, 1, z, eigvec, 1E-5 );

            double * a = new double [n*n];
            for(int i=0;i<n;i++){
                    for(int j=0;j<n;j++){
                            a[i*n+j] = A[i][j];
                    }
            }

            char alpha='v', beta='L';
            int lwork = 3 * n, info=0;
            double *e = new double [ n ];
            double *work = new double [ 3*n ];

            if( n <=1000 ){
                    dsyev_(&alpha, &beta, &n, a, &n, e, work, &lwork, &info);
                    cout<<" The lowest 20 dsyev eigenvalues"<<endl;
                    for(int i=0;i<(n>20?20:n);i++)
                            cout<<"\t"<<e[i]<<endl;

                    double error;
                    for(int i=0;i<(n>20?20:n);i++)
                            error += ( e[i] - z[i] ) * ( e[i] - z[i] );
                    cout<<"Lanczos RMSD = "<< sqrt(error/n)<<endl;
            }

            return 0;
    }

我使用 [0,1) 的随机数构成的实对称矩阵，进行测试。如果只求最低本征值，要求最后一次迭代前后结果只差小于 1E-3，则
n = 1E2 时，26 次迭代即收敛。
n = 1E3 时，35 次迭代即收敛。
n = 1E4 时，48 次迭代即收敛。
n = 1E5 时，内存溢出。
所以，对于大矩阵极端本征值，Lanczos还是很管用的。