这是上次一个小文献笔记()里一个定理的实践。

1. 实数反对称矩阵 \(M\)

所有矩阵元为实数，并且有反对称性 \(M^\top = - M\)。

2. 反对称阵的正则形式

如果反对称矩阵 \(M\) 是块对角的，且每块都是 2x2 的矩阵 \([\begin{smallmatrix} 0 & -\mu \\ \mu & 0 \end{smallmatrix}]\)，则称这种形式为正则形式。

3. 二度简并实数反对称矩阵正则化

3.1 理论来源

文献笔记()中有记录这个定理：任意反对称复数矩阵 \(M\) 都可以变换为正则形式 \(M = U F U^\top\)，其中 \(U\) 是人为构造的一个幺正矩阵，\(F\) 是正则矩阵。
实数反对称矩阵是复数反对称矩阵的特例，这里我使用文献中证明过程的思想，记下二度简并实数反对称矩阵正则化计算过程、代码、效果。

3.2 计算公式

首先，拿着实数反对称阵 \(M\)，构造 \(H = M^\top M\)，容易证明 \(H\) 是厄米半正定的，所以可以通过实数正交相似变换，对角化为实数对角阵 \(H_d\)：

另外构造矩阵 \(M_1 = V M V^\top\)，则有

因为 \(M, M^\top\) 可交换，所以有 \(M_1 H_d = H_d M_1\)，而 \(H_d\) 是对角阵，所以有

所以当 \(d_i \neq d_j\) 时，有 \((M_1)_{ij} = 0\)，\(M_1\) 是分块矩阵。
\(M\) 为二度简并　即 \(H_d\) 的本征值为一些２重简并的值，那就意味着 \(M_1\) 的对角块都是 2x2，而 \(M_1\) 按定义是反对称阵，这就说明了 \(M_1\) 是正则形式。
二度简并这个假定并不特别过分，与之相对的是 3 度或者更高度的简并，想必是非常少见的。而二度简并是这个形式自带的特征，如果没有更高度的简并，就一定有二度简并，我不知道怎么证明（所以也不是100%确定），但实践了几个例子，似乎都是如此。

3.3 c++代码

下面这个函数对于给定的２度简并的实数反对称矩阵 \(M\)，进行正则化，得到 \(M1\)，相似变换信息储存在 \(V\) 中。

// RealSkewCanonical: construct a orthogonal similar transformation, turn a given real skew matrix into canonical form
// input: n, M; output: V, M1
// H = M^\top M, H V^\top = V^\top \Lambda
// M1 = V M V^\top is a canonical matrix
void RealSkewCanonical( int n, double * M, double * V, double * M1 ){
        double * H = new double [ n*n ];
        cblas_dgemm( CblasRowMajor, CblasTrans, CblasNoTrans, n, n, n, 1, M, n, M, n, 0, H, n ); // H = M^\dagger M
        double *e = new double [n]; LapackDsyev( H, V, e, n ); // diagonalize H, get V
        double * C = new double [ n*n ];
        cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1, V, n, M, n, 0, C, n ); //C = VM
        cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasTrans, n, n, n, 1, C, n, V, n, 0, M1, n ); // M1 = CV^\top = VMV^\top
        delete [] C; delete [] e; delete [] H;
}

3.4 运行结果

然后写了个主函数，进行测试，干脆一起贴在下面吧

#include
using namespace std;
#include
#include"mkl.h"
#include
#include

#include"dsyev.h"

void printmtx( int n, double * A ){
        for(int i=0;i>A[i];
        fin.close();
}

void randmtx( int n, double * A ){
        for(int i=0;i

0  1  2  5
-1 0  -9 10
-2 9  0  7
-5 -10 -7 0

M1:
5.55112e-17, 3.69536, 1.72085e-15, -2.10942e-15, 
-3.69536, -5.55112e-17, 2.33147e-15, 6.66134e-16, 
0, -2.44249e-15, 8.88178e-16, -15.6954, 
3.55271e-15, -2.22045e-16, 15.6954, -8.88178e-16,