Inversing a matrix

[debee]

Hi there,

I got a little problem with finding the inverse of a matrix when using the BigDecimal type. I mean, I have got a method that finds the inverse of a matrix perfectly with the primitive type double. Here is the method as shown below:

import java.lang.*;
public class Inverse {
  public static void main(String arg[]) {

    double a[][]= {{0.0038053780181931084, 0.003729663276441426, 0.0064344762767650256, 0.0010346528952388523, -0.026667616352678624, -0.029400456101594908},

{0.003729663276441426, 0.0038541077016533546, 0.00647664820019812, 0.001041434069152789, -0.026842397426785514, -0.029593148362827063},

{0.0064344762767650256, 0.00647664820019812, 0.011273619736763435, 0.0017966991428170357, -0.046308944441478844, -0.05105458508765971},

{0.0010346528952388523, 0.001041434069152789, 0.0017966991428170357, 3.889061813315775E-4, -0.007446399890360643, -0.008209490874480767},

{-0.026667616352678624, -0.026842397426785514, -0.046308944441478844, -0.007446399890360643, 0.19202691229934035, 0.21159516790500957},

{-0.029400456101594908, -0.029593148362827063, -0.05105458508765971, -0.008209490874480767, 0.21159516790500957, 0.23337898388173603},

};
    int n = a.length;
    double d[][] = invert(a);
    for (int i=0; i<a.length; ++i)
      for (int j=0; j<a[i].length; ++j)
        {
			System.out.print(d[i][j]+ "  ");

			  if(j == a[i].length - 1)
			  {
			   System.out.println();
			   System.out.println();
		     }
		}
  }

  public static double[][] invert(double a[][]) {
    int n = a.length;
    double x[][] = new double[n][n];
    double b[][] = new double[n][n];
    int index[] = new int[n];
    for (int i=0; i<n; ++i) b[i][i] = 1;

 // Transform the matrix into an upper triangle
    gaussian(a, index);

 // Update the matrix b[i][j] with the ratios stored
    for (int i=0; i<n-1; ++i)
      for (int j=i+1; j<n; ++j)
        for (int k=0; k<n; ++k)
          b[index[j]][k]
            -= a[index[j]][i]*b[index[i]][k];

 // Perform backward substitutions
    for (int i=0; i<n; ++i) {
      x[n-1][i] = b[index[n-1]][i]/a[index[n-1]][n-1];
      for (int j=n-2; j>=0; --j) {
        x[j][i] = b[index[j]][i];
        for (int k=j+1; k<n; ++k) {
          x[j][i] -= a[index[j]][k]*x[k][i];
        }
        x[j][i] /= a[index[j]][j];
      }
    }
  return x;
  }

// Method to carry out the partial-pivoting Gaussian
// elimination.  Here index[] stores pivoting order.

  public static void gaussian(double a[][],
    int index[]) {
    int n = index.length;
    double c[] = new double[n];

 // Initialize the index
    for (int i=0; i<n; ++i) index[i] = i;

 // Find the rescaling factors, one from each row
    for (int i=0; i<n; ++i) {
      double c1 = 0;
      for (int j=0; j<n; ++j) {
        double c0 = Math.abs(a[i][j]);
        if (c0 > c1) c1 = c0;
      }
      c[i] = c1;
    }

 // Search the pivoting element from each column
    int k = 0;
    for (int j=0; j<n-1; ++j) {
      double pi1 = 0;
      for (int i=j; i<n; ++i) {
        double pi0 = Math.abs(a[index[i]][j]);
        pi0 /= c[index[i]];
        if (pi0 > pi1) {
          pi1 = pi0;
          k = i;
        }
      }

   // Interchange rows according to the pivoting order
      int itmp = index[j];
      index[j] = index[k];
      index[k] = itmp;
      for (int i=j+1; i<n; ++i) {
        double pj = a[index[i]][j]/a[index[j]][j];

     // Record pivoting ratios below the diagonal
        a[index[i]][j] = pj;

     // Modify other elements accordingly
        for (int l=j+1; l<n; ++l)
          a[index[i]][l] -= pj*a[index[j]][l];
      }
    }
  }
}

But when I changed the types to be that of BigDecimal and i called the invert method, as shown below:

public static BigDecimal[][] invert(BigDecimal a[][]) {
     int n = a.length;
     BigDecimal x[][] =new BigDecimal[n][n];
     BigDecimal b[][] = new BigDecimal[n][n];
     int index[] = new int[n];
     for (int i=0; i<n; ++i) b[i][i] = new BigDecimal("1");

  // Transform the matrix into an upper triangle
     gaussian(a, index);

  // Update the matrix b[i][j] with the ratios stored
     for (int i=0; i<n-1; ++i)
       for (int j=i+1; j<n; ++j)
         for (int k=0; k<n; ++k)
           b[index[j]][k]
            =b[index[j]][k].subtract( a[index[j]][i].multiply (b[index[i]][k]));

  // Perform backward substitutions
     for (int i=0; i<n; ++i) {
       x[n-1][i] = b[index[n-1]][i].divide(a[index[n-1]][n-1],20,RoundingMode.HALF_UP);
       for (int j=n-2; j>=0; --j) {
         x[j][i] = b[index[j]][i];
         for (int k=j+1; k<n; ++k) {
           x[j][i] =x[j][i].subtract( a[index[j]][k].multiply(x[k][i]));
         }
         x[j][i] = x[j][i].divide(a[index[j]][j],20,RoundingMode.HALF_UP);
       }
     }
   return x;
   }

 // Method to carry out the partial-pivoting Gaussian
 // elimination.  Here index[] stores pivoting order.

   public static void gaussian(BigDecimal a[][],
     int index[]) {
     int n = index.length;
     BigDecimal c[] = new BigDecimal[n];

  // Initialize the index
     for (int i=0; i<n; ++i) index[i] = i;

  // Find the rescaling factors, one from each row
     for (int i=0; i<n; ++i) {
       BigDecimal c1 = new BigDecimal("0");
       for (int j=0; j<n; ++j) {
         BigDecimal c0 = (a[i][j]).abs();
         if (c0.compareTo(c1)>0) c1 = c0;
       }
       c[i] = c1;
     }

  // Search the pivoting element from each column
     int k = 0;
     for (int j=0; j<n-1; ++j) {
       BigDecimal pi1 = new BigDecimal("0");
       for (int i=j; i<n; ++i) {
         BigDecimal pi0 =a[index[i]][j].abs();
         pi0 =pi0.divide(c[index[i]],20,RoundingMode.HALF_UP);
         if (pi0.compareTo (pi1) >0) {
           pi1 = pi0;
           k = i;
         }
       }

    // Interchange rows according to the pivoting order
       int itmp = index[j];
       index[j] = index[k];
       index[k] = itmp;
       for (int i=j+1; i<n; ++i) {
         BigDecimal pj = a[index[i]][j].divide(a[index[j]][j],20,RoundingMode.HALF_UP);

      // Record pivoting ratios below the diagonal
         a[index[i]][j] = pj;

      // Modify other elements accordingly
         for (int l=j+1; l<n; ++l)
           a[index[i]][l] =a[index[i]][l].subtract( pj.multiply(a[index[j]][l]));
       }
     }
  }
  public void inverse()

  {
  	int i, j, k;

  	System.out.println();
  	System.out.println("The weightMatrix after inversion is:");
  	System.out.println();


  	   weightMatrix1 = invert(weightMatrix);
  			     for (i=0; i<weightMatrix1.length; ++i)
  			       for (j=0; j<weightMatrix1[i].length; ++j)
  			         {
  			 			System.out.print(weightMatrix1[i][j]+ "  ");

  			 			  if(j == weightMatrix1[i].length - 1)
  			 			  {
  			 			   System.out.println();
  			 			   System.out.println();
  
			 		      }

}
}

it gives me an error saying:

Exception in thread "main" java.lang.NullPointerException
        at newtonsMethod1.invert(newtonsMethod1.java:282)
        at newtonsMethod1.inverse(newtonsMethod1.java:359)

the line 282 at newtonsMethod1.invert is this bit:

b[index[j]][k]

in the following part of the code:

Update the matrix b[i][j] with the ratios stored
     for (int i=0; i<n-1; ++i)
       for (int j=i+1; j<n; ++j)
         for (int k=0; k<n; ++k)
           b[index[j]][k]
            =b[index[j]][k].subtract( a[index[j]][i].multiply (b[index[i]][k]));

Can you figure out what is going on here, cos i have only this problem when I changed the types from double to BigDecimal, and it worked well with type double.
Thanks

[Ezzaral]

That line is actually the entirety of this statement

b[index[j]][k] = b[index[j]][k].subtract( a[index[j]][i].multiply (b[index[i]][k]));

You are performing an operation on a value in b[][], but you haven't initialized all of the elements in b[][]. You have only initialized the diagonal of that matrix here

for (int i=0; i<n; ++i) b[i][i] = new BigDecimal("1");

The double array used previously had default values of 0, but with the BigDecimal version you have nulls in all but the diagonal. You need to first initialize the values in b[][] to valid BigDecimal instances for 0.

If your matrix sizes are constant in size, you may be able to make that initialization easier and faster with System.arraycopy() from a single BigDecimal(0) template array, but then that is only a consideration after you get the thing working properly first.

[debee]

Thanks for your help,

I made the initialization with the following code:

int index[] = new int[n];
     for (int i=0; i<a.length; ++i)
        for(int j =0; j<a[i].length; j++)
         {
			 if(i==j)
			   b[i][i] = new BigDecimal("1");
			 else
			   b[i][j] = new BigDecimal("0");
	     }

and I did not have that problem again. I did not really understand what you meant by using System.ArrayCopy() to do the initialization, but I believe what i used is just as good, is it?

[Ezzaral]

Yes, unless you find a reason to speed up that array initialization, that will certainly work just fine. My point was that the identity matrix is constant for a given matrix size and could be cached for reuse instead of creating it each time the inversion routine is called.

If you are coding all of your matrix operations from scratch, you may want to consider using a package like Apache Commons Math.

[debee]

Hi there,

I have got another problem which i can't figure out why. In my algorithm, i have this method that takes in a BigDecimal value and returns a BigDecimal value as shown below:

public BigDecimal sigmoidActivationFunction(BigDecimal sum)
       {
   	  Big_pi test2 = new Big_pi();
         //sum = sum.multiply(new BigDecimal(-1.0));
     	  BigDecimal Set =new BigDecimal("1").add( test2.exp(sum.negate()));
     	  BigDecimal Set1 = new BigDecimal("1").divide(Set,20,RoundingMode.HALF_UP);

   	  return Set1;
     }

In the first iteration of the algorithm, which takes in the value -0.486449196566058673, it returns the right value as 0.41394292596487729849 but in the second iteration of the algorithm where the value
-450.342245610786369217517941611560860293742635964248089208715900796008344070875 is produced and passed into the method, it returned the value 0E-20 which is wrong. Could you tell me why, please? or is it because this valued was not rounded to a scale?