Calculating the determinant of a matrix

Question

I wanted to do some exercise and came up with the idea of a good challenge (for my level of course). I tried to implement Laplace's algorithm for computing the determinant, recursively.

#include <math.h>
#include <stdlib.h>
#include <vector>
#include <iostream>
using namespace std;

int getMinimoCount = 0; //ignore these, just to keep track of the recursion.
int calcDetCount = 0;

void printMatrix ( vector< vector<double> > M) {
  //just does what it means
  int size = M.size();
  for( int i = 0; i < size; i++ ) {
    cout << "\t";
    for( int j = 0; j < size; j++ ) {
      cout << M[i][j] << "\t"; }
    cout << endl << endl << endl; }
  cout << endl;
}

vector< vector<double> > getMinimo( vector< vector<double> > src, int I, int J, int ordSrc ) {
  // Compute and return the minimum of the element I J
  // If the element is not in the Ith row or Jth column it will get copied to the minimum matrix
  getMinimoCount++;

  vector< vector<double> > minimo( ordSrc-1, vector<double> (ordSrc-1,0));

  int rowCont = 0;
  for( int i=0; i < ordSrc; i++)
  {
    int colCont = 0;
    if ( i != I ) { 

      for ( int j=0; j < ordSrc; j++)
      { 
        if ( j != J ) { 
          minimo[rowCont][colCont] = src[i][j];
          colCont++; }
      };

      rowCont++; }
  };
  return minimo;
}


double calcDet( vector< vector<double> > src, int ord) {
  // Here be recursion. 
  calcDetCount++;

  if ( ord == 2 ) {

    double mainDiag = src[0][0] * src[1][1];
    double negDiag = src[1][0] * src[0][1];

    return mainDiag - negDiag; }
  else {
    double det = 0;

    for( int J = 0; J < ord; J++) 
    {
      vector< vector<double> > min = getMinimo( src, 0, J, ord);

      if ( (J % 2) == 0 ) { det += src[0][J] * calcDet( min, ord-1); }
      else { det -= src[0][J] * calcDet( min, ord-1); }

    };

    return det;
  }
}

int main() {

  // Just some UI to gather the matrix. not really convinced of this.
  int ord;
  cout << "############## MATRIX DET ##############" << endl << endl;
  cout << " Matrix order: "; cin >> ord; cout << endl;
  vector <vector<double> > mainMatrix( ord, vector<double> (ord, 0));

  cout << """ insert values one row at time. Top to bottom:\n\n""";
  for ( int countY = 0; countY < ord; countY++) {
    for ( int countX = 0; countX < ord; countX++) {
      cin >> mainMatrix[countY][countX];};
  };

  system("CLS");
  cout << "############## MATRIX DET ##############" << endl << endl;
  cout << endl << endl << " This is the input matrix:" << endl << endl << endl;
  printMatrix( mainMatrix );

  system("PAUSE");

  system("CLS");
  cout << "############## MATRIX DET ##############" << endl << endl;
  cout << " Working...!" << endl;
  double det = calcDet( mainMatrix, ord );
  system("CLS");

  cout << endl << endl << "############## MATRIX DET ##############" << endl << endl;
  cout << " Det =\t" << det << endl << endl;
  cout << " getMinimo() chiamata: " << getMinimoCount << " volte" << endl;
  cout << " calcDet() chiamata: " << calcDetCount << " volte" << endl << endl;

  return 0;
}

The concept is simple: you have a matrix of order n. While doing this by hand you'd prefer chosing a row that's particularly math friendly; since it's a computer doing the dirty work it really doesn't care about what number he's multiplying.

Every element a_IJ of a matrix has a minor. A minor is the determinant of the matrix without the I-th row and the J-th column. With this we can define the det of a matrix like so:

Sum (-1)^i+j * a_ij * M_ij

(where M_ij is the minimum of the element a_ij)

Once a matrix reach the order == 2 it just computes the determinant since is just a simple multiplication between 4 elements.

At first I had problems finding something that could be used as a matrix object and be passed around from one function to another. I tried bi-dimensional arrays, but they aren't dynamic and couldn't understand how I could pass an array object to a function. I was looking for some hidden matrix type or class to use but i had no luck what so ever.

I came up with "a vector of vectors" which worked but I'm not completely sure it's a really good idea. Not even counting that vector <vector<double> > looks awful.

Secondly, running some recursion tracking I found out that the time it takes ramps up so damn quickly:

• Order 3 = 4 calls
• Order 4 = 17 calls
• Order 5 = 86 calls
• ...
• Order 10 = 2.606.501 calls

Although, even in the latter case, it doesn't take too much: 3 to 4 seconds. Is there a way to reduce the steepness of the curve? This way it gets out of hand way too soon.

I don't have a huge programming experience so I know almost nothing on optimization or good practice either. Since I have some fresh code to work with I'd like to know any error I might be doing and way to optimize this algorithm.

How costly is to cast a type? The (type) x type of cast to be clear.

What could be a better way to input a matrix from the user?

Answer 1

Includes

#include <math.h>
#include <stdlib.h>

You're using C++. You really ought to #include the C++ versions of the headers instead of the C ones:

#include <cmath>
#include <cstdlib>

Namespace

using namespace std;

This is a very bad idea. It pollutes the global namespace with everything from std. If std gets updated and includes a symbol that conflicts with something in your project, you're in big trouble. It's not much work to prepend std:: to your objects, but if you're really that lazy, just import what you need:

using std::vector;
using std::cin;
using std::cout;
using std::endl;

System calls

system("CLS");

Please don't do this. It makes your code completely non-portable. If you really must use an OS-dependent feature like this, at the very least isolate it in its own function, so you only have to change it one place instead of all through your code.

Algorithm

The way that's taught in high school for calculating the determinant of the matrix is rather inefficient (though simple to apply). The time is proportional to \$n!\$ -- that's right, factorial.

Especially on a computer, it's much better to get the determinant by one of the decomposition methods (basically, extended Gaussian elimination). This is an \$O(n^3)\$ process instead. Check out Wikipedia for some ideas.