Вычисление обратной матрицы Java

Экспоненциально? Нет, я считаю, что инверсия матрицы - это 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);
    }
}

Адаптируйте их к вашей ситуации.

Я бы рекомендовал использовать для этого Apache Commons Math 2.0. JAMA — мертвый проект. ACM 2.0 фактически взял линейную алгебру из JAMA и развил ее дальше.

Библиотека 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

Инверсия матриц требует больших вычислительных ресурсов. Как ответил duffymo, LU — хороший алгоритм, и есть и другие варианты (например, QR).

К сожалению, вы не можете избавиться от тяжелых вычислений... и, возможно, узким местом является метод getSubmatrix, если вы не используете оптимизированную библиотеку.

Кроме того, специальные матричные структуры (полосная матрица, симметрия, диагональность, разреженность) оказывают большое влияние на производительность, если учитывать их в расчетах. Ваш пробег может отличаться...

Вы НИКОГДА не хотите вычислять обратную матрицу таким образом. Хорошо, вычисления самого обратного следует избегать, так как почти всегда лучше использовать факторизацию, такую ​​как LU.

Вычисление определителя с помощью рекурсивных вычислений — численно непристойное занятие. Оказывается, лучшим выбором является использование факторизации LU для вычисления определителя. Но если вы собираетесь вычислять коэффициенты LU, то зачем вам вычислять обратное? Вы уже проделали сложную работу по вычислению коэффициентов LU.

Когда у вас есть факторы LU, вы можете использовать их для выполнения обратной и прямой замены.

Поскольку матрица 19x19 большая, она даже близко не соответствует тому, что я считаю большой.

Поскольку библиотека 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);
    }
}

Для тех, кто ищет инверсию матриц (не быстро), см. 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));
    }

}

Ваш алгоритм вычисления определителя действительно экспоненциальный. Основная проблема заключается в том, что вы вычисляете по определению, а прямое определение приводит к экспоненциальному количеству поддетерминантов для вычисления. Вам действительно нужно сначала преобразовать матрицу, прежде чем вычислять ее определитель или ее обратную. (Я думал объяснить о динамическом программировании, но эта проблема не может быть решена с помощью динамического программирования, так как количество подзадач также экспоненциально.)

Разложение LU, как рекомендуют другие, является хорошим выбором. Если вы новичок в вычислении матриц, вы также можете взглянуть на исключение Гаусса для вычисления определителей и инверсий, так как это может быть немного легче понять поначалу.

И одна вещь, которую следует помнить при инверсии матриц, — это числовая стабильность, поскольку вы имеете дело с числами с плавающей запятой. Все хорошие алгоритмы включают перестановки строк и/или столбцов для выбора подходящих опорных точек, как они называются. По крайней мере, при исключении Гаусса вы хотите на каждом шаге переставлять столбцы так, чтобы элемент с наибольшим абсолютным значением выбирался в качестве опорного, поскольку это самый стабильный выбор.

Трудно победить Matlab в их игре. Они также осторожно относятся к точности. Если у вас есть опорные значения 2.0 и 2.00001 — будьте осторожны! Ваш ответ может оказаться очень неточным. Кроме того, проверьте реализацию Python (где-то она находится в numpy/scipy...)

У вас должно быть точное решение? Приближенный решатель (Gauss-Seidel) очень эффективен и прост в реализации. ), скорее всего, сработает для вас и очень быстро сойдется. 19х19 это совсем маленькая матрица. Я думаю, что код Гаусса-Зейделя, который я использовал, мог решить матрицу 128x128 в мгновение ока (но не цитируйте меня, это было давно).

С помощью ACM вы можете сделать это без правой стороны:

RealMatrix A = new Array2DRowRealMatrix(ARR);
System.out.println(new LUDecomposition(A).getSolver().getInverse());