有没有快速的矩阵指数计算方法？

Question

除了简单的分治算法，有没有更快的矩阵指数计算方法来计算Mⁿ（其中M是一个矩阵，n是一个整数）？

- Akashdeep Saluja

---

嘿，我在stackoverflow上找到了一个链接，你可以看看：https://dev59.com/FGjWa4cB1Zd3GeqPmw-q - user1206604

Expokit是一个广为人知的用于执行矩阵指数运算的软件包。http://fortranwiki.org/fortran/show/Expokit - Sayan

Answer 1

你可以将矩阵分解为特征值和特征向量。这样，您就可以获得

M = V * D * V^-1

其中 V 是特征向量矩阵，D 是对角矩阵。将其提高至 N 次方，则可以得到如下结果：

M^n = (V * D * V^-1) * (V * D * V^-1) * ... * (V * D * V^-1)
    = V * D^n * V^-1

因为所有的V和V^-1项都会抵消。
由于D是对角线矩阵，你只需要将一堆（实数）数值提高到n次幂，而不是完整的矩阵。您可以在n的对数时间内完成此操作。 计算特征值和特征向量是r^3（其中r是M的行/列数）。根据r和n的相对大小，这可能更快或更慢。

- Michael

---

@AkashdeepSaluja：这比平方取幂要快。 这需要O（r ^ 3）时间，而平方取幂需要O（r ^ 3 logn）时间。 - Keith Randall

不必要，但足够。 - Valentin Waeselynck

@SinByCos 是的，但矩阵大小的对数不是吗？平方在指数上也是对数级别的，所以你不能直接比较两者。 - Stefan Dragnev

即使对于有缺陷的矩阵，您始终可以找到Jordan标准形式。然后，D不是对角线矩阵，而是对角线矩阵和一个幂零矩阵的和，您仍然可以非常高效地使用它。 - WorldSEnder

@WorldSEnder：不幸的是，Jordan标准形式在数值上并不稳定（标准形式是矩阵的不连续函数），因此在计算矩阵时产生的小舍入误差可能会导致结果出现大误差。 - Nate Eldredge

Answer 2

使用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--;
        }
    }
}

请查看以下数字的等效表：

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

- Gloomcore

Answer 3

平方取幂法 经常用于获取矩阵的高次幂。

- MBo

---

我知道这个方法，但需要进一步加快速度。 - Akashdeep Saluja

你最好在问题中添加该算法的名称，以避免类似的答案 :) - MBo

更快的算法通常更加复杂。 - Ari

Answer 4

我建议使用计算斐波那契数列的矩阵形式。据我所知，它的效率为O(log(n))。

- Rail Suleymanov

---

你需要将那个数字乘以矩阵乘法的成本。总运行时间为O(n^3 log n)。 - mrk