Matrix Inversion

Although it doesn't happen often, there may be times when you will need to determine the inverse of a matrix on its own rather than in the context of solving a system of equations. An example would be if you were using an implicit solution technique to perform a certain analysis. The gaussian() and luDecomp() methods we developed previously could be modified to compute a matrix inverse. Since the gaussJordan() method already calculates the matrix inverse, we will use it to create our matrix inverter method.

The modification is quite simple. All we do is remove the b vector as an input parameter and remove any operations performed on the b vector inside the method. The Gauss-Jordan elimination technique is used to convert the A matrix into its inverse. Remember that when the original A matrix is pivoted the computed inverse matrix will have its columns in the wrong order. We can unscramble the inverse matrix by interchanging its columns in the opposite order that the rows of the original A matrix were swapped. The resulting method is called invertMatrix() and its code listing is shown next .

public static void invertMatrix(double a[][]) {
  int i,j,k,m;
  double temp;

  int numRows = a.length;
  int numCols = a[0].length;
  int index[][] = new int[numRows][2];

  //  Perform an implicit partial pivoting of the
  //  a[][] array. We will provide a dummy b  array
  //  to the partialPivot() method.

  partialPivot(a, new double[numRows], index);

  //  Perform the elimination row by row. First dividing
  //  the current row by a[i][i]

  for(i=0; i<numRows; ++i) {
    temp = a[i][i];
    for(j=0; j<numCols; ++j) {
      a[i][j] /= temp;
    }
    a[i][i] = 1.0/temp;

    //  Reduce the other rows by subtracting a multiple
    //  of the current row from them. Don't reduce the
    //  current row. As each column of the a[][] matrix
    //  is reduced its elements are replaced with the
    //  inverse a[][] matrix.

    for(k=0; k<numRows; ++k) {
      if (k != i) {
        temp = a[k][i];
        for(j=0; j<numCols; ++j) {
          a[k][j] -= temp*a[i][j];
        }
        a[k][i] = -temp*a[i][i];
      }
    }
  }

  //  Unscramble the inverse a[][] matrix.
  //  The columns are swapped in the opposite order
  //  that the rows were during the pivoting.

  for(j=numCols-1; j>=0; j) {
    k = index[j][0];
    m = index[j][1];
    if (k != m) {
      for(i=0; i<numRows; ++i) {
        temp = a[i][m];
        a[i][m] = a[i][k];
        a[i][k] = temp;
      }
    }
  }

  return;
}

Testing the `EqnSolver` Class Methods

The methods we wrote look very nice, but we need to test them to make sure everything is working properly. First we have to compile the EqnSolver.java source code. Because the EqnSolver class was placed in the TechJava .MathLib package, you will have to copy the EqnSolver.java source code to the TechJava\MathLib directory and compile the code from the package root directory as we have done previously. If you followed the examples in Chapter 17, your CLASSPATH should already include the package root directory.

We will then write a driver program named SolverDemo.java . This program simply loads the A matrix and b vector from our test case and calls the gaussJordan() , gaussian() , or luDecomp() method. The resulting solution vector is written out. The SolverDemo source code is shown here.

import TechJava.MathLib.EqnSolver;

public class SolverDemo
{
  public static void main(String args[]) {
    double a[][] = { {0.0, 1.0, 5.0, -2.0},
                     {7.0, -6.0, 3.0, 1.0},
                     {1.0, 4.0, 2.0, -3.0},
                     {5.0, -2.0, 1.0, 4.0} };
    double b[] = { 9.0, 8.0, 3.0, 20.0 };

    EqnSolver.gaussJordan(a,b);
//    EqnSolver.gaussian(a,b);
//    EqnSolver.luDecomp(a,b);

    for (int i=0; i<a.length; ++i) {
      System.out.println("b["+i+"] = "+b[i]);
    }
  }
}

Output ”

b[0] = 1.000000000004
b[1] = 2.0
b[2] = 3.0
b[3] = 4.0

In this case the gaussJordan() method was called. There was a slight round-off error in the first unknown value, but the other three values came out cleanly. You will find that there is usually a small round-off error when using floating-point math. If you try calling the gaussian() and luDecomp() methods, you will see that they also return the correct solution vector.

To test the invertMatrix() method, we will create another driver program named InvertDemo . The main() method from this class computes the inverse of the A matrix from our test case using the invertMatrix() method. It then computes the product of the original A matrix with its computed inverse to see if it equals the identity matrix.

import TechJava.MathLib.EqnSolver;

public class InvertDemo
{
  public static void main(String args[]) {
    double a[][] = { {0.0, 1.0, 5.0, -2.0},
                     {7.0, -6.0, 3.0, 1.0},
                     {1.0, 4.0, 2.0, -3.0},
                     {5.0, -2.0, 1.0, 4.0} };
    int i,j,k;

    //  Compute the inverse of the A matrix.

    EqnSolver.invertMatrix(a);

    //  Re-create the original A matrix, compute
    //  the product of it and its inverse, and
    //  see if it equals the identity matrix.

    double A[][] = { {0.0, 1.0, 5.0, -2.0},
                     {7.0, -6.0, 3.0, 1.0},
                     {1.0, 4.0, 2.0, -3.0},
                     {5.0, -2.0, 1.0, 4.0} };
    double I[][] = { {0.0, 0.0, 0.0, 0.0},
                     {0.0, 0.0, 0.0, 0.0},
                     {0.0, 0.0, 0.0, 0.0},
                     {0.0, 0.0, 0.0, 0.0} };

    int numRows = A.length;
    int numCols = A[0].length;

    for(i=0; i<numRows; ++i) {
      for(j=0; j<numCols; ++j) {
        for(k=0; k<numRows; ++k) {
          I[i][j] += A[i][k]*a[k][j];
        }
      }
    }

    for(i=0; i<numRows; ++i) {
      for(j=0; j<numCols; ++j) {
        System.out.println("I["+i+"]["+j+"] = "+
                               I[i][j]);
      }
    }
  }
}

Output ”

I[0][0] = 0.9999999999
I[0][1] = 0.0
I[0][2] = 0.0
I[0][3] = 0.0
I[1][0] = 1.5356534878e-16
I[1][1] = 0.9999999999
I[1][2] = 0.0
I[1][3] = 1.1022334581e-16
I[2][0] = 0.0
I[2][1] = 0.0
I[2][2] = 1.0
I[2][3] = 1.1022334581e-16
I[3][0] = 2.2204460492e-16
I[3][1] = 0.0
I[3][2] = 1.1022334581e-16
I[3][3] = 1.0

While the SolverDemo and InvertDemo classes only test the EqnSolver class methods with one 4x4 system of equations, the three methods were extensively tested with a variety of A matrices and b vectors during the development of this chapter. You are encouraged to modify the SolverDemo class to perform your own testing of the EqnSolver class methods.