1. 程式人生 > >矩陣的特征值分解

矩陣的特征值分解

向上 過大 部分 tro 開始 namespace sin += 通過

引入問題:給定一個對角線非零的上三角矩陣\(M\),求\(M^k\),滿足\(M\)的階\(\le 500\)\(k\le 10^9\)

對998244353取模。

一個顯而易見的算法是矩陣快速冪,然而是\(O(N^3\log k)\)的,無法通過本題。

一開始我想,既然是上三角矩陣,那麽特征多項式一定不難求,那麽是用CH定理+FFT多項式取模啥搞搞?

然而我naive了。

這題我們可以把\(M\)特征值分解為\(Q^{-1}AQ\)形式,其中\(A\)是一個對角矩陣。

那麽\(M^k=(Q^{-1}AQ)^k=Q^{-1}A^kQ\)

對一個對角矩陣進行冪的復雜度是\(N\log C\)的,矩陣乘法的復雜度是\(O(N^3)\)

的,對一個上三角矩陣進行特征值分解可以使用高斯消元,時間復雜度也是\(O(N^3)\),具體怎麽對上三角矩陣進行特征值分解??我tm怎麽知道,這個得好好研究一下

upd

自己手推沒推出來。觀察了下std,手跑了下樣例,得出來一些性質。

矩陣\(Q^{-1}\)的第\(i\)列,即為矩陣\(M\)對應第\(i\)行第\(i\)列特征值的特征向量。

這個性質通過特征值分解那套理論也不難得到--因為特征向量是\(M\)所對應“方向不變”的向量,而\(Q\)\(Q^{-1}\)就是在這些旋轉方向上的向量,通過線性變換把它們旋轉過去//線代那套理論太玄學

std裏在\(Q^{-1}\)上遞推的,沒有看得非常透徹(不過大致也觀察出了一些什麽)

,目前已經觀察得比較透徹了。

求出\(Q^{-1}\)後直接上矩陣求逆板子求\(Q\),然後直接矩陣乘法就行了。

代碼如下:

#include <cstdio>
using namespace std;

const int p = 998244353;

int n, k, a[500][500], l[500][500], r[500][500], v[500];

int qpow(int x, int y)
{
    int res = 1;
    for (x %= p; y > 0; y >>= 1, x = x * (long long)x % p)
        if (y & 1) res = res * (long long)x % p;
    return res;
}

int main()
{
    scanf("%d%d", &n, &k);
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            scanf("%d", &a[i][j]);
    
    //求l
    for (int i = 0; i < n; i++) //枚舉l矩陣的第i列,為a矩陣對應於aii的特征向量
    {
        l[i][i] = 1; //欽定的
        for (int j = i - 1; j >= 0; j--) //求這個特征向量的第j行的值
        {
            int sum = 0; //這裏需要滿足的是\sum_{k=0}^{n-1}a_{j,k}*l{k,i}=0,此時求值為$b_{ji}$
            for (int k = j + 1; k <= i; k++)
                sum = (sum + a[j][k] * (long long)l[k][i]) % p;
            l[j][i] = sum * (long long)qpow((a[i][i] - a[j][j] + p) % p, p - 2) % p;
            //註意這裏是a[i][i] - a[j][j], 相當於乘了個-1,就是我們要求的值了
        }
    }
    
    //求l的逆矩陣r,註意到l是上三角矩陣
    for (int i = n - 1; i >= 0; i--)
    {
        r[i][i] = qpow(l[i][i], p - 2);
        for (int j = 0; j < i; j++)
        {
            int e = l[j][i] * (long long)qpow(l[j][j], p - 2) % p;
            for (int k = i; k < n; k++)
                r[j][k] = ((r[j][k] - r[i][k] * (long long)e % p) % p + p) % p;
        }
    }
    
    //收集答案
    for (int i = 0; i < n; i++) v[i] = qpow(a[i][i], k);
    
    long long ans1 = 0, ans2 = 0;
    
    for (int i = 0; i < n; i++)
        for (int j = i; j < n; j++)
        {
            int sb = 0;
            for (int k = i; k <= j; k++)
                sb = (sb + l[i][k] * (long long)v[k] % p * r[k][j] % p) % p;
            ans1 += sb, ans2 ^= sb;
        }
    printf("%lld %lld\n", ans1, ans2);
    return 0;
}

以後再研究下一般矩陣的特征值分解,就可以弄圖片壓縮啥的了。

這個代碼常數略大,本來可以弄小一點的

把最後收集答案時候先讓對角矩陣和l乘一下再收集、以及優化一下(x%p+p)%p那部分即可。

矩陣的特征值分解