模意义下矩阵对角化的暴力尝试
矩阵对角化的一些简单尝试
给定 \(n\) 阶方阵 \(A\)，求两个矩阵 \(D,P\)，满足 \(A=PDP^{-1}\)。题中给出两个具体方阵，分别为
\[\begin{bmatrix} 1 & 2 & 4 \\ 8 & 6 & 6 \\ 4 & 2 & 3 \end{bmatrix}\quad \begin{bmatrix} 1 & 4 & 8 & 9 \\ 2 & 6 & 7 & 4 \\ 8 & 7 & 4 & 4 \\ 9 & 2 & 9 & 2 \end{bmatrix} \]
如果仅需要求出一组合理的数值解，只需利用已有的工具（如 Matlab，python）或在线计算器（如 Wolfram Alpha）即可。以下是由 AI 生成的 python 代码所计算出的数值结果。
结果经过验算，虚数部分剩余 \(\leq 10^{-15}i\) ，可以接受。
当然这么长一串数值结果显然不适合人眼来阅读，并且无法体现认真做了此项作业，因此我们需要找到一种更简单的方式来表示计算结果，也可以在此过程中理解矩阵对角化的计算过程。
作为线性代数的菜鸟，目前尚未系统阅读任何线代教材，我们仅从 \(A=PDP^{-1}\) 入手来人类本能地推出计算方法。
首先原式同时右乘 \(P\)，得
记 \(D=\text{diag}\{\lambda_1,\cdots,\lambda_n\}\)，\(P=[\vec{p_1},\cdots,\vec{p_n}]\)，（\(\vec{p_i}\) 为列阵）可以得
由于我们的目标，\(\vec{p_i}\) 不能为 \(\vec{0}\)。则我们需要对 \((A-\lambda_iI)\) 的性质进行探讨，不妨设 \(B=A-\lambda_iI\)，即 \(B\vec{p_i}=\vec{0}\)。
假设 \(B\) 矩阵可逆，则
又要求 \(\vec{p_i}\neq \vec{0}\)，所以原假设矛盾。因此 \(B\) 不可逆，等价于 \(\det(B)=0\)，即
因此我们解 \(n\) 次方程 \(\det(A-\lambda I)=0\)，若得到 \(n\) 个不等根，即可通过方程 \((A-\lambda_iI)\vec{p_i}=\vec{0}\) 解得 \(n\) 个不线性相关的 \(\vec{p}\) 构成 \(P\) 矩阵。
同时注意到 \(A-\lambda I=0\) 并不满秩，因此解上述方程时需将若干个未知数用一个钦定的未知数表示。
通过上面完全不严谨的表述，我们可以概括流程如下：
- 解 \(n\) 次方程 \(\det(A-\lambda I)=0\) 的共 \(n\) 个根。
- 将根 \(\lambda_i\) 带入，解方程 \((A-\lambda_iI)\vec{p_i}=\vec{0}\) 的一组合法解 \(\vec{p_i}\)。
- 通过 \(D=\text{diag}\{\lambda_1,\cdots,\lambda_n\}\)，\(P=[\vec{p_1},\cdots,\vec{p_n}]\)，计算 \(P^{-1}\)，并验算。
问题的关键在求解 \(n\) 个特征值 \(\lambda\)。
对于题目 \(3\times 3\) 的方阵，其特征多项式方程为
尚且可以得到 \(3\) 个复杂程度还能接受的解，并将 \(P,D\) 矩阵用分数与根式的形式表示。
而题目中 \(4\) 阶多项式的方程拥有 \(2\) 个实数解 \(2\) 个虚数解，同时每个解都只能表示成较为复杂的分数形式，因此解得的分数形式过于复杂，只能通过实根隔离和牛顿迭代法求解数值解（实际上我目前仍然也不是很会）。当然也是可以使用一元四次方程求根公式的，但是这样做我们需要首先求出特征多项式（这个过程也较为繁琐），再套一系列复杂的公式。
但是，如果我们将其看成一个不具备实际意义的数字游戏，这个对角化的问题却能找到若干个简单而直观的解。
在对结果取模 \(p\) 的意义下，加减乘运算均不用改变，除法运算可以转化为乘一个数 \(a\) 的逆元 \(x\)，即 $ax\equiv 1\pmod p $ 的解。
而对于根式，\(\sqrt a\) 的值相当于找到 \(x^2\equiv a\pmod p\) 的解，特别的，虚数单位 \(i=\sqrt{-1}\) 相当于找到 \(x^2\equiv p-1\pmod p\) 的解。
当然并不是所有的模数我们都能找到合适的根式解，例如对于 \(p=11\)，我们就无法找到对应虚数单位 \(i\) 的解。对于寻找根式解的问题，已经有一系列成熟、高效且普及的算法，但是我们的目标是简化问题而不是复杂化问题。我们只需要找到一个模数 \(p\)，使得在其剩余系中特征值方程有 \(n\) 个不同的解即可。
这个解当然也可以通过求根方式求解，但是我们将问题转化到模意义下的最大好处就是所有数变成了 \([0,p-1]\) 内的整数，我们直接枚举这 \(p\) 个整数带入特征值方程看看是否为 \(0\) 即可。
至此，我们可以将解方程拆解成以下 \(3\) 步
- 枚举质数 \(p\)。
- 枚举 \(0\sim p-1\) 的整数 $\lambda $，求 \(\det(A-\lambda I)\pmod p\) 是否为 \(0\)。
- 如果 $\lambda $ 的解有 \(n\) 个，则可以选择 \(p\) 作为模数，进行剩余步骤。
虽然我们必须承认拥有 \(4\) 个实根的模数 \(p\) 分布就已经非常稀疏了（\(10000\) 以内的素数共 \(1229\) 个，其中只有 \(59\) 个满足），但是幸运的是只需要取 \(p=11\) 即可满足题中 \(3\) 阶行列式的要求，取 \(p=19\) 即可满足 \(4\) 阶行列式的要求。以下是得出的答案
于是我们就在一行写出了答案。
算法时间复杂度 \(O(\frac{p^2 n^3}{\ln p})\)，其中 \(p\) 为最小的合法质数，\(\ln p\) 来自于 \(1\sim n\) 质数分布约为 \(\frac{1}{\ln n}\)，\(n^3\) 来自行列式求值。
后记
正文部分没有加上小标题，因为我每次最后都会发现我写的小标题一般都挺差的，干脆用分割线隔开算了，如果显得非常不专业抱歉。
虽然最后也是勉强把结果算出来了，但是明显整个过程的表述都不够专业而且由于基础太差，有些小问题想了好久才想明白，然后有些潜在的问题直接忽略了。
本来还想利用写好的程序多写些功能并验证一些东西，但是由于就这个现在看来非常简陋的东西已经耗费了太多精力，加上还有若干门课作业没写，于是就搁置了，不如等以后数学基础更牢再做。
最重要的是废话很多，因为知道的很少，只能重复一些简单的东西。
代码
python
python 计算数值解，并验证
import numpy as np
# 定义一个矩阵
A = np.array([[1, 4, 8, 9],
[2, 6, 7, 4],
[8, 7, 4, 4],
[9, 2, 9, 2]])
# 计算特征值和特征向量
eigenvalues, eigenvectors = np.linalg.eig(A)
# 对角矩阵
D = np.diag(eigenvalues)
# 矩阵P（特征向量）
P = eigenvectors
# 验证对角化：P * D * P_inv 应该等于原矩阵 A
P_inv = np.linalg.inv(P)
A_diagonalized = P @ D @ P_inv
# 此处为了方便查看结果对 AI 代码进行了小幅修改
print("A = " , A )
print("P = " , P )
print("P_inv = ", P_inv )
print("D = ", D )
print("A_diagonalized=", A_diagonalized)
c++
又臭又长而且非主流
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef vector<vector<ll>> matrix;
const int N=10;
inline ll inv(ll b,ll p){//费马小定理求除法逆元
ll t=p-2ll,s=1ll;
for(;t;t>>=1,b=b*b%p)
if(t&1ll)s=s*b%p;
return s;
}
vector<ll>prime;
void getprime(int R){//线性筛法求质数
vector<bool>isp(R+1,1);
for(ll x=2;x<=R;++x){
if(isp[x])prime.push_back(x);
for(int i=0;i<(int)prime.size()&&prime[i]*x<=R;++i){
isp[prime[i]*x]=0;
if(x%prime[i]==0)break;
}
}
}
int n;
void input(matrix &a){
for(int i=0;i<n;++i){
for(int j=0;j<n;++j)
scanf("%lld",&a[i][j]);
}
}
void output(matrix a){
for(int i=0;i<n;++i){
for(int j=0;j<n;++j)
printf("%lld ",a[i][j]);
puts("");
}
}
matrix times(matrix a,matrix b,ll p){//矩阵乘法
matrix c(n,vector<ll>(n,0));
for(int i=0;i<n;++i){
for(int j=0;j<n;++j){
for(int k=0;k<n;++k)
c[i][j]=(c[i][j]+a[i][k]*b[k][j]%p)%p;
}
}
return c;
}
ll det(matrix a,ll p){//高斯消元求行列式，利用辗转相除实现
ll f=1,val=1ll;
for(int i=0;i<n;++i){
for(int j=i+1;j<n;++j){
while(a[i][i]){
ll tim=a[j][i]/a[i][i];
for(int k=i;k<n;++k)
a[j][k]=(a[j][k]-tim*a[i][k]%p+p)%p;
swap(a[i],a[j]);f=-f;
}
swap(a[i],a[j]);f=-f;
}
}
for(int i=0;i<n;++i)val=val*a[i][i]%p;
return (val*f+p)%p;
}
matrix matrix_inv(matrix a,ll p){//矩阵求逆
matrix b(n,vector<ll>(2*n,0));
for(int i=0;i<n;++i){
for(int j=0;j<n;++j)
b[i][j]=a[i][j];
b[i][i+n]=1;
}
for(int i=0;i<n;++i){
if(!b[i][i]){
for(int k=i+1;k<n;++k){
if(b[k][i]){
swap(b[i],b[k]);
break;
}
}
}
ll d=inv(b[i][i],p);
for(int k=0;k<n*2;++k)
b[i][k]=b[i][k]*d%p;
for(int j=0;j<n;++j){
if(j==i||b[j][i]==0)continue;
ll bas=b[j][i];
for(int k=0;k<n*2;++k)
b[j][k]=(b[j][k]-bas*b[i][k]%p+p)%p;
}
}
matrix aa(n,vector<ll>(n));
for(int i=0;i<n;++i){
for(int j=0;j<n;++j)
aa[i][j]=b[i][j+n];
}
return aa;
}
vector<ll> gauss(matrix a,ll p){//求解特殊方程一组解
for(int i=0;i<n-1;++i){
if(!a[i][i]){
for(int k=i+1;k<n;++k){
if(a[k][i]){
swap(a[i],a[k]);
break;
}
}
}
ll d=inv(a[i][i],p);
for(int k=0;k<n;++k)
a[i][k]=a[i][k]*d%p;
for(int j=0;j<n;++j){
if(i==j||a[j][i]==0)continue;
ll bas=a[j][i];
for(int k=0;k<n;++k)
a[j][k]=(a[j][k]-bas*a[i][k]%p+p)%p;
}
}
vector<ll>vec(n);
vec[n-1]=1;
for(int i=0;i<n-1;++i)vec[i]=(p-a[i][n-1])%p;
return vec;
}
bool calc(matrix A,ll p){//主体步骤
vector<ll>ans;
matrix B=A;
for(ll x=0;x<p;++x){
for(int i=0;i<n;++i)
B[i][i]=(A[i][i]-x+p)%p;
if(det(B,p)==0)
ans.push_back(x);
}
if((int)ans.size()!=n)return 0;
matrix D(n,vector<ll>(n,0));
for(int i=0;i<n;++i)D[i][i]=ans[i];
matrix P(n,vector<ll>(n,0));
for(int j=0;j<n;++j){
B=A;
for(int i=0;i<n;++i)
B[i][i]=(A[i][i]-ans[j]+p)%p;
vector<ll> vec=gauss(B,p);//vec 是一个列阵
for(int i=0;i<n;++i)
P[i][j]=vec[i];
}
matrix P_1=matrix_inv(P,p);
puts("P: ");output(P);puts("");
puts("D: ");output(D);puts("");
puts("P_inv: ");output(P_1);puts("");
puts("P*D*P_inv: ");output(times(times(P,D,p),P_1,p));//结果验算
return 1;
}
int main(){
getprime(10000000);//预处理 10^7 内质数
scanf("%d",&n);
matrix A(n,vector<ll>(n));
input(A);
for(auto p:prime){
if(p==2)continue;
if(calc(A,p)){
printf("prime: %lld\n",p);
return 0;
}
}
puts("Can't find a prime smaller than 1e7 to solve it.");
return 0;
}
最后也是自己又编了几组更大的数据数据跑了一下
5
3 1 6 1 4
6 4 2 9 3
2 7 3 3 1
9 2 1 8 4
2 4 6 7 9
P:
7 94 21 46 63
19 6 2 72 87
8 18 30 38 90
77 103 102 33 95
1 1 1 1 1
D:
8 0 0 0 0
0 33 0 0 0
0 0 49 0 0
0 0 0 57 0
0 0 0 0 106
P_inv:
11 86 51 30 92
25 76 39 35 7
4 7 36 111 44
26 64 100 63 70
47 106 0 100 14
P*D*P_inv:
3 1 6 1 4
6 4 2 9 3
2 7 3 3 1
9 2 1 8 4
2 4 6 7 9
prime: 113
\(n=8\) 跑了 \(1711\) 秒
8
23 13 62 13 26 12 15 13
11 19 13 29 24 3 71 13
9 13 14 4 32 12 42 23
12 23 2 25 23 3 9 8
18 23 53 34 63 69 43 64
12 14 32 86 39 37 43 61
31 67 28 73 25 73 11 86
6 72 43 36 42 63 23 97
P:
21389 45759 41022 27771 26780 35083 56499 40329
26569 15431 38776 8597 58345 6260 18825 5061
3579 37585 52686 45570 17275 30444 16549 23301
38407 18293 11709 4324 21807 17577 12236 49964
6859 26625 56730 35374 47152 47 51604 12834
33773 60366 61065 57465 13940 24026 33455 52638
27214 1285 5889 38397 3384 32068 30122 7906
1 1 1 1 1 1 1 1
D:
8180 0 0 0 0 0 0 0
0 9923 0 0 0 0 0 0
0 0 22793 0 0 0 0 0
0 0 0 31933 0 0 0 0
0 0 0 0 32037 0 0 0
0 0 0 0 0 42180 0 0
0 0 0 0 0 0 47564 0
0 0 0 0 0 0 0 62435
P_inv:
22878 6150 53285 11838 55517 63229 60687 19807
45321 4277 48554 43254 17533 9346 21640 23623
20819 10559 6819 57298 33432 54726 13292 18163
56227 22513 702 22425 58520 47661 7730 18005
60314 24075 19115 38796 46449 2748 27268 6343
36424 55554 15333 59643 48888 37383 5997 37911
54821 2665 26989 31575 44138 39478 56324 48922
24141 2585 21770 56116 16468 2185 63818 19794
P*D*P_inv:
23 13 62 13 26 12 15 13
11 19 13 29 24 3 71 13
9 13 14 4 32 12 42 23
12 23 2 25 23 3 9 8
18 23 53 34 63 69 43 64
12 14 32 86 39 37 43 61
31 67 28 73 25 73 11 86
6 72 43 36 42 63 23 97
prime: 64189
通过理论时间复杂度计算需要 \(190,572,292,114\) 单位次计算。根据每秒运行 \(10^8\) 次的经验，与实测数据较为吻合。
浙公网安备 33010602011771号