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

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

Rayleigh商定义：

对于nxn实对称阵A，任意非零n维矢量(vec{x} eq 0)，定义Rayleigh商为
(r_A( vec{x} ) = frac{ vec{x}^ op A vec{x} }{ vec{x}^ op 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)，有(Avec{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})的梯度
( abla r_A(vec{x}) = frac{2}{vec{x}^ op vec{x}}( A vec{x} - r_A(vec{x})vec{x} ).)
取( abla 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)个（(kleq n)）正交归一的n维矢量({ vec{q}_1, vec{q}_1, cdots vec{q}_k })，构成子空间，记nxk的矩阵(Q_k equiv egin{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})^ op A (Q_k vec{y}) }{ (Q_k vec{y})^ op (Q_k vec{y}) } = frac{ vec{y}^ op (Q^ op_k A Q_k) vec{y} }{ vec{y}^ op vec{y} } = r_{Q^ op_k A Q_k}(vec{y}))。
转化为(Q^ op_k A Q_k)的Rayleigh商。

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

根据上面的分析，要想构造很小的子空间，使得该子空间内存在最大/最小本征值对应的本征矢(vec{v}_1, vec{v}_n)的近似矢量，则需要构造(Q_k)，使得(r_{Q^ op_k A Q_k})的最小值尽量小，最大值尽量大。
我们可以从选定的(vec{q}_1)开始，逐渐增加(vec{q}_2, vec{q}_3, ...)，逐渐逼近目标。
一个有用的性质是：( abla r_A(vec{q}_1) = frac{2}{vec{q}^ op_1 vec{q}_1}( A vec{q}_1 - r_A(vec{q}_1)vec{q}_1 ))，这个梯度是(Avec{q}_1, vec{q}_1)的线性叠加。
于是，美妙的想法产生了：何不如此构造子空间呢：(kappa(A, vec{q}_1, k) = { vec{q}_1, Avec{q}_1, A^2vec{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, Avec{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^ op A vec{q}_1 ) vec{q}_1, vec{q}_2 = vec{r}_2 / | vec{r}_2 |)。即从(Avec{q}_1)中剔除(vec{q}_1)的分量，然后归一化。
对于(kgeq 2)：(vec{r}_{k+1} = A vec{q}_k - (vec{q}^ op_k A vec{q}_k) vec{q}_k - (vec{q}^ op_{k-1} A vec{q}_k) vec{q}_{k-1}, vec{q}_{k+1} = vec{r}_{k+1} / | vec{r}_{k+1} |)。即从(Avec{q}_k)中剔除(vec{q}_{k-1}, vec{q}_k)的分量，然后归一化。

正交化吗？归纳法证明

显然，如上构造的子空间({vec{q}_1, vec{q}_2} = {vec{q}_1, Avec{q}_1 }, {vec{q}_1, vec{q}_2, vec{q}_3 } = { vec{q}_1, Avec{q}_1, A^2vec{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, Avec{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, Avec{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, Avec{q}_1, cdots, A^{k}vec{q}_1 })

由于(Avec{q}_1, Avec{q}_2, cdots, Avec{q}_{k-2})都是子空间({ vec{q}_1, vec{q}_2, cdots, vec{q}_{k-1} })中的矢量，所以都垂直于(vec{q}_k)，即
(vec{q}_k^ op A vec{q}_{k-2} = vec{q}^ op_k A vec{q}_{k-3} = cdots vec{q}^ op_k A vec{q}_1 = 0)。
构造(vec{r}_{k+1} = A vec{q}_k - (vec{q}^ op_k A vec{q}_k) vec{q}_k - (vec{q}^ op_{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, Avec{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})是(Avec{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)与(Avec{q}_1, Avec{q}_2, cdots, Avec{q}_{k-2})都正交，即
(vec{q}^ op_k A vec{q}_1 = vec{q}^ op_k A vec{q}_2 = cdots = vec{q}^ op_k A vec{q}_{k-2} = 0)。
利用A的对称性，上式即
(vec{q}^ op_1 A vec{q}_k = vec{q}^ op_2 A vec{q}_k = cdots = vec{q}^ op_{k-2} A vec{q}_{k} = 0)。

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

三对角化矩阵

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

截断误差

因为截断误差的存在，执行上面的策略以后，也不能保证(vec{q}_1, vec{q}_2, cdots)的正交性，所以每次计算出(Avec{q}_{k+1})，都要依次与(vec{q}_i, i=1,cdots,k-1)正交化：
(vec{q}_{k+1} = vec{q}_{k+1} - (vec{q}^ op_i vec{q}_{k+1}) vec{q}_i)
总结下来，Lanczos的做法即：构造(Avec{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{k+1} | < 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]<<"	";cout<<endl;
                    //cout<<"ee: "; for(int i=0;i<k;i++)cout<<ee[i]<<"	";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<<"	 lowest 10 lanczos eigen values:"<<endl;
                            for(int i=0;i<(k>10?10:k);i++)
                                    cout<<"	"<<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]<<"	"; 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<<"	 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<<"	"<<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还是很管用的。