快速矩阵指数

Aka*_*uja 19 algorithm

有没有比简单的分治算法更快的矩阵求幂方法来计算M ^ n(其中M是矩阵,n是整数).

小智 25

您可以将矩阵分解为特征值和特征向量.然后你得到

M = V^-1 * D * V
Run Code Online (Sandbox Code Playgroud)

其中V是特征向量矩阵,D是对角矩阵.为了将此提高到N次幂,您会得到类似的结果:

M^n = (V^-1 * D * V) * (V^-1 * D * V) * ... * (V^-1 * D * V)
    = V^-1 * D^n * V
Run Code Online (Sandbox Code Playgroud)

因为所有V和V ^ -1项都取消了.

由于D是对角线,你只需要将一堆(实数)数字提升到第n个幂,而不是完整的矩阵.你可以在n的对数时间内做到这一点.

计算特征值和特征向量是r ^ 3(其中r是M的行/列的数量).根据r和n的相对大小,这可能更快或更快.

  • @WorldSEnder:不幸的是,Jordan 范式在数值上不稳定(范式是矩阵的不连续函数),因此计算矩阵时的小舍入误差可能会导致结果出现大误差。 (2认同)

Glo*_*ore 8

使用Euler快速功率算法非常简单.使用下一个算法.

#define SIZE 10

//It's simple E matrix
// 1 0 ... 0
// 0 1 ... 0
// ....
// 0 0 ... 1
void one(long a[SIZE][SIZE])
{
    for (int i = 0; i < SIZE; i++)
        for (int j = 0; j < SIZE; j++)
            a[i][j] = (i == j);
}

//Multiply matrix a to matrix b and print result into a
void mul(long a[SIZE][SIZE], long b[SIZE][SIZE])
{
    long res[SIZE][SIZE] = {{0}};

    for (int i = 0; i < SIZE; i++)
        for (int j = 0; j < SIZE; j++)
            for (int k = 0; k < SIZE; k++)
            {
                res[i][j] += a[i][k] * b[k][j];
            }

    for (int i = 0; i < SIZE; i++)
        for (int j = 0; j < SIZE; j++)
            a[i][j] = res[i][j];
}

//Caluclate a^n and print result into matrix res
void pow(long a[SIZE][SIZE], long n, long res[SIZE][SIZE])
{
    one(res);

    while (n > 0) {
        if (n % 2 == 0)
        {
            mul(a, a);
            n /= 2;
        }
        else {
            mul(res, a);
            n--;
        }
    }
}
Run Code Online (Sandbox Code Playgroud)

请在下面找到相应的数字:

long power(long num, long pow)
{
    if (pow == 0) return 1;
    if (pow % 2 == 0)
        return power(num*num, pow / 2);
    else
        return power(num, pow - 1) * num;
}
Run Code Online (Sandbox Code Playgroud)