Вычисление определителя матрицы

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

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

Пока все хорошо, и это дает мне правильный результат. Теперь я хотел бы оптимизировать свой код, используя несколько потоков для вычисления этого определяющего значения. Я попытался распараллелить его, используя модель 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);
}

Однако после сравнения производительности я обнаружил, что производительность подхода Fork/Join очень плохая, и чем выше размерность матрицы, тем медленнее она становится (по сравнению с первым подходом). Где накладные расходы? Может ли кто-нибудь пролить свет на то, как это улучшить?


person all_by_grace    schedule 17.05.2013    source источник
comment
Прежде чем добавлять потоки, я бы прекратил выделение в цикле. Одним из вариантов может быть два параметра массива, определяющие, какие столбцы и строки следует вычислять вместо N.   -  person Stefan Haustein    schedule 17.05.2013
comment
Я бы посоветовал вам взглянуть на некоторые алгоритмы, разработанные для параллельной работы. Я не изучал ваш алгоритм, но по моему опыту можно найти много умных оптимизаций для распространенных проблем, выполнив поиск.   -  person Osama Javed    schedule 17.05.2013


Ответы (5)


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

Этот класс использует множество различных методов, чтобы сделать матрицу треугольной, а затем вычислить ее определитель. Его можно использовать для матриц больших размеров, таких как 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];
    }
} }

Этот класс получает от пользователя матрицу размера n x n, а затем вычисляет ее определитель. Он также показывает решение и окончательную треугольную матрицу.

 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;
}}
person Seyyed Mohsen Mousavi    schedule 12.05.2018
comment
Он не просил алгоритм/решение, чтобы найти определитель! Он ясно понимал многопоточность в своем вопросе! - person Yahya; 13.05.2018
comment
Этот алгоритм решает проблему производительности и настолько быстр в результатах, особенно для очень больших размеров, что я решил, что его стоит протестировать. - person Seyyed Mohsen Mousavi; 13.05.2018
comment
@SeyyedMohsenMousavi, чувак, этот алгоритм очень быстрый! Есть ли его описание где-нибудь в Интернете? - person Václav; 16.06.2021
comment
@ Вацлав, спасибо, что имеется в виду под описанием? Сам алгоритм или скорость? - person Seyyed Mohsen Mousavi; 19.06.2021
comment
@SeyyedMohsenMousavi сам алгоритм. Может быть, глава книги или статья. Эти треугольники - я никогда раньше не видел - person Václav; 19.06.2021
comment
@Václav Если вы ищете книгу, взгляните на это книгу у Адамса, раздел Матрица, и если вы просто хотите понять, на чем основан этот алгоритм, посмотрите этот видео на YouTube - person Seyyed Mohsen Mousavi; 20.06.2021

Основная причина того, что код ForkJoin работает медленнее, заключается в том, что он фактически сериализуется с добавлением некоторых накладных расходов на потоки. Чтобы извлечь выгоду из fork/join, вам нужно 1) сначала разветвить все экземпляры, а затем 2) дождаться результатов. Разделите свой цикл в «вычислении» на два цикла: один для разветвления (сохранение экземпляров ParallelDeterminants, скажем, в массиве), а другой для сбора результатов.

Кроме того, я предлагаю разветвляться только на самом внешнем уровне, а не на каком-либо из внутренних. Вы не хотите создавать потоки O (N ^ 2).

person Denis Dmitriev    schedule 30.08.2013

Существует новый метод вычисления определителя матрицы, о котором вы можете прочитать здесь

и я реализовал простую версию этого без каких-либо причудливых методов оптимизации или библиотеки в простой простой Java, и я протестировал методы, описанные ранее, и это было быстрее в среднем в 10 раз.

public class Test {
public static double[][] reduce(int row , int column , double[][] mat){
    int n=mat.length;
    double[][] res = new double[n- 1][n- 1];
    int r=0,c=0;
    for (int i = 0; i < n; i++) {
        c=0;
        if(i==row)
            continue;
        for (int j = 0; j < n; j++) {
            if(j==column)
                continue;
            res[r][c] = mat[i][j];

            c++;
        }
        r++;
    }
    return res;
}

public static double det(double mat[][]){
    int n = mat.length;
    if(n==1)
        return mat[0][0];
    if(n==2)
        return mat[0][0]*mat[1][1] - (mat[0][1]*mat[1][0]);
    //TODO : do reduce more efficiently
    double[][] m11 = reduce(0,0,mat);
    double[][] m1n = reduce(0,n-1,mat);
    double[][] mn1 = reduce(n-1 , 0 , mat);
    double[][] mnn = reduce(n-1,n-1,mat);
    double[][] m11nn = reduce(0,0,reduce(n-1,n-1,mat));
    return (det(m11)*det(mnn) - det(m1n)*det(mn1))/det(m11nn);
}

public static double[][] randomMatrix(int n , int range){
    double[][] mat = new double[n][n];
    for (int i=0; i<mat.length; i++) {
        for (int j=0; j<mat[i].length; j++) {
            mat[i][j] = (Math.random()*range);
        }
    }
    return mat;
}

public static void main(String[] args) {
    double[][] mat = randomMatrix(10,100);
    System.out.println(det(mat));
}
}

в случае с определителем m11nn есть небольшая ошибка, если он равен нулю, он взорвется, и вы должны это проверить. Я тестировал 100 случайных выборок, это случается редко, но я думаю, что об этом стоит упомянуть, а также использование лучшей схемы индексации также может повысить эффективность.

person sajjad    schedule 04.11.2020

Это часть моего класса Matrix, который использует переменную-член double[][] с именем data для хранения данных матрицы. Функция _determinant_recursivetask_impl() использует объект RecursiveTask<Double> с ForkJoinPool, чтобы попытаться использовать несколько потоков для вычислений.

Этот метод выполняется очень медленно по сравнению с матричными операциями для получения верхней/нижней треугольной матрицы. Например, попробуйте вычислить определитель матрицы 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];
    }
}
person Jaideep Heer    schedule 23.01.2020

person    schedule
comment
Не могли бы вы улучшить свой ответ, добавив комментарии к коду и, возможно, уточнить проблемы с производительностью, упомянутые в вопросе? - person Jindra Helcl; 08.12.2018