這是上次一個小文獻筆記(https://www.cnblogs.com/luyi07/p/15442971.html)裡一個定理的實踐。

1. 實數反對稱矩陣 \(M\)

所有矩陣元為實數，並且有反對稱性 \(M^\top = - M\)。

2. 反對稱陣的正則形式

如果反對稱矩陣 \(M\) 是塊對角的，且每塊都是 2x2 的矩陣 \([\begin{smallmatrix} 0 & -\mu \\ \mu & 0 \end{smallmatrix}]\)，則稱這種形式為正則形式。

3. 二度簡併實數反對稱矩陣正則化

3.1 理論來源

文獻筆記(https://www.cnblogs.com/luyi07/p/15442971.html

)中有記錄這個定理：任意反對稱複數矩陣 \(M\) 都可以變換為正則形式 \(M = U F U^\top\)，其中 \(U\) 是人為構造的一個么正矩陣，\(F\) 是正則矩陣。
實數反對稱矩陣是複數反對稱矩陣的特例，這裡我使用文獻中證明過程的思想，記下二度簡併實數反對稱矩陣正則化計算過程、程式碼、效果。

3.2 計算公式

首先，拿著實數反對稱陣 \(M\)，構造 \(H = M^\top M\)，容易證明 \(H\) 是厄米半正定的，所以可以通過實數正交相似變換，對角化為實數對角陣 \(H_d\)：

\[H V^\top = V^\top H_d, ~~ V H V^\top = H_d. \]

另外構造矩陣 \(M_1 = V M V^\top\)

，則有

\[M_1 H_d = VMV^\top V M^\top M V^\top = V M M^\top M V^\top, \\ H_d M_1 = V M^\top M V^\top VMV^\top = V M^\top M M V^\top, \]

因為 \(M, M^\top\) 可交換，所以有 \(M_1 H_d = H_d M_1\)，而 \(H_d\) 是對角陣，所以有

\[(M_1 H_d)_{ij} = (M_1)_{ij} d_j = (H_d M_1)_{ij} = d_i (M_1)_{ij}, \]

所以當 \(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<iostream>
using namespace std;
#include<cmath>
#include"mkl.h"
#include<complex>
#include<fstream>

#include"dsyev.h"

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

void readmtx( int n, double * A, string file ){
        ifstream fin(file);
        for(int i=0;i<n*n;i++) fin>>A[i];
        fin.close();
}

void randmtx( int n, double * A ){
        for(int i=0;i<n*n;i++) A[i] = ((double)rand())/RAND_MAX;
}

// 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;
}

int main(){

        int n = 4;
        double * M = new double [ n*n ];
        readmtx( n, M, "M.txt" );//M is skew-symmetric
        double * V = new double [ n*n ];
        double * M1 = new double [ n*n ];
        RealSkewCanonical( n, M, V, M1 );
        cout<<"M1:"<<endl; printmtx( n, M1 );
        delete [] M; delete [] V; delete [] M1;
        return 0;
}

再自己編個檔案 M.txt：

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, 

可以看到，變成了正則形式。