Java 计算矩阵行列式

Calculating matrix determinant

提问by all_by_grace

I am trying to calculate the determinant of a matrix (of any size), for self coding / interview practice. My first attempt is using recursion and that leads me to the following implementation:

我正在尝试计算矩阵（任何大小）的行列式，用于自我编码/面试练习。我的第一次尝试是使用递归，这导致我进行以下实现：

import java.util.Scanner.*;
public class Determinant {

    double A[][];
    double m[][];
    int N;
    int start;
    int last;

    public Determinant (double A[][], int N, int start, int last){
            this.A = A;
            this.N = N;
            this.start = start;
            this.last = last;
    }

    public double[][] generateSubArray (double A[][], int N, int j1){
            m = new double[N-1][];
            for (int k=0; k<(N-1); k++)
                    m[k] = new double[N-1];

            for (int i=1; i<N; i++){
                  int j2=0;
                  for (int j=0; j<N; j++){
                      if(j == j1)
                            continue;
                      m[i-1][j2] = A[i][j];
                      j2++;
                  }
              }
            return m;
    }
    /*
     * Calculate determinant recursively
     */
    public double determinant(double A[][], int N){
        double res;

        // Trivial 1x1 matrix
        if (N == 1) res = A[0][0];
        // Trivial 2x2 matrix
        else if (N == 2) res = A[0][0]*A[1][1] - A[1][0]*A[0][1];
        // NxN matrix
        else{
            res=0;
            for (int j1=0; j1<N; j1++){
                 m = generateSubArray (A, N, j1);
                 res += Math.pow(-1.0, 1.0+j1+1.0) * A[0][j1] * determinant(m, N-1);
            }
        }
        return res;
    }
}

So far it is all good and it gives me a correct result. Now I would like to optimise my code by making use of multiple threads to calculate this determinant value. I tried to parallelize it using the Java Fork/Join model. This is my approach:

到目前为止，一切都很好，它给了我一个正确的结果。现在我想通过使用多个线程来计算这个行列式值来优化我的代码。我尝试使用 Java Fork/Join 模型对其进行并行化。这是我的方法：

@Override
protected Double compute() {
     if (N < THRESHOLD) {
         result = computeDeterminant(A, N);
         return result;
     }

     for (int j1 = 0; j1 < N; j1++){
          m = generateSubArray (A, N, j1);
          ParallelDeterminants d = new ParallelDeterminants (m, N-1);
          d.fork();
          result += Math.pow(-1.0, 1.0+j1+1.0) * A[0][j1] * d.join();
     }

     return result;
}

public double computeDeterminant(double A[][], int N){
    double res;

    // Trivial 1x1 matrix
    if (N == 1) res = A[0][0];
    // Trivial 2x2 matrix
    else if (N == 2) res = A[0][0]*A[1][1] - A[1][0]*A[0][1];
    // NxN matrix
    else{
        res=0;
        for (int j1=0; j1<N; j1++){
             m = generateSubArray (A, N, j1);
             res += Math.pow(-1.0, 1.0+j1+1.0) * A[0][j1] * computeDeterminant(m, N-1);
        }
    }
    return res;
}

/*
 * Main function
 */
public static void main(String args[]){
    double res;
    ForkJoinPool pool = new ForkJoinPool();
    ParallelDeterminants d = new ParallelDeterminants();
    d.inputData();
    long starttime=System.nanoTime();
    res = pool.invoke (d);
    long EndTime=System.nanoTime();

    System.out.println("Seq Run = "+ (EndTime-starttime)/100000);
    System.out.println("the determinant valaue is  " + res);
}

However after comparing the performance, I found that the performance of the Fork/Join approach is very bad, and the higher the matrix dimension, the slower it becomes (as compared to the first approach). Where is the overhead? Can anyone shed a light on how to improve this?

但是经过性能对比，我发现 Fork/Join 方法的性能很差，矩阵维数越高，它变得越慢（与第一种方法相比）。开销在哪里？任何人都可以阐明如何改进这一点吗？

回答by Denis Dmitriev

The main reason the ForkJoin code is slower is that it's actually serialized with some thread overhead thrown in. To benefit from fork/join, you need to 1) fork all instances first, then 2) wait for the results. Split your loop in "compute" into two loops: one to fork (storing instances of ParallelDeterminants in, say, an array) and another to collect the results.

ForkJoin 代码较慢的主要原因是它实际上是序列化的，并且会带来一些线程开销。为了从 fork/join 中受益，您需要 1) 首先 fork 所有实例，然后 2) 等待结果。将“计算”中的循环拆分为两个循环：一个用于分叉（将 ParallelDeterminants 的实例存储在数组中），另一个用于收集结果。

Also, I suggest to only fork at the outermost level and not in any of the inner ones. You don't want to be creating O(N^2) threads.

另外，我建议只在最外层分叉，不要在任何内层分叉。您不想创建 O(N^2) 个线程。

回答by Seyyed Mohsen Mousavi

Using This class you can calculate the determinant of a matrix with any dimension

使用此类，您可以计算任意维度矩阵的行列式

This class uses many different methods to make the matrix triangularand then, calculates the determinant of it. It can be used for matrix of high dimension like 500 x 500 or even more. the bright side of the this class is that you can get the result in BigDecimalso there is no infinity and you'll have always the accurate answer. By the way, using many various methods and avoiding recursion resulted in much faster way with higher performance to the answer. hope it would be helpful.

这个类使用许多不同的方法使矩阵成为三角形，然后计算它的行​​列式。它可以用于高维矩阵，如 500 x 500 甚至更大。这门课的好处是你可以得到BigDecimal的结果，所以没有无穷大，你总是有准确的答案。顺便说一句，使用许多不同的方法并避免递归导致了更快的答案和更高的性能。希望它会有所帮助。

import java.math.BigDecimal;


public class DeterminantCalc {

private double[][] matrix;
private int sign = 1;


DeterminantCalc(double[][] matrix) {
    this.matrix = matrix;
}

public int getSign() {
    return sign;
}

public BigDecimal determinant() {

    BigDecimal deter;
    if (isUpperTriangular() || isLowerTriangular())
        deter = multiplyDiameter().multiply(BigDecimal.valueOf(sign));

    else {
        makeTriangular();
        deter = multiplyDiameter().multiply(BigDecimal.valueOf(sign));

    }
    return deter;
}


/*  receives a matrix and makes it triangular using allowed operations
    on columns and rows
*/
public void makeTriangular() {

    for (int j = 0; j < matrix.length; j++) {
        sortCol(j);
        for (int i = matrix.length - 1; i > j; i--) {
            if (matrix[i][j] == 0)
                continue;

            double x = matrix[i][j];
            double y = matrix[i - 1][j];
            multiplyRow(i, (-y / x));
            addRow(i, i - 1);
            multiplyRow(i, (-x / y));
        }
    }
}


public boolean isUpperTriangular() {

    if (matrix.length < 2)
        return false;

    for (int i = 0; i < matrix.length; i++) {
        for (int j = 0; j < i; j++) {
            if (matrix[i][j] != 0)
                return false;

        }

    }
    return true;
}


public boolean isLowerTriangular() {

    if (matrix.length < 2)
        return false;

    for (int j = 0; j < matrix.length; j++) {
        for (int i = 0; j > i; i++) {
            if (matrix[i][j] != 0)
                return false;

        }

    }
    return true;
}


public BigDecimal multiplyDiameter() {

    BigDecimal result = BigDecimal.ONE;
    for (int i = 0; i < matrix.length; i++) {
        for (int j = 0; j < matrix.length; j++) {
            if (i == j)
                result = result.multiply(BigDecimal.valueOf(matrix[i][j]));

        }

    }
    return result;
}


// when matrix[i][j] = 0 it makes it's value non-zero
public void makeNonZero(int rowPos, int colPos) {

    int len = matrix.length;

    outer:
    for (int i = 0; i < len; i++) {
        for (int j = 0; j < len; j++) {
            if (matrix[i][j] != 0) {
                if (i == rowPos) { // found "!= 0" in it's own row, so cols must be added
                    addCol(colPos, j);
                    break outer;

                }
                if (j == colPos) { // found "!= 0" in it's own col, so rows must be added
                    addRow(rowPos, i);
                    break outer;
                }
            }
        }
    }
}


//add row1 to row2 and store in row1
public void addRow(int row1, int row2) {

    for (int j = 0; j < matrix.length; j++)
        matrix[row1][j] += matrix[row2][j];
}


//add col1 to col2 and store in col1
public void addCol(int col1, int col2) {

    for (int i = 0; i < matrix.length; i++)
        matrix[i][col1] += matrix[i][col2];
}


//multiply the whole row by num
public void multiplyRow(int row, double num) {

    if (num < 0)
        sign *= -1;


    for (int j = 0; j < matrix.length; j++) {
        matrix[row][j] *= num;
    }
}


//multiply the whole column by num
public void multiplyCol(int col, double num) {

    if (num < 0)
        sign *= -1;

    for (int i = 0; i < matrix.length; i++)
        matrix[i][col] *= num;

}


// sort the cols from the biggest to the lowest value
public void sortCol(int col) {

    for (int i = matrix.length - 1; i >= col; i--) {
        for (int k = matrix.length - 1; k >= col; k--) {
            double tmp1 = matrix[i][col];
            double tmp2 = matrix[k][col];

            if (Math.abs(tmp1) < Math.abs(tmp2))
                replaceRow(i, k);
        }
    }
}


//replace row1 with row2
public void replaceRow(int row1, int row2) {

    if (row1 != row2)
        sign *= -1;

    double[] tempRow = new double[matrix.length];

    for (int j = 0; j < matrix.length; j++) {
        tempRow[j] = matrix[row1][j];
        matrix[row1][j] = matrix[row2][j];
        matrix[row2][j] = tempRow[j];
    }
}


//replace col1 with col2
public void replaceCol(int col1, int col2) {

    if (col1 != col2)
        sign *= -1;

    System.out.printf("replace col%d with col%d, sign = %d%n", col1, col2, sign);
    double[][] tempCol = new double[matrix.length][1];

    for (int i = 0; i < matrix.length; i++) {
        tempCol[i][0] = matrix[i][col1];
        matrix[i][col1] = matrix[i][col2];
        matrix[i][col2] = tempCol[i][0];
    }
} }

This Class Receives a matrix of n x n from the user then calculates it's determinant. It also shows the solution and the final triangular matrix.

这个类从用户那里接收一个 nxn 的矩阵，然后计算它的行​​列式。它还显示了解决方案和最终的三角矩阵。

 import java.math.BigDecimal;
 import java.text.NumberFormat;
 import java.util.Scanner;


public class DeterminantTest {

public static void main(String[] args) {

    String determinant;

    //generating random numbers
    /*int len = 300;
    SecureRandom random = new SecureRandom();
    double[][] matrix = new double[len][len];

    for (int i = 0; i < len; i++) {
        for (int j = 0; j < len; j++) {
            matrix[i][j] = random.nextInt(500);
            System.out.printf("%15.2f", matrix[i][j]);
        }
    }
    System.out.println();*/

    /*double[][] matrix = {
        {1, 5, 2, -2, 3, 2, 5, 1, 0, 5},
        {4, 6, 0, -2, -2, 0, 1, 1, -2, 1},
        {0, 5, 1, 0, 1, -5, -9, 0, 4, 1},
        {2, 3, 5, -1, 2, 2, 0, 4, 5, -1},
        {1, 0, 3, -1, 5, 1, 0, 2, 0, 2},
        {1, 1, 0, -2, 5, 1, 2, 1, 1, 6},
        {1, 0, 1, -1, 1, 1, 0, 1, 1, 1},
        {1, 5, 5, 0, 3, 5, 5, 0, 0, 6},
        {1, -5, 2, -2, 3, 2, 5, 1, 1, 5},
        {1, 5, -2, -2, 3, 1, 5, 0, 0, 1}
    };
    */

    double[][] matrix = menu();

    DeterminantCalc deter = new DeterminantCalc(matrix);

    BigDecimal det = deter.determinant();

    determinant = NumberFormat.getInstance().format(det);

    for (int i = 0; i < matrix.length; i++) {
        for (int j = 0; j < matrix.length; j++) {
            System.out.printf("%15.2f", matrix[i][j]);
        }
        System.out.println();
    }

    System.out.println();
    System.out.printf("%s%s%n", "Determinant: ", determinant);
    System.out.printf("%s%d", "sign: ", deter.getSign());

}


public static double[][] menu() {

    Scanner scanner = new Scanner(System.in);

    System.out.print("Matrix Dimension: ");
    int dim = scanner.nextInt();

    double[][] inputMatrix = new double[dim][dim];

    System.out.println("Set the Matrix: ");
    for (int i = 0; i < dim; i++) {
        System.out.printf("%5s%d%n", "row", i + 1);
        for (int j = 0; j < dim; j++) {

            System.out.printf("M[%d][%d] = ", i + 1, j + 1);
            inputMatrix[i][j] = scanner.nextDouble();
        }
        System.out.println();
    }
    scanner.close();

    return inputMatrix;
}}

回答by M.Sey

int det(int[][] mat) {
    if (mat.length == 1)
        return mat[0][0];
    if (mat.length == 2)
        return mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
    int sum = 0, sign = 1;
    int newN = mat.length - 1;
    int[][] temp = new int[newN][newN];
    for (int t = 0; t < newN; t++) {
        int q = 0;
        for (int i = 0; i < newN; i++) {
            for (int j = 0; j < newN; j++) {
                temp[i][j] = mat[1 + i][q + j];
            }
            if (q == i)
                q = 1;
        }
        sum += sign * mat[0][t] * det(temp);
        sign *= -1;
    }
    return sum;
}

回答by Jaideep Heer

This is a part of my Matrix class which uses a double[][]member variable called datato store the matrix data. The _determinant_recursivetask_impl()function uses a RecursiveTask<Double>object with the ForkJoinPoolto try to use multiple threads for calculation.

这是我的 Matrix 类的一部分，它使用一个double[][]称为data存储矩阵数据的成员变量。该_determinant_recursivetask_impl()函数使用带有 的RecursiveTask<Double>对象ForkJoinPool来尝试使用多个线程进行计算。

This method performs very slow compared to matrix operations to get an upper/lower triangular matrix. Try to compute the determinant of a 13x13 matrix for example.

与获取上/下三角矩阵的矩阵运算相比，此方法执行速度非常慢。例如，尝试计算 13x13 矩阵的行列式。

public class Matrix
{
    // Dimensions
    private final int I,J;
    private final double[][] data;
    private Double determinant = null;
    static class MatrixEntry
    {
        public final int I,J;
        public final double value;
        private MatrixEntry(int i, int j, double value) {
            I = i;
            J = j;
            this.value = value;
        }
    }

    /**
     * Calculates determinant of this Matrix recursively and caches it for future use.
     * @return determinant
     */
    public double determinant()
    {
        if(I!=J)
            throw new IllegalStateException(String.format("Can't calculate determinant of (%d,%d) matrix, not a square matrix.", I,J));
        if(determinant==null)
            determinant = _determinant_recursivetask_impl(this);
        return determinant;
    }
    private static double _determinant_recursivetask_impl(Matrix m)
    {
        class determinant_recurse extends RecursiveTask<Double>
        {
            private final Matrix m;
            determinant_recurse(Matrix m) {
                this.m = m;
            }

            @Override
            protected Double compute() {
                // Base cases
                if(m.I==1 && m.J==1)
                    return m.data[0][0];
                else if(m.I==2 && m.J==2)
                    return m.data[0][0]*m.data[1][1] - m.data[0][1]*m.data[1][0];
                else
                {
                    determinant_recurse[] tasks = new determinant_recurse[m.I];
                    for (int i = 0; i <m.I ; i++) {
                        tasks[i] = new determinant_recurse(m.getSubmatrix(0, i));
                    }
                    for (int i = 1; i <m.I ; i++) {
                        tasks[i].fork();
                    }
                    double ret = m.data[0][0]*tasks[0].compute();
                    for (int i = 1; i < m.I; i++) {
                        if(i%2==0)
                            ret += m.data[0][i]*tasks[i].join();
                        else
                            ret -= m.data[0][i]*tasks[i].join();
                    }
                    return ret;
                }
            }
        }
        return ForkJoinPool.commonPool().invoke(new determinant_recurse(m));
    }

    private static void _map_impl(Matrix ret, Function<Matrix.MatrixEntry, Double> operator)
    {
        for (int i = 0; i <ret.I ; i++) {
            for (int j = 0; j <ret.J ; j++) {
                ret.data[i][j] = operator.apply(new Matrix.MatrixEntry(i,j,ret.data[i][j]));
            }
        }
    }
    /**
     * Returns a new Matrix that is sub-matrix without the given row and column.
     * @param removeI row to remove
     * @param removeJ col. to remove
     * @return new Matrix.
     */
    public Matrix getSubmatrix(int removeI, int removeJ)
    {
        if(removeI<0 || removeJ<0 || removeI>=this.I || removeJ>=this.J)
            throw new IllegalArgumentException(String.format("Invalid element position (%d,%d) for matrix(%d,%d).", removeI,removeJ,this.I,this.J));
        Matrix m = new Matrix(this.I-1, this.J-1);
        _map_impl(m, (e)->{
            int i = e.I, j = e.J;
            if(e.I >= removeI) ++i;
            if(e.J >= removeJ) ++j;
            return this.data[i][j];
        });
        return m;
    }
    // Constructors
    public Matrix(int i, int j) {
        if(i<1 || j<1)
            throw new IllegalArgumentException(String.format("Invalid array dimensions: (%d,%d)", i, j));
        I = i;
        J = j;
        data = new double[I][J];
    }
}