Determinant of a matrix of doubles

Question

I have been working with some matrix functions and need the determinant. I have looked at several examples, but they all seem to use a hard coded size for the matrix, or calculate for a matrix of int, etc.

My data is stored in a vector<vector<double>> and can be of any size. This is a typical square matrix of a matrix and its transpose MT*M.

This is my current code and is based in part on this post.

#include <vector>
#include <iostream>
#include <cstdlib>

// add namespace
using namespace std;

// declare, dynamically size, and return a 2D array of double as pointer to pointer
double** matrix_allocate_size(const unsigned int rows, const unsigned int columns) {

   // declare a pointer to a 2 dimensiona array of double
   // and allocate the number of rows
   double **matrix = (double **) calloc( rows, sizeof(double) );

   // allocate space for each column of each row
   // using calloc here initalizes array elements to 0
   for(unsigned int i=0; i<columns; i++) {
      matrix[i] = (double *) calloc(columns, sizeof(double));
   }

return matrix;
};


// calculate the determinant of a square matrix
// this function is called recursively
double matrix_determinant(unsigned int size, double **matrix) {

   // 2D array version of cofactor matrix
   // declare pointer to pointer and size by function call
   double **matrix_copy = matrix_allocate_size(size, size);

   // return value
   double determinant = 0.0;

   // sign for each product
   double sign = 1.0;

   // loop iterator
   unsigned int position = 0;

   // current location in array_matrix
   unsigned int matrix_row = 0, matrix_column = 0;

   // current location in matrix copy
   unsigned int copy_row = 0, copy_column = 0;

   // for matrix with one element
   if(size == 1) { return (matrix[0][0]); }

   // for a matrix with 4 elements, use the formula
   // this is essentially length*width = area
   else if(size == 2) { return ( (matrix[0][0]*matrix[1][1])-(matrix[1][0]*matrix[0][1]) ); }

   // for a matrix with 9 or more elements
   else {

      // reinitialize local determinant value
      determinant = 0.0;

      // loop on size of current array_matrix
      for(position=0; position<size; position++) {

         // zero out current matrix position
         copy_row = 0;  copy_column = 0;

         // create determinant products by skipping each row and column in sequence
         // iterate through array_matrix up to size,size
         // size decrements on each recursive call
         for(matrix_row=0; matrix_row<size; matrix_row++) {
            for(matrix_column=0; matrix_column<size; matrix_column++) {

               // skip current row and column when making copy
               if(matrix_row != 0 && matrix_column != position) {

                  // copy matrix elements not part of current row and column
                  matrix_copy[copy_row][copy_column] = matrix[matrix_row][matrix_column];

                  // increment column location in matrix copy up to size-2
                  if( copy_column<(size-2) ) { copy_column++; }
                  // move to next row in matrix copy, reset column number to 0
                  else { copy_column=0;  copy_row++; }

               }

            }
         }

         // calculate current element of determinant and add to sum
         // pass matrix copy in all subsequent recursive calls
         determinant += sign * ( matrix[0][position] * matrix_determinant(size-1, matrix_copy) );

         // reverse sign for next iteration
         sign = -1 * sign;

      } // for(position=0; position<size; position++) endbrace

   } // else endbrace

return (determinant);
};


// re-cast data from vector<vector<double>> to 2D array
// creates a matrix as a 2D array ot pointer to pointer
// populates the matrix, and passed it to a function to calculate the determinant
// returns the calculated determinant
double matrix_determinant_caller(const vector< vector<double> >& M1) {

   // return value
   double M1_determinant = 0.0;

   // store size of passed vector<vector<double>> matrix
   unsigned int M1_dimension = M1.size();

   // call function to size array version of matrix[rows][columns]
   // this is a square matrix, so use the same dimension for rows and columns
   double **array_matrix = matrix_allocate_size(M1_dimension, M1_dimension);

   // copy M1 to a_matrix
   for(unsigned int i=0; i<M1_dimension; i++) {
      for(unsigned int j=0; j<M1[0].size(); j++) { array_matrix[i][j] = M1[i][j]; }
   }

   // call function to get determinant
   M1_determinant = matrix_determinant(M1_dimension, array_matrix);

return M1_determinant;
};


// program entry point
int main() {

   // vector version of matrix
   vector< vector<double> > matrix_vector_version;

   // number of matrix rows and matrix columns
   unsigned int num_rows = 4, num_columns = 4;

   // temp matrix row for insert
   vector<double> temp_row;

   // size temp row vector
   temp_row.resize(num_columns);

   // add rows and columns to matrix
   for(unsigned int i=0; i<num_rows; i++) { matrix_vector_version.push_back(temp_row); }

   // add test data to matrix row 0
   matrix_vector_version[0][0] = 1248.14;
   matrix_vector_version[0][1] = 1408.68;
   matrix_vector_version[0][2] = -828.282;
   matrix_vector_version[0][3] =  -53.0927;

   // add test data to matrix row 1
   matrix_vector_version[1][0] = 1408.68 ;
   matrix_vector_version[1][1] = 1623.81;
   matrix_vector_version[1][2] = -952.41;
   matrix_vector_version[1][3] =  -67.2946 ;

   // add test data to matrix row 2
   matrix_vector_version[2][0] = -828.282;
   matrix_vector_version[2][1] = -952.41;
   matrix_vector_version[2][2] =  559.421;
   matrix_vector_version[2][3] =   38.0384;

   // add test data to matrix row 3
   matrix_vector_version[3][0] =  -53.0927;
   matrix_vector_version[3][1] =  -67.2946;
   matrix_vector_version[3][2] =   38.0384;
   matrix_vector_version[3][3] =    5.46328;

   // call function to calculate determinant
   double determinant = matrix_determinant_caller(matrix_vector_version);

   // print determinant
   // should be 44.1164 for this data
   cout << "determinant: " << determinant << endl;
   cout << endl;

return 0;
} // main EOF brace

This code compiles and runs and gives the correct value for the determinant of 44.1164 which is confirmed by an online calculator.

It would be nice to have some input on ways to improve the code. In particular, I think there is no exception handling to speak of and I am sure that performance and stability can be improved. The test matrix here is a 4x4 but I will need to use this on matrices that are at least 50x50. I guess in that case I might want to use a long double for the return value.

Current run time for this data is,

real    0m0.000s
user    0m0.030s
sys     0m0.000s

I can post a more realistic version of the data, or a version that reads an input file if that would be helpful.

This is more or less written in c++98. I don't object to more recent versions, but at times I work on c and c++ that is married to very old F77 code and I can't always get the newer gfortran compiler versions to be happy with my very old fortran code. I generally have good luck if I stick with c++98 standards and older c++ compilers.

Answer 1

Are you sure you need to do this?

Be aware that determinants for any reasonable sized matrix are very expensive to calculate, so if you're actually doing something else (linear solve, eigen-solve etc) there is probably a better way without calculating the determinant. And there are many external libraries which will solve this stuff for you.

In a real code review, I would be asking why you are implementing this yourself.

Getting an idea of the complexity

The best known algorithms are approx O(n^2.4), and the best actually implemented algorithms are O(n³) (as discussed in this SO question). That would put a rank 50 matrix determinant about 4600x slower than a 3x3. So if you are going to need determinants of large matrices, make sure your method will permit that to calculate in an acceptable time frame.

This method, if I understand it correctly, calculates the determinants of n minor matrices each of rank n-1. And that happens recursively. So y(n)=n.y(n-1) -> y(n)/y(n-1) = n. This is a method that will be O(n!)... at 50x50 it will take of order 10⁶³ times longer than the 3x3.

Critique of the code itself

Aside from the questions over the method, we can look at the code itself.

Data structure

An array of pointers to arrays of doubles has some issues in terms of performance; every access of a matrix element requires 2 dereferences; one to the array of pointers, and then from there to an array of doubles.

Although it looks the same, this is less efficient than a multidimensional C array like double[4][4] - being a fixed size makes it possible for the lookup to index directly into the correct member without the additional dereference step; element [2][3] is always the same fixed distance from the beginning, and the whole array can be allocated as a single contiguous block.

The equivalent thing with vectors would not be a vector of vectors (which again requires 2 dereferences) but a single vector of size n² which you index into via a function. That could be put together as a matrix class:

class SquareMatrix
{
    vector<double> _elements;
    size_t _rank;
  public:
    SquareMatrix(size_t rank) : _rank(rank), _elements(rank*rank)
    { }
    double & elem(size_t i, size_t j)
    {
       return *(_elements.begin() + i*_rank + j);
    }
};

Usage:

SquareMatrix a(10);
a.elem(2,3) = 5;
std::cout << a.elem(2,3);

The improvement from this sort of change should be quite significant, partly from the reduction in dereferencing, but most significantly because the CPU will be able to usually cache the right values - much of the time the next value is adjacent to the prior one. This kind of predictable memory usage is also friendly to branch prediction.

Copying copying

The code is creating a copy of each minor in order to calculate its determinant; so for each minor you copy (n-1)² doubles in. That small part itself could be improved by reusing the matrix; instead of copying columns b,c,d in, then copying a,c,d in, you could just replace the relevant values:

a b c d    f g h      a c  d
e f g h -> j k l then i k* l*
i j k l    n o p      m o* p*
m n o p 

See that the marked elements have not changed. So at the very least this part could change to only copy in 2n-1 elements rather than (n-1)². At 50x50 that means 96% fewer copies.

Small things

I would prefer pre-increment (++value) over post-increment (value++) because it is always unambiguous and thus a preferable default even when either would work.

You can do sign = -sign instead of sign = -1 * sign; I don't know whether the compiler will optimise it away, but a negation is much faster than a multiplication.

Memory management in C++ is not like in Fortran

Your code leaks almost all the memory it allocates. For a 50x50 matrix it would leak 50 matrices of size 49x49, and recurse to 50 rank-49 determinants which each leak 49 matrices of size 48x48, which is something like O((n+1)!) leaking.

calloc is a C thing; the C++ equivalent is to call new. Either way, if you calloc some memory you must call free on it to release it, and if you new something you must delete it.

Generally in modern C++ you avoid calling either, and instead rely on structures that manage the memory for you, like vector<double> - this will free the memory it allocates when the vector goes out of scope, which is probably what you're expecting if you're used to Fortran 95/2003.

If you still want to allocate memory via new, then you can use a smart pointer to ensure that the memory is freed when the array goes out of scope. In C++11 that would likely be a unique_ptr or a shared_ptr. Note that you couldn't just wrap your array of pointers in the smart pointer, or it would only free that array and not all the arrays those pointers point to. It would need to be an array of smart pointers, so that each deletes its array when the top level array deconstructor runs.

Summary

Whilst this is a good initial attempt, and not a bad coding exercise, this code will not achieve what you want to achieve. Fixing the leaks will be an educational exercise, but will not result in code that can handle the sizes of matrices that you state you need to cope with. In general anything that needs to cope with matrices of any significant size requires a different kind of method - generally matrix methods that can be done by hand do not scale well.