利用 lapack 函数 dsyev 去线性相关

我们经常遇到线性相关的情况，如果

\[c_1 |1\rangle + c_2 |2\rangle + \cdots + c_n |n\rangle = 0, \]

有 \(\{c_1, c_2, \cdots, c_n\}\) 的非零解，则称 \(\{ |1\rangle, \cdots, |n\rangle\}\) 线性相关。否则就叫线性无关。线性代数里有相关阐释，这个概念与矩阵的秩、特征值、方程求根等等都是相互联系的，所以还挺重要挺基本。

在量子力学里的近似计算里，我们常需要构造基矢，要是这些基矢线性相关，就会有伪态，即没有意义的态。
如果 \(n\) 个基矢可以用 \(m\) 个线性无关的基矢表示，那就需要去除多余的 \(n-m\) 个基矢，这就叫去线性相关。俗话说“筛基”。

下面的代码中，EliminateLinear 函数可以去线性相关，主函数中有个简单的展示例子。

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

extern "C" void dsygvx_(int * ITYPE, char * JOBS, char * RANGE, char * UPLO, int * N, double * A, int * LDA, double * B, int * LDB, double * VL, double * VU, int * IL, int * IU, double * ABSTOL, int * M, double * W, double * Z, int * LDZ, double * WORK, int * LWORK, int * IWORK, int * IFAIL, int * INFO);

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

/*
 * Eliminate linear bases
 * B: Overlap matrix of the kets, on exit it is written over by overlap matrix of only those linearly independent kets. A does also shrink as B does, e.g. A can be the hamiltonian matrix.
 * index: indices of linearly independent kets
 * zero_trun:  if the overlap matrix has eigenvalue < zero_trun * Max|eig|, then the kets are linearly dependent.
 *
 *
 */
void EliminateLinear( int & n, double * A, double * B, int * index,  double zero_trun ){

	double * C = new double [ n*n ];
	double * tempE = new double [ n ];

	int count = 0, coord1, coord2;
	bool flag;

	char a='v';
	char b='L';
	int lwork=n*3;
	int info=0;
	double *work=new double [n*3];   //temporary array

	int ntemp;

	// **** define seq[], so that the seq[i]-th basis has the i-th largest norm

	int * seq = new int [ n ];
	for(int i=0; i<n; i++) seq[i] = i; // initial values of seq[i]: 0,1,2,3,...
	double * ove = new double [ n ];
	for(int i=0; i<n; i++) ove[i] = B[ i*n+i ];
	//cout<<" originally ove: "; for(int i=0;i<n;i++)cout<<ove[i]<<",  "; cout<<endl;
	// bubble sort
	for(int i=n-2;i>=0;i--){
		for(int j=0;j<=i;j++){ // put the j-th smallest at n-1-j th position
			if( ove[seq[j]] < ove[seq[j+1]] ){
				ntemp = seq[j]; seq[j] = seq[j+1]; seq[j+1] = ntemp;
			}
		}
		//cout<<" ove: "; for(int i=0;i<n;i++)cout<< ove[ seq[i] ] <<",  "; cout<<endl;
	}

	//cout<<" seq: "; for(int i=0;i<n;i++)cout<<seq[i]<<",  "; cout<<endl;

	delete [] ove;

	double max_fabs_eig = 0;
	// 按 overlap 从大到小 sort basis，如果 前 k 个 基的 overlap 矩阵本征值里有特别小的，就扔掉第 k 个基
	for(int p=0;p<n;p++){ // basis with the p-th largest norm

		int k = seq[p];
		//cout<<" overlap of the new basis = " << B[ k*n + k ] <<endl;

		// C is the overlap matrix of { index[0], ..., index[count], k }
		for(int i=0;i< count+1;i++){ // include one more row and line

			if( i==count ) coord1 = k; // the new basis
			else coord1 = index[i]; // index[i]: the i-th basis that survived

			for(int j=0;j< count+1;j++){

				if( j==count ) coord2 = k;
				else coord2 = index[j];

				C[ i* (count+1) + j ] = B[ coord1 * n + coord2 ];
			}
		}

		ntemp = count + 1;

		dsyev_( &a, &b, &ntemp, C, &ntemp, tempE, work, &lwork, &info );

		max_fabs_eig = 0; // find the maximum |eig| of the overlap matrix C
		for(int i=0;i<count+1;i++) 
			if( max_fabs_eig < fabs(tempE[i]) ) max_fabs_eig = fabs(tempE[i]);

		flag = true;
		for(int i=0;i< count+1;i++){//if there's 0 or smaller eigenvalue, kick out the bad ket
			//cout<<tempE[i]<<", ";
			if( tempE[i] < zero_trun * max_fabs_eig ){
				flag = false; break;
			}
		}
		cout<<" --------------------------------------------- \n";
		cout<<"sort: p = "<<p<<" \t eig = "; for(int i=0;i<ntemp;i++)cout<<tempE[i]<<"  "; cout<<endl;
		if( flag ) cout<<" \t p-th ket is saved. \n";
		else cout<<" \t p-th ket is abandoned. \n";

		if(flag){ // otherwise add up this new basis
			index[ count ] = k; count ++;
		}
		//cout<<"count = "<<count<<"\t eigenvalues: ";
		//for(int i=0;i< count;i++) cout<< tempE[i]<<"\t"; cout<<endl;
	}
	delete [] seq;

	// B shrinks from original dimension to smaller (linearly independent) dimension
	for(int i=0;i<count;i++)
		for(int j=0;j<count;j++)
			C[ i*count + j ] = B[ index[i]*n + index[j] ];
	for(int i=0;i<count;i++)
		for(int j=0;j<count;j++)
			B[ i * count + j ] = C[ i*count + j];

	// A shrinks from original dimension to smaller (linearly independent) dimension
	for(int i=0;i<count;i++)
		for(int j=0;j<count;j++)
			C[ i*count + j ] = A[ index[i]*n + index[j] ];
	for(int i=0;i<count;i++)
		for(int j=0;j<count;j++)
			A[ i * count + j ] = C[ i*count + j];

	delete [] tempE; delete [] C; delete [] work;

	n = count;
};


double InnerProduct(int n, double *a, double *b){
	double y = 0;
	for(int i=0;i<n;i++) y += a[i] * b[i];
	return y;
}

int main(){

	double ket[4][3]={ {1,2,3}, {2,3,4}, {3,4,5}, {4,5,6} };
	double H[16];
	double Ove[4*4];
	for(int i=0;i<4;i++)
		for(int j=0;j<4;j++)
			Ove[i*4+j] = InnerProduct(3, ket[i], ket[j]);
	int n = 4, index[4];
	EliminateLinear( n, H, Ove, index, 1E-5 );
	cout<<" n = "<<n<<endl;
	cout<<" index: "; for(int i=0;i<n;i++)cout<<index[i]<<"  "; cout<<endl;
	
	return 0;
}

这文件我命名为 EliminateLinear.cpp，运行：

icc EliminateLinear.cpp -qmkl
./a.out

运行结果：

 --------------------------------------------- 
sort: p = 0 	 eig = 77  
 	 p-th ket is saved. 
 --------------------------------------------- 
sort: p = 1 	 eig = 0.0472617  126.953  
 	 p-th ket is saved. 
 --------------------------------------------- 
sort: p = 2 	 eig = 7.60943e-15  0.231112  155.769  
 	 p-th ket is abandoned. 
 --------------------------------------------- 
sort: p = 3 	 eig = 1.61181e-14  0.598283  140.402  
 	 p-th ket is abandoned. 
 n = 2
 index: 3  2  

这个例子里，我们给了 4 个 3 维矢量：

\[\left( \begin{matrix} 1 \\ 2\\ 3 \end{matrix} \right) \left( \begin{matrix} 2 \\ 3\\ 4 \end{matrix} \right) \left( \begin{matrix} 3 \\ 4\\ 5 \end{matrix} \right) \left( \begin{matrix} 4 \\ 5\\ 6 \end{matrix} \right). \]

很显然，这些矢量的秩为 2，因为只需要 2 个就可以表示其他。相邻的任意两个相减得到

\[\left( \begin{matrix} 1 \\ 1\\ 1 \end{matrix} \right). \]

筛出来结果确实也是 \(n=2\)，存活的基矢是第 \(3,2\) 个，也就是最后 2 个（c++序号从0开始）。
取最后 2 个，是因为我代码里优先挑模方比较大的基矢。