I'm trying to calculate the inverse ma开发者_StackOverflow社区trix in Java.
I'm following the adjoint method (first calculation of the adjoint matrix, then transpose this matrix and finally, multiply it for the inverse of the value of the determinant).
It works when the matrix is not too big. I've checked that for matrixes up to a size of 12x12 the result is quickly provided. However, when the matrix is bigger than 12x12 the time it needs to complete the calculation increases exponentially.
The matrix I need to invert is 19x19, and it takes too much time. The method that more time consumes is the method used for the calculation of the determinant.
The code I'm using is:
public static double determinant(double[][] input) {
int rows = nRows(input); //number of rows in the matrix
int columns = nColumns(input); //number of columns in the matrix
double determinant = 0;
if ((rows== 1) && (columns == 1)) return input[0][0];
int sign = 1;
for (int column = 0; column < columns; column++) {
double[][] submatrix = getSubmatrix(input, rows, columns,column);
determinant = determinant + sign*input[0][column]*determinant(submatrix);
sign*=-1;
}
return determinant;
}
Does anybody know how to calculate the determinant of a large matrix more efficiently? If not, does anyone knows how to calcultate the inverse of a large matrix using other algorithm?
Thanks
Exponentially? No, I believe matrix inversion is O(N^3).
I would recommend using LU decomposition to solve a matrix equation. You don't have to solve for the determinant when you use it.
Better yet, look into a package to help you. JAMA comes to mind.
12x12 or 19x19 are not large matricies. It's common to solve problems with tens or hundreds of thousands of degrees of freedom.
Here's a working example of how to use JAMA. You have to have the JAMA JAR in your CLASSPATH when you compile and run:
package linearalgebra;
import Jama.LUDecomposition;
import Jama.Matrix;
public class JamaDemo
{
public static void main(String[] args)
{
double [][] values = {{1, 1, 2}, {2, 4, -3}, {3, 6, -5}}; // each array is a row in the matrix
double [] rhs = { 9, 1, 0 }; // rhs vector
double [] answer = { 1, 2, 3 }; // this is the answer that you should get.
Matrix a = new Matrix(values);
a.print(10, 2);
LUDecomposition luDecomposition = new LUDecomposition(a);
luDecomposition.getL().print(10, 2); // lower matrix
luDecomposition.getU().print(10, 2); // upper matrix
Matrix b = new Matrix(rhs, rhs.length);
Matrix x = luDecomposition.solve(b); // solve Ax = b for the unknown vector x
x.print(10, 2); // print the solution
Matrix residual = a.times(x).minus(b); // calculate the residual error
double rnorm = residual.normInf(); // get the max error (yes, it's very small)
System.out.println("residual: " + rnorm);
}
}
Here's the same problem solved using Apache Commons Math, per quant_dev's recommendation:
package linearalgebra;
import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.linear.ArrayRealVector;
import org.apache.commons.math.linear.DecompositionSolver;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealVector;
public class LinearAlgebraDemo
{
public static void main(String[] args)
{
double [][] values = {{1, 1, 2}, {2, 4, -3}, {3, 6, -5}};
double [] rhs = { 9, 1, 0 };
RealMatrix a = new Array2DRowRealMatrix(values);
System.out.println("a matrix: " + a);
DecompositionSolver solver = new LUDecompositionImpl(a).getSolver();
RealVector b = new ArrayRealVector(rhs);
RealVector x = solver.solve(b);
System.out.println("solution x: " + x);;
RealVector residual = a.operate(x).subtract(b);
double rnorm = residual.getLInfNorm();
System.out.println("residual: " + rnorm);
}
}
Adapt these to your situation.
I would recommend using Apache Commons Math 2.0 for this. JAMA is a dead project. ACM 2.0 actually took linear algebra from JAMA and developed it further.
The la4j (Linear Algebra for Java) library supports matrix inversion. Here is the brief example:
Matrix a = new Basic2DMatrix(new double[][]{
{ 1.0, 2.0, 3.0 },
{ 4.0, 5.0, 6.0 },
{ 7.0, 8.0. 9.0 }
});
Matrix b = a.invert(Matrices.DEFAULT_INVERTOR); // uses Gaussian Elimination
For those who seeks matrix inversion (not fast), see https://github.com/rchen8/Algorithms/blob/master/Matrix.java.
import java.util.Arrays;
public class Matrix {
private static double determinant(double[][] matrix) {
if (matrix.length != matrix[0].length)
throw new IllegalStateException("invalid dimensions");
if (matrix.length == 2)
return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
double det = 0;
for (int i = 0; i < matrix[0].length; i++)
det += Math.pow(-1, i) * matrix[0][i]
* determinant(submatrix(matrix, 0, i));
return det;
}
private static double[][] inverse(double[][] matrix) {
double[][] inverse = new double[matrix.length][matrix.length];
// minors and cofactors
for (int i = 0; i < matrix.length; i++)
for (int j = 0; j < matrix[i].length; j++)
inverse[i][j] = Math.pow(-1, i + j)
* determinant(submatrix(matrix, i, j));
// adjugate and determinant
double det = 1.0 / determinant(matrix);
for (int i = 0; i < inverse.length; i++) {
for (int j = 0; j <= i; j++) {
double temp = inverse[i][j];
inverse[i][j] = inverse[j][i] * det;
inverse[j][i] = temp * det;
}
}
return inverse;
}
private static double[][] submatrix(double[][] matrix, int row, int column) {
double[][] submatrix = new double[matrix.length - 1][matrix.length - 1];
for (int i = 0; i < matrix.length; i++)
for (int j = 0; i != row && j < matrix[i].length; j++)
if (j != column)
submatrix[i < row ? i : i - 1][j < column ? j : j - 1] = matrix[i][j];
return submatrix;
}
private static double[][] multiply(double[][] a, double[][] b) {
if (a[0].length != b.length)
throw new IllegalStateException("invalid dimensions");
double[][] matrix = new double[a.length][b[0].length];
for (int i = 0; i < a.length; i++) {
for (int j = 0; j < b[0].length; j++) {
double sum = 0;
for (int k = 0; k < a[i].length; k++)
sum += a[i][k] * b[k][j];
matrix[i][j] = sum;
}
}
return matrix;
}
private static double[][] rref(double[][] matrix) {
double[][] rref = new double[matrix.length][];
for (int i = 0; i < matrix.length; i++)
rref[i] = Arrays.copyOf(matrix[i], matrix[i].length);
int r = 0;
for (int c = 0; c < rref[0].length && r < rref.length; c++) {
int j = r;
for (int i = r + 1; i < rref.length; i++)
if (Math.abs(rref[i][c]) > Math.abs(rref[j][c]))
j = i;
if (Math.abs(rref[j][c]) < 0.00001)
continue;
double[] temp = rref[j];
rref[j] = rref[r];
rref[r] = temp;
double s = 1.0 / rref[r][c];
for (j = 0; j < rref[0].length; j++)
rref[r][j] *= s;
for (int i = 0; i < rref.length; i++) {
if (i != r) {
double t = rref[i][c];
for (j = 0; j < rref[0].length; j++)
rref[i][j] -= t * rref[r][j];
}
}
r++;
}
return rref;
}
private static double[][] transpose(double[][] matrix) {
double[][] transpose = new double[matrix[0].length][matrix.length];
for (int i = 0; i < matrix.length; i++)
for (int j = 0; j < matrix[i].length; j++)
transpose[j][i] = matrix[i][j];
return transpose;
}
public static void main(String[] args) {
// example 1 - solving a system of equations
double[][] a = { { 1, 1, 1 }, { 0, 2, 5 }, { 2, 5, -1 } };
double[][] b = { { 6 }, { -4 }, { 27 } };
double[][] matrix = multiply(inverse(a), b);
for (double[] i : matrix)
System.out.println(Arrays.toString(i));
System.out.println();
// example 2 - example 1 using reduced row echelon form
a = new double[][]{ { 1, 1, 1, 6 }, { 0, 2, 5, -4 }, { 2, 5, -1, 27 } };
matrix = rref(a);
for (double[] i : matrix)
System.out.println(Arrays.toString(i));
System.out.println();
// example 3 - solving a normal equation for linear regression
double[][] x = { { 2104, 5, 1, 45 }, { 1416, 3, 2, 40 },
{ 1534, 3, 2, 30 }, { 852, 2, 1, 36 } };
double[][] y = { { 460 }, { 232 }, { 315 }, { 178 } };
matrix = multiply(
multiply(inverse(multiply(transpose(x), x)), transpose(x)), y);
for (double[] i : matrix)
System.out.println(Arrays.toString(i));
}
}
Matrix inversion is computationally very intensive. As duffymo answered LU is a good algorithm, and there are other variants (QR, for instance).
Unfortunately you can't get rid of the heavy calculations... and maybe the bottelneck is the getSubmatrix method if you are not using an optimized library.
Also, special matrix structures (band-matricity, symmetry, diagonality, sparsity) have a great impact in performance if considered in the calculations. Your mileage may vary...
You NEVER want to compute an inverse matrix this way. Ok, computation of the inverse itself is to be avoided, as it is almost always better to use a factorization such as an LU.
Computation of the determinant using recursive computations is a numerically obscene thing to do. It turns out that a better choice is to use an LU factorization to compute a determinant. But, if you are going to bother to compute LU factors, then why would you possibly want to compute the inverse? You have already done the difficult work by computing the LU factors.
Once you have LU factors, you can use them to do back and forward substitution.
As far as a 19x19 matrix being big, it is not even close to what I'd think of as big.
Since ACM library has updated over the years, here is the implementation conforming to latest definition for matrix inversion.
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.DecompositionSolver;
import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
public class LinearAlgebraDemo
{
public static void main(String[] args)
{
double [][] values = {{1, 1, 2}, {2, 4, -3}, {3, 6, -5}};
double [][] rhs = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
// Solving AB = I for given A
RealMatrix A = new Array2DRowRealMatrix(values);
System.out.println("Input A: " + A);
DecompositionSolver solver = new LUDecomposition(A).getSolver();
RealMatrix I = new Array2DRowRealMatrix(rhs);
RealMatrix B = solver.solve(I);
System.out.println("Inverse B: " + B);
}
}
Your algorithm to compute a determinant is indeed exponential. The basic problem is that you are computing from the definition, and the straight definition leads to an exponential amount of subdeterminants to compute. You really need to transform the matrix first before computing either its determinant or its inverse. (I thought of explaining about dynamic programming, but this problem cannot be solved by dynamic programming as the number of subproblems is exponential too.)
LU decomposition, as recommended by others, is a good choice. If you are new to matrix calculation, you might also want to look at Gaussian elimination to compute determinants and inverses, as that might be a bit easier to comprehend at first.
And one thing to remember in matrix inversion is numerical stability, since you are dealing with floating point numbers. All the good algorithm include permutations of rows and/or columns to choose the suitable pivots, as they are called. At least in Gaussian elimination, you want to, at each step, to permute the columns so that the element largest in absolute value is chosen as the pivot, as this is the stablest choice.
It is tough to beat Matlab at their game. They also are cautiou about precision. If you have 2.0 and 2.00001 as pivots - watch out! Your answer might end up being very imprecise. Also, check out Python's implementation (it is in numpy / scipy somewhere ...)
Do you have to have an exact solution? An approximate solver (Gauss-Seidel is very performant and easy to implement) will likely work for you, and will converge very quickly. 19x19 is quite a small matrix. I think the Gauss-Seidel code that I used could solve a 128x128 matrix in the blink of an eye (but don't quote me on that, it's been a while).
By ACM you can do this without rhs:
RealMatrix A = new Array2DRowRealMatrix(ARR);
System.out.println(new LUDecomposition(A).getSolver().getInverse());
![Interactive visualization of a graph in python [closed]](https://www.devze.com/res/2023/04-10/09/92d32fe8c0d22fb96bd6f6e8b7d1f457.gif)
加载中,请稍侯......
精彩评论