Determinant of a matrix of doubles

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>

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

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

// 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[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.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 = 1248.14;
matrix_vector_version = 1408.68;
matrix_vector_version = -828.282;
matrix_vector_version =  -53.0927;

// add test data to matrix row 1
matrix_vector_version = 1408.68 ;
matrix_vector_version = 1623.81;
matrix_vector_version = -952.41;
matrix_vector_version =  -67.2946 ;

// add test data to matrix row 2
matrix_vector_version = -828.282;
matrix_vector_version = -952.41;
matrix_vector_version =  559.421;
matrix_vector_version =   38.0384;

// add test data to matrix row 3
matrix_vector_version =  -53.0927;
matrix_vector_version =  -67.2946;
matrix_vector_version =   38.0384;
matrix_vector_version =    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.

• Shouldn't your first calloc() call in matrix_allocate_size() use sizeof (double*) rather than sizeof(double)? They may be the same on your platform, so it may work, but it doesn't properly coney the meaning. Aug 2 '18 at 5:40
• A better choice is sizeof *matrix, which clearly connects the allocation to the variable. Aug 2 '18 at 8:12
• Do you have any choice in how the matrix is represented? The array-of-row-pointers can be an unwieldy way to represent a 2D grid, but if you're stuck with that, there's no point suggesting improvements in that respect. Aug 2 '18 at 8:14
• Yes, the call to calloc() should have used size of (double*) because it is sizing an array of pointers to double not an array of double. I have already made some modifications based on the post of Phil H below using a matrix class instead of the dynamically sized array double**. Thanks for the input so far. Aug 2 '18 at 18:38

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(n2.4), and the best actually implemented algorithms are O(n3) (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 1063 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 - being a fixed size makes it possible for the lookup to index directly into the correct member without the additional dereference step; element  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 n2 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)2 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)2. 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.

• I am working on a computational geometry and linear algebra algorithm that I hope can provide some information about complex data. Part of the solution involves the inverse of a matrix and the only way I know how to do that is with a determinant. I have used libraries like Eigen/Dense, but I find that I don't learn much by including template definitions and header files. I think my code demonstrates I have allot to learn. I have incorporated your suggestion for a matrix class and am working on the copying efficiency suggestion you made. How do I go about posting incremental revised code here? Aug 2 '18 at 18:42
• @LMHmedchem: I strongly advise you to analyse the linear algebra problem more before coding, if you need your solution to work at any reasonable size. As an example, LAPACK is a well known library for solving Ax=b problems (without inverting A), eigenvalue Ax=bx problems (without inverting A). Please read johndcook.com/blog/2010/01/19/dont-invert-that-matrix . If you do need to invert the matrix, use a library to invert it; determinant-based methods are for doing 3x3 matrices by hand, no more. Aug 3 '18 at 9:11