определитель матрицы n*n

Я пытаюсь реализовать код и алгоритмы, найденные здесь:

определитель матрицы и здесь: Как вычислить определитель матрицы? n*n или просто 5*5

Но я застрял с этим.

Мой первый вопрос заключается в том, какое правило на самом деле используют эти алгоритмы (поскольку в математике, очевидно, есть несколько правил, по которым кто-то может вычислить определитель), поэтому я хотел бы сначала проверить, правильно ли применяется алгоритм.

Мой второй вопрос, что я сделал не так (я имею в виду с реализацией) или что не так с самим алгоритмом, так как похоже, что для 3x3 и 4x4 он работает нормально, а для 5x5 выдает NaN. Результаты проверяются с помощью нескольких онлайн-калькуляторов определителя матриц, и они в порядке, за исключением 5x5.

Это мой код:

using System;

public class Matrix
{
    private int row_matrix; //number of rows for matrix
    private int column_matrix; //number of colums for matrix
    private double[,] matrix; //holds values of matrix itself

    //create r*c matrix and fill it with data passed to this constructor
    public Matrix(double[,] double_array)
    {
        matrix = double_array;
        row_matrix = matrix.GetLength(0);
        column_matrix = matrix.GetLength(1);
        Console.WriteLine("Contructor which sets matrix size {0}*{1} and fill it with initial data executed.", row_matrix, column_matrix);
    }

    //returns total number of rows
    public int countRows()
    {
        return row_matrix;
    }

    //returns total number of columns
    public int countColumns()
    {
        return column_matrix;
    }

    //returns value of an element for a given row and column of matrix
    public double readElement(int row, int column)
    {
        return matrix[row, column];
    }

    //sets value of an element for a given row and column of matrix
    public void setElement(double value, int row, int column)
    {
        matrix[row, column] = value;
    }

    public double deterMatrix()
    {
        double det = 0;
        double value = 0;
        int i, j, k;

        i = row_matrix;
        j = column_matrix;
        int n = i;

        if (i != j)
        {
            Console.WriteLine("determinant can be calculated only for sqaure matrix!");
            return det;
        }

        for (i = 0; i < n - 1; i++)
        {
            for (j = i + 1; j < n; j++)
            {
                det = (this.readElement(j, i) / this.readElement(i, i));

                //Console.WriteLine("readElement(j, i): " + this.readElement(j, i));
                //Console.WriteLine("readElement(i, i): " + this.readElement(i, i));
                //Console.WriteLine("det is" + det);
                for (k = i; k < n; k++)
                {
                    value = this.readElement(j, k) - det * this.readElement(i, k);

                    //Console.WriteLine("Set value is:" + value);
                    this.setElement(value, j, k);
                }
            }
        }
        det = 1;
        for (i = 0; i < n; i++)
            det = det * this.readElement(i, i);

        return det;
    }
}

internal class Program
{
    private static void Main(string[] args)
    {
        Matrix mat03 = new Matrix(new[,]
        {
            {1.0, 2.0, -1.0},
            {-2.0, -5.0, -1.0},
            {1.0, -1.0, -2.0},
        });

        Matrix mat04 = new Matrix(new[,]
        {
            {1.0, 2.0, 1.0, 3.0},
            {-2.0, -5.0, -2.0, 1.0},
            {1.0, -1.0, -3.0, 2.0},
            {4.0, -1.0, -3.0, 1.0},
        });

        Matrix mat05 = new Matrix(new[,]
        {
            {1.0, 2.0, 1.0, 2.0, 3.0},
            {2.0, 1.0, 2.0, 2.0, 1.0},
            {3.0, 1.0, 3.0, 1.0, 2.0},
            {1.0, 2.0, 4.0, 3.0, 2.0},
            {2.0, 2.0, 1.0, 2.0, 1.0},
        });

        double determinant = mat03.deterMatrix();
        Console.WriteLine("determinant is: {0}", determinant);

        determinant = mat04.deterMatrix();
        Console.WriteLine("determinant is: {0}", determinant);

        determinant = mat05.deterMatrix();
        Console.WriteLine("determinant is: {0}", determinant);
    }
}

person Nenad Bulatovic    schedule 20.03.2013    source источник
comment
Нэн, вероятно, из-за 0.0/0.0.   -  person leppie    schedule 20.03.2013
comment
Я знаю, что это из-за этого или даже из-за n/0.0, но вопрос в том, почему алгоритм вообще попадает в такую ​​ситуацию. Должно быть, это рабочий алгоритм, предложенный в других сообщениях (я упоминал их в самом начале моего вопроса, кроме того, он работает для 3x3 и 4x4.   -  person Nenad Bulatovic    schedule 20.03.2013
comment
Алгоритм, по-видимому, использует исключение Гаусса. Но предположение, что деление справедливо, очевидно, ложно. Мой совет — правильно внедрить исключение Гаусса, если это то, что вы хотите сделать. Вероятно, этот алгоритм верен только при определенных предположениях; может быть, что система уравнений допускает решения? В любом случае, я бы переписал его с нуля на вашем месте.   -  person Eric Lippert    schedule 20.03.2013
comment
Я согласен. Однако это означает, что решения, предложенные в двух других решениях по ссылкам, приведенным выше, на самом деле являются ЛОЖНЫМИ решениями, по крайней мере, для общего случая. Поправьте меня, если я ошибаюсь?   -  person Nenad Bulatovic    schedule 20.03.2013
comment
@nenad: я понятия не имею, подходит ли код, который вы случайно взяли из Интернета, для какой-либо конкретной цели или нет.   -  person Eric Lippert    schedule 21.03.2013


Ответы (2)


Зачем изобретать велосипед? Хорошо известный метод для получения определенного И для инвертирования матрицы состоит в том, чтобы выполнить LU разложение. Фактически, журнал MSDN недавно опубликовал сообщение о том, как сделать это с помощью C#, здесь http://msdn.microsoft.com/en-us/magazine/jj863137.aspx.

Пример кода

Определитель матрицы

Имея в руках метод разложения матрицы, легко написать код для вычисления определителя матрицы:

static double MatrixDeterminant(double[][] matrix)
{
  int[] perm;
  int toggle;
  double[][] lum = MatrixDecompose(matrix, out perm, out toggle);
  if (lum == null)
    throw new Exception("Unable to compute MatrixDeterminant");
  double result = toggle;
  for (int i = 0; i < lum.Length; ++i)
    result *= lum[i][i];
  return result;
}

Определитель вычисляется из произведения диагоналей на разложенной матрице с проверкой знака. Прочтите статью для более подробной информации.

Обратите внимание, что они используют зубчатый массив для матрицы, но вы можете заменить свое собственное хранилище матрицы, переводя lum[i][j] в lum[i,j].

person John Alexiou    schedule 20.03.2013
comment
Хорошо, я собираюсь реализовать это и сообщить вам о результате. - person Nenad Bulatovic; 20.03.2013
comment
Ты был прав, это правильный способ сделать это. Как вы указали, в этой статье используется подход статического метода и дизайн массива массивов для матриц. Из-за других требований в моем полном коде я ДОЛЖЕН использовать многомерный массив С#. До сих пор я делал это довольно хорошо, вот часть, которую я не понимаю: double[] rowPtr = result[pRow]; результат[pRow] = результат[j]; результат[j] = строкаPtr; Я не знаю, как скрыть это, начиная с этого: double[] rowPtr = result[pRow]; Даже если я распечатаю с помощью Console.WriteLine(result[pRow]); пишет: System.Double[] - person Nenad Bulatovic; 22.03.2013
comment
rowPtr[] — это одномерный массив, содержащий все элементы одной строки матрицы. Вам понадобится цикл по столбцам для передачи значений из 2D-массива в 1D-массив rowPtr. for(j=0; j<M; j++) { rowPtr[j] = result[i,j]; } - person John Alexiou; 22.03.2013

@ ja72 Большое спасибо за ваши инструкции. Окончательное решение для вычисления любого определителя n*n выглядит так:

using System;

internal class MatrixDecompositionProgram
{
    private static void Main(string[] args)
    {
        float[,] m = MatrixCreate(4, 4);
        m[0, 0] = 3.0f; m[0, 1] = 7.0f; m[0, 2] = 2.0f; m[0, 3] = 5.0f;
        m[1, 0] = 1.0f; m[1, 1] = 8.0f; m[1, 2] = 4.0f; m[1, 3] = 2.0f;
        m[2, 0] = 2.0f; m[2, 1] = 1.0f; m[2, 2] = 9.0f; m[2, 3] = 3.0f;
        m[3, 0] = 5.0f; m[3, 1] = 4.0f; m[3, 2] = 7.0f; m[3, 3] = 1.0f;

        int[] perm;
        int toggle;

        float[,] luMatrix = MatrixDecompose(m, out perm, out toggle);

        float[,] lower = ExtractLower(luMatrix);
        float[,] upper = ExtractUpper(luMatrix);

        float det = MatrixDeterminant(m);

        Console.WriteLine("Determinant of m computed via decomposition = " + det.ToString("F1"));
    }

    // --------------------------------------------------------------------------------------------------------------
    private static float[,] MatrixCreate(int rows, int cols)
    {
        // allocates/creates a matrix initialized to all 0.0. assume rows and cols > 0
        // do error checking here
        float[,] result = new float[rows, cols];
        return result;
    }

    // --------------------------------------------------------------------------------------------------------------
    private static float[,] MatrixDecompose(float[,] matrix, out int[] perm, out int toggle)
    {
        // Doolittle LUP decomposition with partial pivoting.
        // rerturns: result is L (with 1s on diagonal) and U; perm holds row permutations; toggle is +1 or -1 (even or odd)
        int rows = matrix.GetLength(0);
        int cols = matrix.GetLength(1);

        //Check if matrix is square
        if (rows != cols)
            throw new Exception("Attempt to MatrixDecompose a non-square mattrix");

        float[,] result = MatrixDuplicate(matrix); // make a copy of the input matrix

        perm = new int[rows]; // set up row permutation result
        for (int i = 0; i < rows; ++i) { perm[i] = i; } // i are rows counter

        toggle = 1; // toggle tracks row swaps. +1 -> even, -1 -> odd. used by MatrixDeterminant

        for (int j = 0; j < rows - 1; ++j) // each column, j is counter for coulmns
        {
            float colMax = Math.Abs(result[j, j]); // find largest value in col j
            int pRow = j;
            for (int i = j + 1; i < rows; ++i)
            {
                if (result[i, j] > colMax)
                {
                    colMax = result[i, j];
                    pRow = i;
                }
            }

            if (pRow != j) // if largest value not on pivot, swap rows
            {
                float[] rowPtr = new float[result.GetLength(1)];

                //in order to preserve value of j new variable k for counter is declared
                //rowPtr[] is a 1D array that contains all the elements on a single row of the matrix
                //there has to be a loop over the columns to transfer the values
                //from the 2D array to the 1D rowPtr array.
                //----tranfer 2D array to 1D array BEGIN

                for (int k = 0; k < result.GetLength(1); k++)
                {
                    rowPtr[k] = result[pRow, k];
                }

                for (int k = 0; k < result.GetLength(1); k++)
                {
                    result[pRow, k] = result[j, k];
                }

                for (int k = 0; k < result.GetLength(1); k++)
                {
                    result[j, k] = rowPtr[k];
                }

                //----tranfer 2D array to 1D array END

                int tmp = perm[pRow]; // and swap perm info
                perm[pRow] = perm[j];
                perm[j] = tmp;

                toggle = -toggle; // adjust the row-swap toggle
            }

            if (Math.Abs(result[j, j]) < 1.0E-20) // if diagonal after swap is zero . . .
                return null; // consider a throw

            for (int i = j + 1; i < rows; ++i)
            {
                result[i, j] /= result[j, j];
                for (int k = j + 1; k < rows; ++k)
                {
                    result[i, k] -= result[i, j] * result[j, k];
                }
            }
        } // main j column loop

        return result;
    } // MatrixDecompose

    // --------------------------------------------------------------------------------------------------------------
    private static float MatrixDeterminant(float[,] matrix)
    {
        int[] perm;
        int toggle;
        float[,] lum = MatrixDecompose(matrix, out perm, out toggle);
        if (lum == null)
            throw new Exception("Unable to compute MatrixDeterminant");
        float result = toggle;
        for (int i = 0; i < lum.GetLength(0); ++i)
            result *= lum[i, i];

        return result;
    }

    // --------------------------------------------------------------------------------------------------------------
    private static float[,] MatrixDuplicate(float[,] matrix)
    {
        // allocates/creates a duplicate of a matrix. assumes matrix is not null.
        float[,] result = MatrixCreate(matrix.GetLength(0), matrix.GetLength(1));
        for (int i = 0; i < matrix.GetLength(0); ++i) // copy the values
            for (int j = 0; j < matrix.GetLength(1); ++j)
                result[i, j] = matrix[i, j];
        return result;
    }

    // --------------------------------------------------------------------------------------------------------------
    private static float[,] ExtractLower(float[,] matrix)
    {
        // lower part of a Doolittle decomposition (1.0s on diagonal, 0.0s in upper)
        int rows = matrix.GetLength(0); int cols = matrix.GetLength(1);
        float[,] result = MatrixCreate(rows, cols);
        for (int i = 0; i < rows; ++i)
        {
            for (int j = 0; j < cols; ++j)
            {
                if (i == j)
                    result[i, j] = 1.0f;
                else if (i > j)
                    result[i, j] = matrix[i, j];
            }
        }
        return result;
    }

    // --------------------------------------------------------------------------------------------------------------
    private static float[,] ExtractUpper(float[,] matrix)
    {
        // upper part of a Doolittle decomposition (0.0s in the strictly lower part)
        int rows = matrix.GetLength(0); int cols = matrix.GetLength(1);
        float[,] result = MatrixCreate(rows, cols);
        for (int i = 0; i < rows; ++i)
        {
            for (int j = 0; j < cols; ++j)
            {
                if (i <= j)
                    result[i, j] = matrix[i, j];
            }
        }
        return result;
    }
}
person Nenad Bulatovic    schedule 25.03.2013