Java обратный матричный расчет
Я пытаюсь вычислить обратную матрицу в Java.
Я следую за сопряженным методом (сначала вычисление сопряженной матрицы, затем транспонирование этой матрицы и, наконец, умножение ее на обратное значение определителя).
Работает, когда матрица не слишком большая. Я проверил, что для матриц размером до 12х12 результат предоставляется быстро. Однако, когда матрица больше 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;
}
Кто-нибудь знает, как более эффективно вычислить определитель большой матрицы? Если нет, кто-нибудь знает, как вычислить инверсию большой матрицы, используя другой алгоритм?
Спасибо
10 ответов
Экспоненциально? Нет, я считаю, что инверсия матрицы - это O(N^3).
Я бы рекомендовал использовать разложение LU для решения матричного уравнения. Вам не нужно искать определитель, когда вы его используете.
Еще лучше, посмотрите на пакет, чтобы помочь вам. Джама приходит на ум.
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 и развил ее дальше.
Для тех, кто ищет инверсию матрицы (не быстро), см. 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(minor(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(minor(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[][] minor(double[][] matrix, int row, int column) {
double[][] minor = 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)
minor[i < row ? i : i - 1][j < column ? j : j - 1] = matrix[i][j];
return minor;
}
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));
}
}
Библиотека la4j (Linear Algebra for 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
Поскольку библиотека 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);
}
}
Матричная инверсия в вычислительном отношении очень интенсивна. Как ответил Даффимо, LU - хороший алгоритм, и есть другие варианты (например, QR).
К сожалению, вы не можете избавиться от тяжелых вычислений... и, возможно, bottelneck - это метод getSubmatrix, если вы не используете оптимизированную библиотеку.
Кроме того, специальные матричные структуры (полосовая матрица, симметрия, диагональ, разреженность) оказывают большое влияние на производительность, если учитывать их в расчетах. Ваш пробег может варьироваться...
Вы НИКОГДА не хотите вычислять обратную матрицу таким образом. Хорошо, вычисление самого обратного следует избегать, так как почти всегда лучше использовать факторизацию, такую как LU.
Вычисление детерминанта с использованием рекурсивных вычислений является численно непристойным занятием. Оказывается, что лучший выбор - использовать факторизацию LU для вычисления детерминанта. Но если вы собираетесь потрудиться вычислить коэффициенты LU, то почему вы, возможно, захотите вычислить обратное? Вы уже выполнили сложную работу, рассчитав коэффициенты LU.
Если у вас есть коэффициенты LU, вы можете использовать их для замены вперед и назад.
Поскольку матрица 19x19 является большой, она даже не близка к тому, что я считаю большим.
Ваш алгоритм вычисления определителя действительно экспоненциальный. Основная проблема заключается в том, что вы вычисляете по определению, а прямое определение приводит к экспоненциальному количеству субдетерминантов для вычисления. Вам действительно нужно сначала преобразовать матрицу, прежде чем вычислять ее детерминант или инверсию. (Я думал об объяснении о динамическом программировании, но эта проблема не может быть решена динамическим программированием, поскольку число подзадач также является экспоненциальным.)
Разложение LU, как рекомендовано другими, является хорошим выбором. Если вы новичок в матричных вычислениях, вы также можете захотеть взглянуть на исключение Гаусса, чтобы вычислить детерминанты и инверсии, так как поначалу это будет немного легче понять.
И одну вещь, которую следует помнить при инверсии матриц, это числовая стабильность, поскольку вы имеете дело с числами с плавающей запятой. Весь хороший алгоритм включает в себя перестановки строк и / или столбцов для выбора подходящих опорных точек, как они называются. По крайней мере, при устранении по Гауссу вы хотите на каждом шаге переставлять столбцы, чтобы элемент, наибольший по абсолютной величине, был выбран в качестве основного, поскольку это самый стабильный выбор.
У вас должно быть точное решение? Приблизительный решатель ( Gauss-Seidel очень эффективен и прост в реализации), вероятно, подойдет вам и очень быстро сойдет. 19х19 - довольно маленькая матрица. Я думаю, что код Гаусса-Зейделя, который я использовал, мог бы решить матрицу 128x128 в мгновение ока (но не цитируйте меня об этом, это было давно).
Сложно обыграть Матлаба в их игре. Они также предостерегают о точности. Если у вас есть 2.0 и 2.00001 в качестве опорных точек - будьте осторожны! Ваш ответ может оказаться очень неточным. Кроме того, проверьте реализацию Python (где-то в numpy / scipy...)
С помощью ACM вы можете сделать это без права доступа:
RealMatrix A = new Array2DRowRealMatrix(ARR);
System.out.println(new LUDecomposition(A).getSolver().getInverse());