反对称实矩阵正则化

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

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

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

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

H = i K,

$H = iK,$

则有

H^{†} = - i K^{†} = i K = H,

$H^\dagger = - i K^\dagger = i K = H,$

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

H = U Λ U^{†} ，

$H = U \Lambda U^\dagger，$

那么，

K = - i H = U (- i Λ) U^{†},

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

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

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

假设 $K$ 有特征值 $\lambda$ ，有

K \vec{z} = λ \vec{z},

$K \vec{z} = \lambda \vec{z},$

则有

(K {\vec{z}}^{*})^{†} = (K \vec{z})^{⊤} = λ {\vec{z}}^{⊤},

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

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

K {\vec{z}}^{*} = (λ {\vec{z}}^{⊤})^{†} = λ^{*} {\vec{z}}^{*},

$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}^{*}, \dots] d i a g [λ_{1}, λ_{1}^{*}, \dots, λ_{k}, λ_{k}^{*}, 0, \dots, 0] [{\vec{z}}_{1}, {\vec{z}}_{1}^{*}, \dots]^{†}

$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} = \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.

$\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{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} ⊥ \vec{y} .

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

另外

{{\vec{z}}_{k}, {\vec{z}}_{k}^{*}} ⊥ {{\vec{z}}_{l}, {\vec{z}}_{l}^{*}} \Rightarrow {{\vec{x}}_{k}, {\vec{y}}_{k}} ⊥ {{\vec{x}}_{l}, {\vec{y}}_{l}},

$\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\},$

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

{{\vec{x}}_{1}, {\vec{y}}_{1}, \dots, {\vec{x}}_{k}, {\vec{y}}_{k}, {\vec{w}}_{2 k + 1}, \dots, {\vec{w}}_{n}},

$\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}, \dots, {\vec{x}}_{k}, {\vec{x}}_{y}, {\vec{w}}_{2 k + 1}, \dots, {\vec{w}}_{n}] = [{\vec{x}}_{1}, {\vec{y}}_{1}, \dots, {\vec{x}}_{k}, {\vec{x}}_{y}, {\vec{w}}_{2 k + 1}, \dots, {\vec{w}}_{n}] [[\begin{matrix} 0 & - i \sqrt{2} λ_{1} \\ i \sqrt{2} λ_{1} & 0 \end{matrix}], \dots, [\begin{matrix} 0 & - i \sqrt{2} λ_{k} \\ i \sqrt{2} λ_{k} & 0 \end{matrix}], 0, \dots, 0] .

$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}, \dots, {\vec{x}}_{k}, {\vec{x}}_{y}, {\vec{w}}_{2 k + 1}, \dots, {\vec{w}}_{n}]^{⊤} K [{\vec{x}}_{1}, {\vec{y}}_{1}, \dots, {\vec{x}}_{k}, {\vec{x}}_{y}, {\vec{w}}_{2 k + 1}, \dots, {\vec{w}}_{n}] = [[\begin{matrix} 0 & - i \sqrt{2} λ_{1} \\ i \sqrt{2} λ_{1} & 0 \end{matrix}], \dots, [\begin{matrix} 0 & - i \sqrt{2} λ_{k} \\ i \sqrt{2} λ_{k} & 0 \end{matrix}], 0, \dots, 0] .

$[ \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} = λ {\vec{z}}_{1}, K {\vec{z}}_{2} = λ {\vec{z}}_{2}, {\vec{z}}_{1} ⊥ {\vec{z}}_{2}, | {\vec{z}}_{1} | = | {\vec{z}}_{2} | = 1; K {\vec{z}}_{3} = - λ {\vec{z}}_{3}, K {\vec{z}}_{4} = - λ {\vec{z}}_{3}, {\vec{z}}_{3} ⊥ {\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}] d i a g {λ, λ, - λ, - λ} . {\vec{z}}_{i} = ({\vec{x}}_{i} + i {\vec{y}}_{i}) / \sqrt{2} .

$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} = μ {\vec{z}}_{i}, \Rightarrow K {\vec{z}}_{i}^{*} = - μ {\vec{z}}_{i}^{*}, \Rightarrow {\vec{z}}_{i} ⊥ {\vec{z}}_{i}^{*}

$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} + 2 i {\vec{x}}_{i} \cdot {\vec{y}}_{i}) = 0, \Rightarrow | {\vec{x}}_{i} | = | {\vec{y}}_{i} |, {\vec{x}}_{i} ⊥ {\vec{y}}_{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{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} ⊥ {\vec{y}}_{i} .

$|\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}] [[\begin{matrix} 0 & - i \sqrt{2} λ_{1} \\ i \sqrt{2} λ_{1} & 0 \end{matrix}], [\begin{matrix} 0 & - i \sqrt{2} λ_{2} \\ i \sqrt{2} λ_{2} & 0 \end{matrix}]] .

$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^{†} = H .

$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} .

$\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}}^{†} = \frac{1}{2} \sum_{α β} p_{α β} {\hat{c}}_{α}^{†} {\hat{c}}_{β}^{†},

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

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

{\hat{c}}_{α^{'}}^{†} = \sum_{α} {\hat{c}}_{α}^{†} U_{α^{'} α},

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

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

{\hat{c}}_{α}^{†} = \sum_{α} {\hat{c}}_{α^{'}}^{†} U_{α α^{'}}^{†} .

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

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

{\hat{P}}^{†} = \frac{1}{2} \sum_{α β} p_{α β} \sum_{α^{'} β^{'}} U_{α α^{'}}^{†} U_{β β^{'}}^{†} {\hat{c}}_{α^{'}}^{†} {\hat{c}}_{β^{'}}^{†} = \frac{1}{2} \sum_{α^{'} β^{'}} (\sum_{α β} U_{α α^{'}}^{†} p_{α β} U_{β β^{'}}^{†}) {\hat{c}}_{α^{'}}^{†} {\hat{c}}_{β^{'}}^{†} = \frac{1}{2} \sum_{α^{'} β^{'}} (U^{*} p U^{†})_{α^{'} β^{'}} {\hat{c}}_{α^{'}}^{†} {\hat{c}}_{β^{'}}^{†} .

$\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}}^{†} = {\hat{P}}_{12}^{†} + \dots + {\hat{P}}_{2 n - 1, 2 n}^{†},

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

其中

{\hat{P}}_{12}^{†} = {\hat{c}}_{1^{'}}^{†} {\hat{c}}_{2^{'}}^{†}, \dots, {\hat{P}}_{2 n - 1, 2 n}^{†} = {\hat{c}}_{(2 n - 1)^{'}}^{†} {\hat{c}}_{(2 n)^{'}}^{†},

$\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}}_{α^{'}}^{†} = \sum_{α} {\hat{c}}_{α}^{†} U_{α^{'} α} .

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

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

({\hat{P}}^{†})^{n} = ({\hat{P}}_{12}^{†} + \dots + {\hat{P}}_{2 n - 1, 2 n}^{†})^{n} = n! {\hat{P}}_{12}^{†} \dots {\hat{P}}_{2 n - 1, 2 n}^{†} | 0 ⟩,

$(\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^{'} = [(\begin{matrix} 0 & 1 \\ - 1 & 0 \end{matrix}), \dots, (\begin{matrix} 0 & 1 \\ - 1 & 0 \end{matrix}), 0, \dots, 0] .

$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}}^{†} = \frac{1}{2} \sum_{α β} p_{α β} {\hat{c}}_{α}^{†} {\hat{c}}_{β}^{†} .

$\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 阅读(329) 评论(1) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 群论基础(3)：有限群表示论

· General Seniority 学习笔记(1): BCS

· NumPy-1-26-中文文档-十一-

· 你疑似有点太线代化了

· Note -「A. Algebra 24 Aut.」“还有一束日光正在为你送达”

阅读排行：
· TypeScript + Deepseek 打造卜卦网站：技术与玄学的结合
· Manus的开源复刻OpenManus初探
· 写一个简单的SQL生成工具
· AI 智能体引爆开源社区「GitHub 热点速览」
· C#/.NET/.NET Core技术前沿周刊 | 第 29 期（2025年3.1-3.9）

历史上的今天：
2020-05-13 输入输出重定向

luyi07的笔记

反对称实矩阵正则化

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

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

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

2. 反对称实矩阵正则化

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

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

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

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

3. 数值验证

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

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

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

4. 应用

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

4.2 从 Hartree-Fock 波函数得到 $M$ -scheme 配对

导航

公告

搜索

常用链接

随笔分类

随笔档案

阅读排行榜

评论排行榜

推荐排行榜

最新评论

luyi07的笔记

反对称实矩阵正则化

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

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

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

2. 反对称实矩阵正则化

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

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

2.2.1 →xi,→yix→i,y→i\vec{x}_i, \vec{y}_i 这两个矢量正交归一

2.2.2 →x1,→y1,→x2,→y2x→1,y→1,x→2,y→2\vec{x}_1, \vec{y}_1, \vec{x}_2, \vec{y}_2正交

3. 数值验证

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

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

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

4. 应用

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

4.2 从 Hartree-Fock 波函数得到 MMM-scheme 配对

导航

公告

搜索

常用链接

随笔分类

随笔档案

阅读排行榜

评论排行榜

推荐排行榜

最新评论

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

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

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

4.2 从 Hartree-Fock 波函数得到 $M$ -scheme 配对