How to calculate matrix determinant? n*n or just 5*5

Question

everyone. I need to find matrix n*n (or 5*5) determinant. I have a function translated from Pascal, but there's INDEX OUT OF RANGE EXCEPTION. Could somebody help me?

Here's my code:

public static double DET(double[,] a, int n)
    {
        int i, j, k;
        double det = 0;
        for (i = 0; i < n - 1; i++)
        {   
            for (j = i + 1; j < n + 1; j++)
            {
                det = a[j, i] / a[i, i];
                for (k = i; k < n; k++)
                    a[j, k] = a[j, k] - det * a[i, k]; // Here's exception
            }
        }
        det = 1;
        for (i = 0; i < n; i++)
            det = det * a[i, i];
            return det;
    }

Thanx for any help.

Answer 1

Working solution for calculating n * n determinant looks like:

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

Answer 2

for (j = i + 1; j < n + 1; j++)

Last J value will be bigger than array size. So you must to recheck array sizes and all how was all indexes translated from pascal.

Answer 3

Why bother with a translation when you can download working C# code

• http://www.codeproject.com/KB/cs/CSML.aspx
• http://www.codeproject.com/KB/recipes/AdvancedMatrixLibrary.aspx
• http://www.codeproject.com/KB/recipes/psdotnetmatrix.aspx

Answer 4

For bigger matrices, you might want to run the Bareiss algorithm to calculate the determinant:

• http://en.wikipedia.org/wiki/Bareiss_algorithm

Answer 5

I think that this algorithm is not good, at least for calculation of 5x5 matrices. Even If we correct this

for (j = i + 1; j < n + 1; j++)

to be like this

for (j = i + 1; j < n; j++)

And then write a complete code such as:

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

                for (k = i; k < n; k++)
                {
                    value = this.readElement(j, k) - det * this.readElement(i, k);

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

Result is: determinant is: -8 determinant is: -142 determinant is: -NaN

NaN occurs because of division by zero (I debugged it) It could be possible that for some very specific input this works OK, but in general case this is NOT a good algorithm.

So, it works for 3x3 and 4x4 but NOT for 5x5

I wrote this to anyone who might come across this question to avoid losing few hours in trying to implement or fix something that has wrong algorithm in the first place.