【學習筆記】求矩陣的特征多項式


先膜拜一波神仙yww

給定一個矩陣(沒有任何特殊性質),如何求它的特征多項式?

算法一

直接把\(\lambda\)代入\((n+1)\)個點值,求完行列式之后插值即可。
時間復雜度\(O(n^4)\)

算法二

下面介紹一個更快的做法。
定義 對於矩陣\(\bm A,\bm B\), 若存在可逆矩陣\(\bm\Phi\)滿足\(\bm A=\bm\Phi^{-1}\bm B\bm\Phi\), 則稱\(\bm A\)\(\bm B\)相似矩陣\(\bm\Phi\)稱為橋接矩陣
定理 相似矩陣的特征多項式相同。
證明: 設\(\bm A\)\(\bm B\)為相似矩陣,\(f(\lambda)\)\(g(\lambda)\)分別為其特征多項式,則\(f(\lambda)=|\lambda \bm I-\bm A|=|\lambda \bm \Phi^{-1}\bm\Phi-\bm\Phi^{-1}\bm B\bm\Phi|=|\bm \Phi^{-1}(\lambda\bm I-\bm B)\bm\Phi|=|\bm\Phi^{-1}||\lambda \bm I-\bm B| |\bm\Phi|=|\bm\Phi^{-1}\bm\Phi|g(\lambda)=g(\lambda)\). 證畢。
於是我們可以考慮把沒有特殊性質的矩陣化成有特殊性質的與之相似的矩陣來加速運算。
我們知道,對矩陣進行初等行(列)變換相當於將其左(右)乘一個可逆的初等矩陣,那么我們根據相似矩陣的定義,對\(\bm A\)進行初等行變換左乘初等矩陣的同時給它右乘同一初等矩陣的逆,不就構造出相似矩陣了嗎?
簡單推導可知,若左乘某一初等矩陣等價於把原矩陣第\(i\)行乘以\(x\)加到第\(j\)行上,則右乘該初等矩陣的逆矩陣等價於把原矩陣第\(j\)行乘以\(-x\)加到第\(i\)行上。

現在的問題是,我們對這個沒有特殊性質的矩陣左乘初等矩陣的同時右乘它的逆,能把它變成什么有特殊性質的矩陣呢?
答案是,上海森堡矩陣。
定義 若一\(n\times n\)矩陣滿足對於任意\(1\le j<j+1<i\le n\)\((i,j)\), 矩陣第\(i\)行第\(j\)列值均為\(0\),則其為上海森堡矩陣
至於為什么是這個答案,首先變成上三角矩陣顯然不現實,因為用第\(i\)行來消第\(j\)行 (\(j>i\))的第\(i\)列時對第\(i\)行和第\(j\)行進行了初等行變換,那這對應的初等列變換勢必會影響第\(j\)列和第\(i\)列。那么我們可以考慮在消第\(i\)列的時候,用第\((i+1)\)行來消第\(j\)行 (\(j>i+1\))的第\(i\)列,這樣進行的初等行變換是第\((i+1)\)行和第\(j\)行之間的,影響到的是第\((i+1)\)列和第\(j\)列,而不會影響第\(i\)列了。

現在問題轉化為,如何快速求上海森堡矩陣的特征多項式?
這就非常簡單了,最后一行只有兩列有值,那么枚舉選的是哪一列,如果選第\(n\)列,那么直接就是\((\lambda-a_{n,n})\)乘以左上角\((n-1)\times (n-1)\)的子矩陣的特征多項式;若選第\((n-1)\)列,那么第\((n-1)\)行選\((n-2)\)列,\((n-2)\)行選\((n-3)\)列……一直到某一行為止。假設這一行是第\((j+1)\)行選第\(j\)列 (\(0<j<n\)),那么第\(j\)行一定選第\(n\)列,前\((j-1)\)行選的就是前\((j-1)\)列了,而\(a_{n-1,n-2},a_{n-2,n-3},...,a_{j+1,j},a_{j,n}\)都是常數(都不在對角線上),所以直接遞推即可。

時間復雜度\(O(n^3)\).

代碼

(每一項對\(998244353\)取模, \(n\le 500\))

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cassert>
#include<iostream>
#define llong long long
using namespace std;

inline int read()
{
	int x=0; bool f=1; char c=getchar();
	for(;!isdigit(c);c=getchar()) if(c=='-') f=0;
	for(; isdigit(c);c=getchar()) x=(x<<3)+(x<<1)+(c^'0');
	if(f) return x;
	return -x;
}

const int N = 500;
const int P = 998244353;
llong a[N+3][N+3];
llong c[N+3][N+3];
int n;

llong quickpow(llong x,llong y)
{
	llong cur = x,ret = 1ll;
	for(int i=0; y; i++)
	{
		if(y&(1ll<<i)) {y-=(1ll<<i); ret = ret*cur%P;}
		cur = cur*cur%P;
	}
	return ret;
}
llong mulinv(llong x) {return quickpow(x,P-2);}

void gauss()
{
	for(int i=1; i<=n; i++)
	{
		if(a[i+1][i]==0)
		{
			bool found = false;
			for(int j=i+2; j<=n; j++)
			{
				if(a[j][i]!=0)
				{
					for(int k=i; k<=n; k++) swap(a[i+1][k],a[j][k]);
					for(int k=1; k<=n; k++) swap(a[k][i+1],a[k][j]);
					found = false; break;
				}
			}
			if(found) {continue;}
		}
		for(int j=i+2; j<=n; j++)
		{
			llong coe = P-a[j][i]*mulinv(a[i+1][i])%P;
			for(int k=i; k<=n; k++) a[j][k] = (a[j][k]+coe*a[i+1][k])%P;
			for(int k=1; k<=n; k++) a[k][i+1] = (a[k][i+1]-coe*a[k][j]%P+P)%P;
		}
	}
}

void charpoly()
{
	c[0][0] = 1ll;
	for(int i=1; i<=n; i++)
	{
		for(int j=0; j<=i; j++)
		{
			c[i][j] = (c[i-1][j-1]-a[i][i]*c[i-1][j]%P+P)%P;
		}
		llong coe = P-1,cur = P-a[i][i-1];
		for(int j=i-2; j>=0; j--)
		{
			llong tmp = cur*(P-a[j+1][i])%P;
			tmp = coe*tmp%P;
			for(int k=0; k<=j; k++)
			{
				c[i][k] = (c[i][k]+c[j][k]*tmp)%P;
			}
			cur = cur*(P-a[j+1][j])%P;
			coe = P-coe;
		}
		for(int k=0; k<=i; k++) c[i][k] %= P;
//		printf("%d: ",i); for(int j=0; j<=i; j++) printf("%lld ",c[i][j]); puts("");
	}
}

int main()
{
	scanf("%lld",&n);
	for(int i=1; i<=n; i++) for(int j=1; j<=n; j++) scanf("%lld",&a[i][j]);
	gauss();
//	for(int i=1; i<=n; i++) {for(int j=1; j<=n; j++) printf("%lld ",a[i][j]); puts("");}
	charpoly();
	for(int i=0; i<=n; i++) printf("%lld ",c[n][i]); puts("");
	return 0;
}


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM