Calculating the exponential of a square matrix

Question

I'm trying to write a method that calculates the exponential of a square matrix. In this instance, the matrix is a square array of value:

[1 0]

[0 10]

and the method should return a value of:

[e 0]

[0 e^10]

However, when I run my code, I get a range of values depending on what bits I've rearranged, non particularly close to the expected value.

The way the method works is to utilise the power series for the matrix, so basically for a matrix A, n steps and an identity matrix I:

exp(A) = I + A + 1/2!*AA + 1/3!*AAA + ... +1/n!*AAA..

The code follows here. The method where I'm having the issue is the method exponential(Matrix A, int nSteps). The methods involved are enclosed, and the Matrix objects take the arguments (int m, int n) to create an array of size double[m][n].

    public static Matrix multiply(Matrix m1, Matrix m2){

    if(m1.getN()!=m2.getM()) return null;

    Matrix res = new Matrix(m1.getM(), m2.getN());

    for(int i = 0; i < m1.getM(); i++){
        for(int j = 0; j < m2.getN(); j++){
            res.getArray()[i][j] = 0;
            for(int k = 0; k < m1.getN(); k++){
                res.getArray()[i][j] = res.getArray()[i][j] + m1.getArray()[i][k]*m2.getArray()[k][j];
            }
        }
    }   
    return res;
}


public static Matrix identityMatrix(int M){

    Matrix id = new Matrix(M, M);

    for(int i = 0; i < id.getM(); i++){
        for(int j = 0; j < id.getN(); j++){
            if(i==j) id.getArray()[i][j] = 1;
            else id.getArray()[i][j] = 0;
        }
    }
    return id;
}


public static Matrix addMatrix(Matrix m1, Matrix m2){
    Matrix m3 = new Matrix(m1.getM(), m2.getN());
    
    for(int i = 0; i < m3.getM(); i++){
        for(int j = 0; j < m3.getN(); j++){
            m3.getArray()[i][j] = m1.getArray()[i][j] + m2.getArray()[i][j];
        }
    }       
    return m3;
}


public static Matrix scaleMatrix(Matrix m, double scale){
    Matrix res = new Matrix(m.getM(), m.getN());
    for(int i = 0; i < res.getM(); i++){
        for(int j = 0; j < res.getN(); j++){
            res.getArray()[i][j] = m.getArray()[i][j]*scale;
        }
    }
    return res;
}


public static Matrix exponential(Matrix A, int nSteps){
    
    Matrix runtot = identityMatrix(A.getM());
    Matrix sum = identityMatrix(A.getM());
    
    double factorial = 1.0;
    
    for(int i = 1; i <= nSteps; i++){

        sum = Matrix.multiply(Matrix.scaleMatrix(sum, factorial), A);
        runtot = Matrix.addMatrix(runtot, sum);
        factorial /= (double)i;

    }
    return runtot;
}

So my question is, how should I modify my code, so that I can input a matrix and a number of timesteps to calculate the exponential of said matrix after said timesteps?

Answer 1

My way to go would be to keep two accumulators :

• the sum, which is your approximation of exp(A)
• the nth term of the series M_n, that is A^n/n!

Note that there is a nice recursive relationship with M_n: M_{n+1} = M_n * A / (n+1)

Which yields :

public static Matrix exponential(Matrix A, int nSteps){

    Matrix seriesTerm = identityMatrix(A.getM());
    Matrix sum = identityMatrix(A.getM());

    for(int i = 1; i <= nSteps; i++){
        seriesTerm = Matrix.scaleMatrix(Matrix.multiply(seriesTerm,A),1.0/i);
        sum = Matrix.addMatrix(seriesTerm, sum);
    }
    return sum;
}

I totally understand the sort of thrill that implementing such algorithms can give you. But if this is not a hobby project, I concur that you should that you should use a library for this kind of stuff. Making such computations precise and efficient is really not a trivial matter, and a huge wheel to reinvent.