Calculating n x n determinant

Question

According to instructions and sample code from MSDN Magazine and comments from post here determinant of n x n...

I changed code to be like this:

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

Any comments or ideas on how to improve this code (especially this part):

 //----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

...are welcome.

Answer 1

Optimizing for performance and readability you may want to avoid the 3 executions of the same loop:

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

You may also want to loop backwards (try it and see if it works out better for you):

for (int k = result.GetLength(1)-1; k >= 0; k--)

And finally, something a little non-intuitive. Use double over float for accuracy and speed in .NET as your floats are probably converted to double anyway. CLR via C# covers this quirk but again, test it out as things may have changed.

Answer 2

Unused functions

Removing these 3 lines does not affect the final result.

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

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

Comments

I don't really see the point in // ----------------------------- on top of your methods. You can just type /// which automatically propose a neat documentation, per example, in the MatrixCreate function :

    /// <summary>
    /// 
    /// </summary>
    /// <param name="rows"></param>
    /// <param name="cols"></param>
    /// <returns></returns>
    private static float[,] MatrixCreate(int rows, int cols)

Hardcoded values

They are usually a bad idea (what happens if, per example, you use this value somewhere else, and decide to change it someday?):

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

Another option would be to have a threshold or tolerance value in your class.

Multiple accesses inside the same loop

Here you have to access result[i, j] every time in the second loop. This harms performance and makes this second loop unclear.

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

Could be written:

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

Function naming

ExtractLower and ExtractUpper are kind of misguiding, I would expect that ExtractLower(x)=ExtractUpper(Transpose(x)), but this is not the case as each function treats the diagonal in a different manner.

In order to improve readability, I would write a Transpose method and a SetDiagonalValuesTo(float[,] matrix, double value)

private static float[,] Transpose(float[,] matrix)
{
    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)
        {
            result[i, j] = matrix[j, i];
        }
    }
    return result;
}

private static float[,] ExtractLower(float[,] matrix)
{
    return ExtractUpper(Transpose(matrix));
}

private static float[,] SetDiag(float[,] matrix, float diagonalValue)
{
    int rows = matrix.GetLength(0);
    float[,] result = MatrixCreate(rows, rows);
    for (int i = 0; i < rows; ++i)
    {
        result[i, i] = diagonalValue;
    }
    return result;
}