Java矩阵求逆计算

Question

我正在尝试在Java中计算逆矩阵。

我正在遵循伴随矩阵法（首先计算伴随矩阵，然后转置此矩阵，最后将其乘以行列式值的倒数）。

当矩阵不太大时，它可以正常工作。我已经检查过，对于大小不超过12x12的矩阵，结果可以迅速提供。然而，当矩阵大于12x12时，完成计算所需的时间呈指数增长。

我需要求逆的矩阵大小为19x19，这需要太长时间。消耗更多时间的方法是用于计算行列式的方法。

我正在使用以下代码：

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

有人知道如何更高效地计算大矩阵的行列式吗？如果不知道，有人知道如何使用其他算法计算大矩阵的逆矩阵吗？

谢谢。

- dedalo

---

@duffymo：感谢您的回答。您说得对，不是指数级增长，我的意思是从矩阵大小为12x12开始，计算行列式所需的时间会显著增加。我尝试过Jama，但我无法让它正常工作（我在Java方面还很新）。我也会研究LU分解。谢谢。 - dedalo

不，"指数级"是正确的，因为您的算法确实是指数级的，但是duffymo也是正确的，矩阵求逆或行列式计算并不需要指数级时间。 - JaakkoK

感谢大家。我正在研究Jama，并在Matrix类中找到了可以快速计算行列式的方法'det'。我还发现了计算矩阵L和U（A = L*U）的方法，然后det(A) = det(L)*det(U)。 - dedalo

Answer 1

指数级的算法？不是的，我认为矩阵求逆的时间复杂度是O(N^3)。

我建议使用LU分解来求解矩阵方程，这样就不需要计算行列式了。

最好使用一个帮助你的软件包。比如JAMA。

12x12或19x19的矩阵并不算大。通常需要解决具有数万或数十万自由度的问题。

下面是如何使用JAMA的实例代码。编译和运行时必须将JAMA JAR文件加入到CLASSPATH中：

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

以下是使用Apache Commons Math解决同样问题的代码，根据quant_dev的建议实现：

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

根据您的情况进行相应的调整。

- duffymo

---

谢谢你的回答。你是对的，不是指数级增长，我的意思是从矩阵大小为12x12开始计算行列式所需时间显著增加。我已经尝试了Jama，但我不知道如何使用（我在Java方面还很新）。我也会研究LU分解。谢谢。 - dedalo

我曾经在一个实验室工作，使用数百或数千个处理器，解决了拥有数十亿自由度的矩阵问题。 - davidtbernal

我曾经在一台Unix工作站上遇到一个问题，整整一个月都在解决它。那时候多处理器还不普及，我们每周必须对结果进行检查点以确保如果出现崩溃需要重新启动，不会浪费太多时间。 - duffymo

使用Math 3.6.1，您只需要调用MatrixUtils.inverse(MatrixUtils.createRealMatrix(matrix))。 - Spenhouet

我不知道那个方法是什么时候添加的。 这个答案发布已经超过七年了。 我不能预见 API 的变化。 我发布的内容应该仍然是正确的。它们都在我发布它们的那一天成功运行 - 超过七年了。 - duffymo

Answer 2

我建议使用Apache Commons Math 2.0。JAMA是一个已经停止维护的项目。ACM 2.0实际上是从JAMA中提取了线性代数，并进一步开发了它。

- quant_dev

---

不知道这个，quant_dev。谢谢你提供的信息。 - duffymo

Answer 3

la4j（Java线性代数）库支持矩阵求逆。以下是简要示例：

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

- Vladimir Kostyukov

Answer 4

对于那些寻求矩阵求逆（不考虑速度）的人，请查看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));
    }

}

- CoolMind

---

我猜测方法“minor”应该被称为“submatrix”？ - acala

@acala，当然，它是子矩阵，但我认为它也是次要的，因为我们删除了行和列，请参见https://en.wikipedia.org/wiki/Minor_(linear_algebra)。我很多年前学过，可能忘记了。 - CoolMind

但我认为一个余子式和一个子矩阵是不同的东西，余子式是子矩阵的行列式。请查看维基百科页面进行审核。 - acala

@acala，我同意你的看法，谢谢。在许多文章中，他们称小行列式为determinants。因此，我应该重命名这个方法。 - CoolMind

无法处理2*2矩阵。 - wddd

Answer 5

矩阵求逆是计算密集型的。正如duffymo所说，LU是一个很好的算法，还有其他变体（例如QR）。

不幸的是，你无法摆脱繁重的计算...如果你没有使用优化库，getSubmatrix方法可能是瓶颈。

此外，特殊的矩阵结构（带状、对称、对角线、稀疏性）在计算中考虑到会对性能产生巨大影响。你的表现可能会有所不同...

- Pablo Rodriguez

Answer 6

你永远不应该用这种方式计算逆矩阵。逆矩阵本身的计算应该避免，因为通常最好使用诸如LU之类的分解。

使用递归计算来计算行列式是一件数字上令人不快的事情。结果表明，更好的选择是使用LU分解来计算行列式。但是，如果你要费力地计算LU因子，为什么还要计算逆矩阵呢？通过计算LU因子，你已经完成了困难的工作。

一旦你有了LU因子，就可以使用它们进行前向和后向替换。

至于19x19矩阵是否很大，它甚至没有接近我认为的大。

- user85109

Answer 7

自从ACM库更新以来，这里提供符合最新矩阵求逆定义的实现。

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

- Soumya Kanti

---

这个错误意味着 Exception in thread "AWT-EventQueue-0" org.apache.commons.math3.linear.SingularMatrixException: matrix is singular 必须使用相同大小的矩阵吗？ - zygimantus

Answer 8

你计算行列式的算法确实是指数级的。基本问题在于你从定义开始计算，而直接定义会导致需要计算指数级数量的子行列式。在计算行列式或其逆矩阵之前，你真的需要先转换矩阵。 （我考虑过解释动态规划，但这个问题也无法通过动态规划来解决，因为子问题的数量也是指数级的。）

正如其他人建议的那样，LU分解是一个不错的选择。如果你是矩阵计算的新手，你可能还想看看高斯消元来计算行列式和逆矩阵，因为这可能更容易理解。

在矩阵求逆中要记住一件事就是数值稳定性，因为你正在处理浮点数。所有好的算法都包括行和/或列的置换来选择合适的主元，它们被称为“枢轴”。至少在高斯消元中，你希望在每一步中对列进行置换，以便选择绝对值最大的元素作为主元，因为这是最稳定的选择。

- JaakkoK

Answer 9

很难在它们的游戏中击败Matlab。他们对精度也非常谨慎。如果您有2.0和2.00001作为枢轴，请小心！您的答案可能会非常不精确。还要查看Python的实现（它在numpy/scipy中的某个地方...）

- Hamish Grubijan

Answer 10

你是否必须要有一个精确的解决方案？ 一个近似求解器（高斯-塞德尔）非常高效且易于实现，可能适合您，并且收敛速度非常快。19x19是相当小的矩阵。我认为我使用的高斯-塞德尔代码可以在眨眼之间解决128x128矩阵（但不要引用我说的话，这已经过了一段时间）。

- davidtbernal