Is there any fast method of matrix exponentiation?

Question

Is there any faster method of matrix exponentiation to calculate Mⁿ (where M is a matrix and n is an integer) than the simple divide and conquer algorithm?

Solution 1

You could factor the matrix into eigenvalues and eigenvectors. Then you get

M = V^-1 * D * V

Where V is the eigenvector matrix and D is a diagonal matrix. To raise this to the Nth power, you get something like:

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

Because all the V and V^-1 terms cancel.

Since D is diagonal, you just have to raise a bunch of (real) numbers to the nth power, rather than full matrices. You can do that in logarithmic time in n.

Calculating eigenvalues and eigenvectors is r^3 (where r is the number of rows/columns of M). Depending on the relative sizes of r and n, this might be faster or not.

Solution 2

It's quite simple to use Euler fast power algorith. Use next algorith.

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

Below please find equivalent for numbers:

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

Solution 3

Exponentiation by squaring is frequently used to get high powers of matrices.