反对称实矩阵正则化

之前构造了一个办法来正则化反对称实矩阵：https://www.cnblogs.com/luyi07/p/15578164.html
那个办法只是初期尝试做的，它有漏洞。只要反对称实矩阵的特征值有简并，那个办法就失败。
所以现在重新做一个办法来正则化反对称实矩阵。

1. 反对称实矩阵特征值的特点

1.1 反对称实矩阵的非０特征值一定是纯虚数

设有反厄米矩阵 \(K\)，反厄米即 \(K^\dagger = -K\)，可以定义

\[H = iK, \]

则有

\[H^\dagger = - i K^\dagger = i K = H, \]

即 \(H\) 是厄米矩阵。众所周知，厄米矩阵的特征值是实数，即存在幺正矩阵 \(U\) 和实数对角阵 \(\Lambda\), 使得

\[H = U \Lambda U^\dagger， \]

那么，

\[K = -i H = U(-i\Lambda) U^\dagger, \]

这说明 \((-i\Lambda)\) 的对角元是 \(K\) 的特征值，\(U\) 的列向量是相应的本征向量。所以 \(K\) 的非０本征值都是纯虚数。
反对称实矩阵是反厄米矩阵的特殊情况，所以反对称实矩阵的非０特征值也是纯虚数。

1.2 反对称实矩阵的非０特征值一定是两两一组，每组两个特征值互为共轭，两个特征向量也互为共轭

假设 \(K\) 有特征值 \(\lambda\)，有

\[K \vec{z} = \lambda \vec{z}, \]

则有

\[( K \vec{z}^* )^\dagger = (K \vec{z})^\top = \lambda \vec{z}^\top, \]

再将上式左右取厄米共轭，得到

\[K \vec{z}^* = (\lambda \vec{z}^\top)^\dagger = \lambda^* \vec{z}^*, \]

这说明，\((\lambda^*, \vec{z}^*)\) 也是一个本正对（即本征值+本征向量）。
前面证明了 \(\lambda\) 是纯虚数，所以 \(\lambda \neq 0\) 时，\(\lambda^* \neq \lambda\)。
因此，反对称实矩阵的非零特征值一定是两两一组，每组两个特征值互为共轭。

2. 反对称实矩阵正则化

2.1 非０本征值无简并的情况

我们先讨论简单的情况，即反对称实矩阵的非０本征值没有简并，即下面的式子中 \(\lambda_i \neq \lambda_j, i\neq j\).
综合前面的分析，任意反对称实矩阵 \(K\) 可以写作

\[K = [ \vec{z}_1, \vec{z}^*_1, \cdots] diag[ \lambda_1, \lambda^*_1, \cdots, \lambda_k, \lambda^*_k, 0, \cdots, 0 ] [ \vec{z}_1, \vec{z}^*_1, \cdots]^\dagger \]

其中 \(\vec{z}_1, \vec{z}^*_1, \cdots\) 为正交归一基，且 \(K \vec{z} = \lambda \vec{z},~~ K\vec{z}^* = \lambda^* \vec{z}^*\)，且 \(\lambda \neq 0\)，记

\[\vec{z} = \frac{\vec{x} + i \vec{y}}{\sqrt{2}},~~~~~~~ \vec{x}, \vec{y} \in R^n \]

由归一性，

\[\vec{z}^* \cdot \vec{z} = (\vec{x} - i \vec{y}) \cdot (\vec{x} + i \vec{y})/2 = \frac{ |\vec{x}|^2 + |\vec{y}|^2 }{2} = 1, \Rightarrow |\vec{x}|^2 + |\vec{y}|^2 = 2. \]

由于 \(\lambda \neq \lambda^*\), 所以 \(\vec{z}, \vec{z}^*\) 必须正交——注意，复矢量内积定义为\((\vec{\alpha}, \vec{\beta}) = \vec{\alpha}^* \cdot \vec{\beta}\)——即

\[\vec{z} \cdot \vec{z} = 0 \Rightarrow (\vec{x} - i \vec{y}) \cdot (\vec{x} + i \vec{y}) = 0 \Rightarrow \{ \vec{x} \cdot \vec{y} = 0, \vec{x}\cdot \vec{x} - \vec{y} \cdot \vec{y} = 0 \}, \]

所以有

\[|\vec{x}| = |\vec{y}| = 1, ~~~ \vec{x} \perp \vec{y}. \]

另外

\[\left\{\vec{z}_k, \vec{z}^*_k\right\} \perp \left\{ \vec{z}_l, \vec{z}^*_l \right\} \Rightarrow \left\{ \vec{x}_k, \vec{y}_k \right\} \perp \left\{ \vec{x}_l, \vec{y}_l \right\}, \]

所以，以下基矢是正交归一基：

\[\left\{ \vec{x}_1, \vec{y}_1, \cdots, \vec{x}_k, \vec{y}_k, \vec{w}_{2k+1}, \cdots, \vec{w}_n \right\}, \]

其中，\(\vec{w}_{2k+1}, \cdots, \vec{w}_n\)　是人为补齐的实数矢量，对应值为 \(0\) 的特征值。
这套正交归一基有如下结构：

\(\left\{ \vec{x}_i, \vec{y}_i \right\}\) 构成 \(K\) 的不变子空间。
\(\left\{ \vec{w}_{2k+1}, \cdots, \vec{w}_n \right\}\) 构成 \(K\) 的不变子空间。

所以，如果用这套正交归一基构造幺正变换，会得到块对角形式

\[K [ \vec{x}_1, \vec{y}_1, \cdots, \vec{x}_k, \vec{x}_y, \vec{w}_{2k+1}, \cdots, \vec{w}_n] \\ =\left[ \vec{x}_1, \vec{y}_1, \cdots, \vec{x}_k, \vec{x}_y, \vec{w}_{2k+1}, \cdots, \vec{w}_n \right] \left[ [\begin{smallmatrix} 0 & -i\sqrt{2}\lambda_1 \\ i\sqrt{2}\lambda_1 & 0 \end{smallmatrix}], \cdots, [\begin{smallmatrix} 0 & -i\sqrt{2}\lambda_k \\ i\sqrt{2}\lambda_k & 0 \end{smallmatrix}], 0, \cdots, 0 \right]. \]

由于 \(\lambda_1, \cdots, \lambda_k\) 都是纯虚数，所以这样就把 \(K\) 正则化了。

\[[ \vec{x}_1, \vec{y}_1, \cdots, \vec{x}_k, \vec{x}_y, \vec{w}_{2k+1}, \cdots, \vec{w}_n]^\top K [ \vec{x}_1, \vec{y}_1, \cdots, \vec{x}_k, \vec{x}_y, \vec{w}_{2k+1}, \cdots, \vec{w}_n] \\ = \left[ [\begin{smallmatrix} 0 & -i\sqrt{2}\lambda_1 \\ i\sqrt{2}\lambda_1 & 0 \end{smallmatrix}], \cdots, [\begin{smallmatrix} 0 & -i\sqrt{2}\lambda_k \\ i\sqrt{2}\lambda_k & 0 \end{smallmatrix}], 0, \cdots, 0 \right]. \]

2.2 非０本征值有简并的情况

非０本征值有简并的时候，就存在子空间。以２重简并为例（更高的简并是同理），假设 \(\lambda\) 有两重简并，那么根据前面的分析，\(-\lambda\) 自然也是２重简并。

\[K \vec{z}_1 = \lambda \vec{z}_1, ~~ K \vec{z}_2 = \lambda \vec{z}_2, ~~ \vec{z}_1 \perp \vec{z}_2, |\vec{z}_1| = |\vec{z}_2| = 1; \\ K \vec{z}_3 = -\lambda \vec{z}_3, ~~ K \vec{z}_4 = - \lambda \vec{z}_3, ~~ \vec{z}_3 \perp \vec{z}_4, |\vec{z}_3| = |\vec{z}_4| = 1.\\ K[\vec{z}_1, \vec{z}_2, \vec{z}_3, \vec{z}_4] = [\vec{z}_1, \vec{z}_2, \vec{z}_3, \vec{z}_4]diag\{\lambda, \lambda, -\lambda, -\lambda\}. \\ \vec{z}_i = (\vec{x}_i + i \vec{y}_i)/\sqrt{2}. \]

2.2.1 \(\vec{x}_i, \vec{y}_i\) 这两个矢量正交归一

假设其对应的特征值是 \(\mu\)，则有

\[K \vec{z}_i = \mu \vec{z}_i, ~~ \Rightarrow ~~ K \vec{z}^*_i = - \mu \vec{z}^*_i, ~~ \Rightarrow ~~ \vec{z}_i \perp \vec{z}^*_i \]

因此有

\[(\vec{z}^*_i, \vec{z}_i) = \frac{1}{2}( |\vec{x}_i|^2 - |\vec{y}_i|^2 + 2i \vec{x}_i \cdot \vec{y}_i ) = 0, ~~\Rightarrow~~ |\vec{x}_i| = |\vec{y}_i|, \vec{x}_i \perp \vec{y}_i. \]

另外有归一化条件

\[(\vec{z}_i, \vec{z}_i) = \frac{1}{2}( |\vec{x}|^2 + |\vec{y}|^2 ) = 1, \]

所以有

\[|\vec{x}_i| = |\vec{y}_i| = 1, \vec{x}_i \perp \vec{y}_i. \]

2.2.2 \(\vec{x}_1, \vec{y}_1, \vec{x}_2, \vec{y}_2\)正交

\(\vec{z}_1 \perp \vec{z}_2\)，这是我们知道的；另外，由于 \(\vec{z}^*_2\) 对应 \(-\lambda\)，所以也有 \(\vec{z}_1 \perp \vec{z}^*_2\)；所以，自然有 \(\vec{z}_1 \perp \{ \vec{x}_2, \vec{y}_2 \}\)。

另外，\(\vec{z}_1 \perp \vec{z}_2 \Rightarrow \vec{z}^*_1 \perp \vec{z}^*_2\)；由于 \(\vec{z}_2\) 对应 \(\lambda\)，所以也有 \(\vec{z}^*_1 \perp \vec{z}_2\)；所以，有 \(\vec{z}^*_1 \perp \{\vec{x}_2, \vec{y}_2\}\)。

而前面已经分析得到 \(\vec{z}_1 \perp \vec{z}^*_1\)，所以有 \(\{x_1, y_1\} \perp \{\vec{x}_2, \vec{y}_2\}\)。

这样可以分析得到 \(\{ \vec{x}_1, \vec{y}_1, \vec{x}_2, \vec{y}_2\}\)是一组正交归一实数基，张成一个不变子空间；类似地，\(\{ \vec{x}_3, \vec{y}_3, \vec{x}_4, \vec{y}_4\}\)也是一组正交归一实数基，也可以张成同一个不变子空间。\(\vec{x}_3\) 不一定等于 \(\vec{x}_1\) 或 \(\vec{x}_2\)。

对于反对称实矩阵的这个不变子空间，取 \(\{ \vec{x}_1, \vec{y}_1, \vec{x}_2, \vec{y}_2\}\) 这组基，即在这个子空间内完成正则化：

\[K[\vec{x}_1, \vec{y}_1, \vec{x}_2, \vec{y}_2] = [\vec{x}_1, \vec{y}_1, \vec{x}_2, \vec{y}_2] \left[ [\begin{smallmatrix} 0 & -i\sqrt{2}\lambda_1 \\ i\sqrt{2}\lambda_1 & 0 \end{smallmatrix}], [\begin{smallmatrix} 0 & -i\sqrt{2}\lambda_2 \\ i\sqrt{2}\lambda_2 & 0 \end{smallmatrix}] \right]. \]

3. 数值验证

根据上面的理论推导，我们得到结论：只要能把反对称实矩阵 \(K\) 对角化，就能通过实正交变换把它正则化，变为实数正则矩阵。
如何对角化反对称实矩阵呢？GSL 有现成的函数，只是要稍微变一下。

3.1 反对称实矩阵 \(\rightarrow\) 厄米复数矩阵

反对称实矩阵 \(K\) 只要乘上 \(i\)，就是一个厄米矩阵：

\[H = i K \Rightarrow H^\dagger = H. \]

3.2 GSL库函数：对角化厄米复数矩阵

通过翻阅 GSL 的文档，可以发现　gsl_eigen_hermv　函数，可以对角化复数厄米矩阵。所以我写了个函数做这个，如下。

#include<gsl/gsl_complex.h>
#include<gsl/gsl_complex_math.h>
#include<gsl/gsl_math.h>
#include<gsl/gsl_eigen.h>

#include<complex>

/*
 * GSL_eigen_hermv: given complex Hermitian matrix H[], get eigenvalues eig[] and eigenvectors eigvec[]
 * in respec to eig[i], its eigenvector is eigvec[i*n + x ], x = 0,...,n-1
 */
void GSL_eigen_hermv( int n, complex<double> * H, double * eig, complex<double> * eigvec ){

        gsl_eigen_hermv_workspace * w = gsl_eigen_hermv_alloc( n );

        gsl_matrix_complex * A = gsl_matrix_complex_alloc( n, n );
        for(int i=0;i<n;i++)
                for(int j=0;j<n;j++)
                        gsl_matrix_complex_set( A, i, j, gsl_complex_rect( H[i*n+j].real(), H[i*n+j].imag() ) );
        gsl_matrix_complex * evec = gsl_matrix_complex_alloc( n, n );
        gsl_vector * eval = gsl_vector_alloc( n );

        cout<< " gsl_eigen_herm returns "<< gsl_eigen_hermv( A, eval, evec, w );
        cout<<", it seems 0 means OK. \n";

        gsl_eigen_hermv_sort( eval, evec, GSL_EIGEN_SORT_ABS_ASC );

        for(int i=0;i<n;i++){

                double eval_i = gsl_vector_get( eval, i ); eig[i] = eval_i;
                cout<<" eigenvalue = "<<eval_i<<endl;
                gsl_vector_complex_view evec_i = gsl_matrix_complex_column( evec, i );
                gsl_vector_complex_fprintf( stdout, &evec_i.vector, "%g");

                for(int j=0;j<n;j++){
                       gsl_complex x = gsl_matrix_complex_get( evec, j, i );
                        eigvec[i*n+j].real( GSL_REAL(x) ); eigvec[i*n+j].imag( GSL_IMAG(x) );
                }

        }

        gsl_eigen_hermv_free( w );
        gsl_matrix_complex_free( A );
        gsl_vector_free( eval );
        gsl_matrix_complex_free( evec );
}

然后写了个测试主函数，使用 GSL 文档上的实数对称阵测试例子，进行测试。

int main(){

        double data[] = { 0.0 , 1/2.0, 1/3.0, 1/4.0,
-1/2.0, 0.0, 1/4.0, 1/5.0,
-1/3.0, -1/4.0, 0.0, 1/6.0,
-1/4.0, -1/5.0, -1/6.0, 0.0 };

        complex<double> * H = new complex<double> [ 4*4 ];
        complex<double> x; x.real(0); x.imag(1);
        cout<<" x = "<< x <<endl;
        for(int i=0;i<16;i++) H[i] = data[i] * x ;

        cout<<" H : "<<endl;
        for(int i=0;i<4;i++){
                for(int j=0;j<4;j++)cout<<H[i*4+j]<<"   "; cout<<endl;
        }

        double * eig = new double [ 4 ];

        complex<double> * eigvec = new complex<double> [ 4*4 ];

        GSL_eigen_hermv( 4, H, eig, eigvec );

        cout<<" eig: "; for(int i=0;i<4;i++)cout<<eig[i]<<","; cout<<endl;
        cout<<" eigvec: "<<endl;
        for(int i=0;i<4;i++){
                for(int j=0;j<4;j++)cout<<eigvec[i*4+j]<<","; cout<<endl;
        }

        return 0;
}

命名为 test_complex_Hermitian.cpp，进行编译运行

g++ test_complex_Hermitian.cpp -lgsl -lgslcblas
./a.out

得到结果

 x = (0,1)
 H : 
(0,0)   (0,0.5)   (0,0.333333)   (0,0.25)   
(-0,-0.5)   (0,0)   (0,0.25)   (0,0.2)   
(-0,-0.333333)   (-0,-0.25)   (0,0)   (0,0.166667)   
(-0,-0.25)   (-0,-0.2)   (-0,-0.166667)   (0,0)   
 gsl_eigen_herm returns 0, it seems 0 means OK. 
 eigenvalue = 0.1075
-0.334421 0
0.375581 0.194331
-0.234738 -0.515822
-0.438178 0.442903
 eigenvalue = -0.1075
-0.334421 0
0.375581 -0.194331
-0.234738 0.515822
-0.438178 -0.442903
 eigenvalue = -0.736432
0.623027 0
0.2016 0.529653
-0.126 0.40367
-0.2352 0.237736
 eigenvalue = 0.736432
0.623027 0
0.2016 -0.529653
-0.126 -0.40367
-0.2352 -0.237736
 eig: 0.1075,-0.1075,-0.736432,0.736432,
 eigvec: 
(-0.334421,0),(0.375581,0.194331),(-0.234738,-0.515822),(-0.438178,0.442903),
(-0.334421,0),(0.375581,-0.194331),(-0.234738,0.515822),(-0.438178,-0.442903),
(0.623027,0),(0.2016,0.529653),(-0.126,0.40367),(-0.2352,0.237736),
(0.623027,0),(0.2016,-0.529653),(-0.126,-0.40367),(-0.2352,-0.237736),

与 GSL 文档上的参考结果一致。

3.3 自编函数：正则化反对称实矩阵

根据前面的理论推断，对于给定的反对称实矩阵 \(K\)，我们可以如下正则化：

对角化 \(H = iK\)
将 \(H\) 的所有本征对——即（本征值，本征矢）这样的对——按本征值的绝对值降序重排
取所有特征值小于０的特征矢量，录用其实部向量、虚部向量：

\[\vec{x}_k = \frac{ \vec{z}_k + \vec{z}^*_k }{\sqrt{2}}, \\ \vec{y}_k = \frac{ \vec{z}_k - \vec{z}^*_k }{ \sqrt{2} i}. \]

　　根据前面的理论分析，这些实部、虚部向量一定是正交归一的。

补齐正交归一基。即补齐 \(w_{2k+1}, \cdots w_n\)。这可以自己写一个函数，给定 \(m\) 个正交归一基，补齐剩下 \(n-m\) 个正交归一基。
最后，得到新的实数正交归一基以后，用这些基做正交变换，得到 \(K\) 的正则形式。

写出来的代码如下：


void printmtx( int n, double * U ){
        for(int i=0;i<n;i++){
                for(int j=0;j<n;j++) cout<<U[i*n+j]<<"  "; cout<<endl;
        }
}

void RandomNormal(int n, double * v){
        double y=0;
        for(int i=0;i<n;i++){
                v[i] = ((double)rand())/RAND_MAX ;
                y += v[i] * v[i];
        }
        for(int i=0;i<n;i++)
                v[i] /= sqrt(y);
};

void PickOut( int n, double * x, double * y ){
        double innerproduct = 0; // inner product
        for(int i=0;i<n;i++) innerproduct += x[i] * y[i];
        for(int i=0;i<n;i++) x[i] -= innerproduct * y[i];
}

void Normalize( int n, double * x ){
        double y = 0; for(int i=0;i<n;i++) y += x[i] * x[i]; y = sqrt(y);
        for(int i=0;i<n;i++) x[i] /= y;
}

// Given U[0], U[1], ..., U[m-1] are ortho-normal, find the other n-m bases, so that U[0], ..., U[n-1] are complete and ortho-normal
void complete_orthonormal_basis( int n, int m, double * U ){
        int count = m; double * y = new double [n];
        while( count < n ){
                for(int i=0;i<n;i++) y[i] = 0; double normy = 0;
                while (normy < 1E-4){
                        RandomNormal( n, y );
                        for(int i=0;i<count;i++) PickOut( n, y, U+i*n);
                        normy = 0; for(int i=0;i<n;i++) normy += y[i] * y[i];
                }
                Normalize(n, y); for(int i=0;i<n;i++) U[count * n + i] = y[i];
                count ++;
        }
}

/*
 * ARMC: transform an Anti-Real-Matrix into Canonical form
 * 
 * M = U^\top Mcanonical U
 */
void ARMC( int n, double * M, double * U, double * Mcanonical ){

        complex<double> * H = new complex<double> [ n*n ];
        complex<double> x; x.real(0); x.imag(1);
        for(int i=0;i<n*n;i++) H[i] = M[i] * x ;

        double * eig = new double [ n ];

        complex<double> * eigvec = new complex<double> [ n*n ];

        GSL_eigen_hermv( n, H, eig, eigvec );

        double * xtemp = new double [n]; double * ytemp = new double [n];

        int count = 0;
        for(int i=0;i<n;i++){
                if( eig[i] < 0 && fabs(eig[i]) > 1E-10 * fabs(eig[0]) ){ // 1E-10 * fabs(eig[0]) is used as zero (small enough)
                        for(int j=0;j<n;j++){
                                xtemp[j] = eigvec[i*n+j].real() * sqrt(2);
                                ytemp[j] = eigvec[i*n+j].imag() * sqrt(2);
                                U[count*n+j] = xtemp[j]; U[(count+1)*n+j] = ytemp[j];
                        }
                        count+=2;
                }
        }

        if( count < n ) complete_orthonormal_basis( n, count, U );

        // UM
        double * UM = new double [ n*n ];
        for(int i=0;i<n;i++)
        for(int j=0;j<n;j++){
                double x = 0;
                for(int k=0;k<n;k++) x += U[i*n+k] * M[k*n+j];
                UM[i*n+j] = x;
        }
        // Mcanonical = (UM) * U^\top 
        for(int i=0;i<n;i++)
        for(int j=0;j<n;j++){
                double x = 0;
                for(int k=0;k<n;k++) x += UM[i*n+k] * U[j*n+k];
                Mcanonical[i*n+j] = x;
        }
        delete [] xtemp; delete [] ytemp;
        delete [] UM; delete [] eigvec; delete [] eig; delete [] H;
}

测试：

加了几行测试输出代码，正则化 12*12 反对称阵：

0 0.0582782 0.0809623 0 0.242372 0.148711 0.0747338 0.1232 0.00109344 0.0352467 0.142309 0.0586949
-0.0582782 0 0 0.0809623 0.130991 0.0181291 0.141333 0.00607454 0.146359 -0.0488136 0.16139 0.00949629
-0.0809623 -0 0 -0.0582782 0.0488136 0.146359 -0.00607454 0.141333 -0.0181291 0.130991 -0.00949629 0.16139
-0 -0.0809623 0.0582782 0 0.0352467 -0.00109344 0.1232 -0.0747338 0.148711 -0.242372 0.0586949 -0.142309
-0.242372 -0.130991 -0.0488136 -0.0352467 0 0.17432 0.165275 0.123782 0.257831 0 0.0425534 0.19093
-0.148711 -0.0181291 -0.146359 0.00109344 -0.17432 0 -0.0727708 0.068801 0 0.257831 -0.098141 0.235799
-0.0747338 -0.141333 0.00607454 -0.1232 -0.165275 0.0727708 0 0 0.068801 -0.123782 -0.151594 -0.00598481
-0.1232 -0.00607454 -0.141333 0.0747338 -0.123782 -0.068801 -0 0 0.0727708 0.165275 -0.00598481 0.151594
-0.00109344 -0.146359 0.0181291 -0.148711 -0.257831 -0 -0.068801 -0.0727708 0 -0.17432 -0.235799 -0.098141
-0.0352467 0.0488136 -0.130991 0.242372 -0 -0.257831 0.123782 -0.165275 0.17432 0 0.19093 -0.0425534
-0.142309 -0.16139 0.00949629 -0.0586949 -0.0425534 0.098141 0.151594 0.00598481 0.235799 -0.19093 0 0
-0.0586949 -0.00949629 -0.16139 0.142309 -0.19093 -0.235799 0.00598481 -0.151594 0.098141 0.0425534 -0 0

得到

  U: 
-0.0303352  -0.323071  0.284728  -0.475271  -0.0864182  0.457386  0.0200367  0.146637  0.0431553  -0.472297  -0.32487  -0.136787  
0  0.0853881  -0.0450322  -0.214972  -0.176651  0.184243  -0.378003  0.260995  -0.432414  0.502859  -0.339272  0.336687  
-0.521628  -0.240866  -0.25917  0.0276394  -0.223088  0.138885  0.241166  0.137526  0.492668  0.151539  0.0141235  0.435819  
0  -0.158371  -0.210943  0.0125017  -0.652812  -0.411771  -0.177369  -0.352668  0.0206277  -0.125337  -0.363139  -0.175236  
0.511228  -0.15273  0.308259  -0.156083  -0.122579  -0.0502583  0.0552624  0.0744275  0.528361  0.515782  -0.0123233  -0.156335  
0  0.3574  0.375394  -0.348846  -0.290007  -0.182507  0.524695  -0.205627  -0.191796  -0.0316432  0.116331  0.351185  
0.382172  0.468555  -0.263858  0.20879  -0.181767  0.0407166  0.1573  0.501511  0.187118  -0.31478  -0.256712  0.101154  
0  0.128064  0.285378  0.466648  -0.483729  0.560618  -0.151918  -0.163847  -0.0286619  0.00655211  0.286129  -0.0587596  
-0.527576  0.533958  0.158684  -0.0908641  0.0293914  -0.00837181  -0.0171707  0.0710162  0.194925  0.195811  -0.241528  -0.512025  
0  -0.281206  -0.169889  0.181609  -0.120986  0.103253  0.630479  0.19511  -0.409262  0.255714  -0.0895545  -0.396615  
-0.203072  -0.222936  0.490404  0.236063  -0.141072  -0.453226  -0.142712  0.571527  -0.0960864  -0.12135  0.125592  0.027982  
0  0.0658956  -0.351698  -0.471816  -0.289534  -0.00879994  -0.150812  0.266751  -0.0387136  -0.02785  0.635374  -0.256072  
 U * U^+ : 
1   0   0   0   0   0   0   0   0   0   0   0   
0   1   0   0   0   0   0   0   0   0   0   0   
0   0   1   0   0   0   0   0   0   0   0   0   
0   0   0   1   0   0   0   0   0   0   0   0   
0   0   0   0   1   0   0   0   0   0   0   0   
0   0   0   0   0   1   0   0   0   0   0   0   
0   0   0   0   0   0   1   0   0   0   0   0   
0   0   0   0   0   0   0   1   0   0   0   0   
0   0   0   0   0   0   0   0   1   0   0   0   
0   0   0   0   0   0   0   0   0   1   0   0   
0   0   0   0   0   0   0   0   0   0   1   0   
0   0   0   0   0   0   0   0   0   0   0   1   
 Mcanonical: 
0  0.707073  0  0  0  0  0  0  0  0  0  0  
-0.707073  0  0  0  0  0  0  0  0  0  0  0  
0  0  0  0.707073  0  0  0  0  0  0  0  0  
0  0  -0.707073  0  0  0  0  0  0  0  0  0  
0  0  0  0  0  0.00687151  0  0  0  0  0  0  
0  0  0  0  -0.00687151  0  0  0  0  0  0  0  
0  0  0  0  0  0  0  0.00687151  0  0  0  0  
0  0  0  0  0  0  -0.00687151  0  0  0  0  0  
0  0  0  0  0  0  0  0  0  0.000784735  0  0  
0  0  0  0  0  0  0  0  -0.000784735  0  0  0  
0  0  0  0  0  0  0  0  0  0  0  0.000784735  
0  0  0  0  0  0  0  0  0  0  -0.000784735  0

另外测试本征值有 0 的反对称实矩阵

0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 -0 0 -0.0582782 0.0488136 0.146359 -0.00607454 0.141333 -0.0181291 0.130991 -0.00949629 0.16139
-0 -0.0 0.0582782 0 0.0352467 -0.00109344 0.1232 -0.0747338 0.148711 -0.242372 0.0586949 -0.142309
-0.0 -0.0 -0.0488136 -0.0352467 0 0.17432 0.165275 0.123782 0.257831 0 0.0425534 0.19093
-0.0 -0.0 -0.146359 0.00109344 -0.17432 0 -0.0727708 0.068801 0 0.257831 -0.098141 0.235799
-0.0 -0.0 0.00607454 -0.1232 -0.165275 0.0727708 0 0 0.068801 -0.123782 -0.151594 -0.00598481
-0.0 -0.0 -0.141333 0.0747338 -0.123782 -0.068801 -0 0 0.0727708 0.165275 -0.00598481 0.151594
-0.0 -0.0 0.0181291 -0.148711 -0.257831 -0 -0.068801 -0.0727708 0 -0.17432 -0.235799 -0.098141
-0.0  0.0 -0.130991 0.242372 -0 -0.257831 0.123782 -0.165275 0.17432 0 0.19093 -0.0425534
-0.0 -0.0 0.00949629 -0.0586949 -0.0425534 0.098141 0.151594 0.00598481 0.235799 -0.19093 0 0
-0.0 -0.0 -0.16139 0.142309 -0.19093 -0.235799 0.00598481 -0.151594 0.098141 0.0425534 -0 0

得到

 U: 
0  0  -0.395163  0.44073  -0.117531  -0.441911  -0.0208166  -0.124056  0.0612319  0.544158  0.214009  0.281114  
0  0  0  0.250065  0.0198452  -0.378119  0.244911  -0.395889  0.300962  -0.484856  0.237043  -0.442196  
0  0  -0.245674  -0.172487  -0.641937  0.0870772  0.146768  -0.022326  0.518727  -0.15687  -0.33935  0.243605  
0  0  0  -0.0782339  -0.39862  -0.335716  -0.428845  -0.264016  -0.472958  -0.00912163  -0.394059  -0.299347  
0  0  -0.552342  0.318859  0.166778  -0.0159522  -0.0533873  0.348302  -0.28037  -0.548004  -0.17742  0.174987  
0  0  0  0.503196  -0.131035  0.432106  -0.493068  0.169037  0.311836  0.101262  0.0183059  -0.404213  
0  0  -0.422101  -0.332344  0.452347  0.132859  -0.41161  -0.499227  0.248035  0.00468517  -0.0690757  0.0643759  
0  0  0  -0.108948  0.342077  -0.437731  0.0463857  0.432066  0.346534  0.211698  -0.490304  -0.292186  
0  0  -0.547941  -0.30591  -0.144001  0.193388  0.320104  0.132951  -0.185184  0.226694  0.229868  -0.537939  
0  0  0  0.369557  0.165698  0.328366  0.462741  -0.392692  -0.145453  0.193199  -0.545663  -0.0598417  
0  1  0  0  0  0  0  0  0  0  0  0  
1  0  0  0  0  0  0  0  0  0  0  0  
 U * U^+ : 
1   0   0   0   0   0   0   0   0   0   0   0   
0   1   0   0   0   0   0   0   0   0   0   0   
0   0   1   0   0   0   0   0   0   0   0   0   
0   0   0   1   0   0   0   0   0   0   0   0   
0   0   0   0   1   0   0   0   0   0   0   0   
0   0   0   0   0   1   0   0   0   0   0   0   
0   0   0   0   0   0   1   0   0   0   0   0   
0   0   0   0   0   0   0   1   0   0   0   0   
0   0   0   0   0   0   0   0   1   0   0   0   
0   0   0   0   0   0   0   0   0   1   0   0   
0   0   0   0   0   0   0   0   0   0   1   0   
0   0   0   0   0   0   0   0   0   0   0   1   
 Mcanonical: 
0  0.680657  0  0  0  0  0  0  0  0  0  0  
-0.680657  0  0  0  0  0  0  0  0  0  0  0  
0  0  0  0.553307  0  0  0  0  0  0  0  0  
0  0  -0.553307  0  0  0  0  0  0  0  0  0  
0  0  0  0  0  0.00614115  0  0  0  0  0  0  
0  0  0  0  -0.00614115  0  0  0  0  0  0  0  
0  0  0  0  0  0  0  0.00304889  0  0  0  0  
0  0  0  0  0  0  -0.00304889  0  0  0  0  0  
0  0  0  0  0  0  0  0  0  0.000424323  0  0  
0  0  0  0  0  0  0  0  -0.000424323  0  0  0  
0  0  0  0  0  0  0  0  0  0  0  0  
0  0  0  0  0  0  0  0  0  0  0  0

通过了测试。

4. 应用

4.1 原子核配对结构系数正则化

核子配对定义为，

\[\hat{P}^\dagger = \frac{1}{2} \sum_{\alpha \beta} p_{\alpha \beta} \hat{c}^\dagger_\alpha \hat{c}^\dagger_\beta, \]

其中 \(p\) 是对结构系数，是一个反对称矩阵。如果对单粒子基做一个幺正变换，

\[\hat{c}^\dagger_{\alpha'} = \sum_\alpha \hat{c}^\dagger_\alpha U_{\alpha' \alpha}, \]

其中 \(UU^\dagger = U^\dagger U = 1\),　则逆变换为

\[\hat{c}^\dagger_{\alpha} = \sum_\alpha \hat{c}^\dagger_{\alpha'} U^\dagger_{\alpha \alpha'}. \]

那么集体对产生算符 \(\hat{P}^\dagger\) 可以写作

\[\hat{P}^\dagger = \frac{1}{2} \sum_{\alpha\beta} p_{\alpha\beta} \sum_{\alpha'\beta'} U^\dagger_{\alpha \alpha'} U^\dagger_{\beta \beta'} \hat{c}^\dagger_{\alpha'} \hat{c}^\dagger_{\beta'} =\frac{1}{2} \sum_{\alpha'\beta'} ( \sum_{\alpha \beta} U^\dagger_{\alpha\alpha'} p_{\alpha\beta} U^\dagger_{\beta \beta'} )\hat{c}^\dagger_{\alpha'} \hat{c}^\dagger_{\beta'} = \frac{1}{2}\sum_{\alpha' \beta'} (U^* p U^\dagger)_{\alpha' \beta'} \hat{c}^\dagger_{\alpha'} \hat{c}^\dagger_{\beta'}. \]

换言之，在新的单粒子基下，新的对结构系数是 \(p' = U^* p U^\top\), 或 \(p = U^\top p' U^*\)。

根据前面的分析，对于给定实数矩阵 \(p\)，总可以找到实数正交阵 \(U\)，使得 \(p' = U p U^\top\) 是正则形式。
\(U[i][]\) 表示第 i 个新基矢在原基矢下的坐标。

4.2 从 Hartree-Fock 波函数得到 \(M\)-scheme 配对

如果体系有 \(2n\) 个同类核子，定义

\[\hat{P}^\dagger = \hat{P}^\dagger_{12} + \cdots + \hat{P}^\dagger_{2n-1, 2n}, \]

其中

\[\hat{P}^\dagger_{12} = \hat{c}^\dagger_{1'} \hat{c}^\dagger_{2'}, ~~~ \cdots, ~~~ \hat{P}^\dagger_{2n-1,2n} = \hat{c}^\dagger_{(2n-1)'} \hat{c}^\dagger_{(2n)'}, \]

\(\hat{c}^\dagger_{i'}\) 代表 Hartree-Fock 的单粒子基，在球形基下的展开为

\[\hat{c}^\dagger_{\alpha'} = \sum_\alpha \hat{c}^\dagger_\alpha U_{\alpha' \alpha}. \]

那么，\(\hat{P}^\dagger\) 的凝聚

\[(\hat{P}^\dagger)^n = ( \hat{P}^\dagger_{12} + \cdots + \hat{P}^\dagger_{2n-1, 2n} )^n = n! \hat{P}^\dagger_{12} \cdots\hat{P}^\dagger_{2n-1, 2n} | 0 \rangle, \]

是一个 Slater Determinant，即 Hartree-Fock 波函数。

所以，如果定义 canonical form 的集体对

\[p' = \left[(\begin{smallmatrix} 0 & 1 \\ -1 & 0 \end{smallmatrix}), ~ \cdots,~ (\begin{smallmatrix} 0 & 1 \\ -1 & 0 \end{smallmatrix}),~ 0, ~ \cdots,~ 0 \right]. \]

然后，根据 4.1 中的论述，这个集体对转回球形基，会得到 \(p = U^\top p' U^*\)，

\[\hat{P}^\dagger = \frac{1}{2} \sum_{\alpha \beta} p_{\alpha \beta} \hat{c}^\dagger_\alpha \hat{c}^\dagger_\beta. \]

posted on 2022-05-13 16:14 luyi07 阅读(315) 评论(1) 编辑收藏举报

刷新页面返回顶部

luyi07的笔记