Merge pull request #1118 from josdejong/determinant-lu-decomp

Improved performance of determinant

lib/function/matrix/det.js

@@ -10,6 +10,7 @@ function factory (type, config, load, typed) {
   var subtract = load(require('../arithmetic/subtract'));
   var multiply = load(require('../arithmetic/multiply'));
   var unaryMinus = load(require('../arithmetic/unaryMinus'));
+  var lup = load(require('../algebra/decomposition/lup'));
 
   /**
    * Calculate the determinant of a matrix.
@@ -116,49 +117,40 @@ function factory (type, config, load, typed) {
       );
     }
     else {
-      // this is an n x n matrix
-      var compute_mu = function (matrix) {
-        var i, j;
-
-        // Compute the matrix with zero lower triangle, same upper triangle,
-        // and diagonals given by the negated sum of the below diagonal
-        // elements.
-        var mu = new Array(matrix.length);
-        var sum = 0;
-        for (i = 1; i < matrix.length; i++) {
-          sum = add(sum, matrix[i][i]);
-        }
-
-        for (i = 0; i < matrix.length; i++) {
-          mu[i] = new Array(matrix.length);
-          mu[i][i] = unaryMinus(sum);
 
-          for (j = 0; j < i; j++) {
-            mu[i][j] = 0; // TODO: make bignumber 0 in case of bignumber computation
-          }
+      // Compute the LU decomposition
+      var decomp = lup(matrix);
 
-          for (j = i + 1; j < matrix.length; j++) {
-            mu[i][j] = matrix[i][j];
-          }
+      // The determinant is the product of the diagonal entries of U (and those of L, but they are all 1)
+      var det = decomp.U[0][0];
+      for(var i=1; i<rows; i++) {
+        det = multiply(det, decomp.U[i][i]);
+      }
 
-          if (i+1 < matrix.length) {
-            sum = subtract(sum, matrix[i + 1][i + 1]);
-          }
+      // The determinant will be multiplied by 1 or -1 depending on the parity of the permutation matrix.
+      // This can be determined by counting the cycles. This is roughly a linear time algorithm.
+      var evenCycles=0;
+      var i=0;
+      var visited=[];
+      while(true) {
+        while(visited[i]) {
+         i++;
+        }
+        if(i >= rows) break;
+        var j=i;
+        var cycleLen = 0;
+        while(!visited[decomp.p[j]]) {
+          visited[decomp.p[j]] = true;
+          j = decomp.p[j];
+          cycleLen++;
+        }
+        if(cycleLen % 2 === 0) {
+          evenCycles++;
         }
-
-        return mu;
-      };
-
-      var fa = matrix;
-      for (var i = 0; i < rows - 1; i++) {
-        fa = multiply(compute_mu(fa), matrix);
       }
 
-      if (rows % 2 == 0) {
-        return unaryMinus(fa[0][0]);
-      } else {
-        return fa[0][0];
-      }
+      return evenCycles % 2 === 0 ? det : unaryMinus(det);
+
     }
   }
 }

test/function/matrix/det.test.js

@@ -5,6 +5,7 @@ var math = require('../../../index');
 var BigNumber = math.type.BigNumber;
 var Complex = math.type.Complex;
 var DenseMatrix = math.type.DenseMatrix;
+var SparseMatrix = math.type.SparseMatrix;
 var det = math.det;
 var diag = math.diag;
 var eye = math.eye;
@@ -32,9 +33,28 @@ describe('det', function() {
       [1,7,5,9,7],
       [2,7,4,3,7]
     ]), -1176);
+    approx.equal(det([
+      [0,7,0,3,7],
+      [1,7,4,3,7],
+      [0,7,4,3,0],
+      [1,7,5,9,7],
+      [2,7,4,3,7]
+    ]), 1176);
     approx.equal(det(diag([4,-5,6])), -120);
   });
 
+  it('should return the determinant of a sparse matrix', function() {
+
+    approx.equal(det(new SparseMatrix([
+      [1,7,4,3,7],
+      [0,7,0,3,7],
+      [0,7,4,3,0],
+      [1,7,5,9,7],
+      [2,7,4,3,7]
+    ])), -1176);
+
+  });
+
   it('should return 1 for the identity matrix',function() {
     assert.equal(det(eye(7)), 1);
     assert.equal(det(eye(2)), 1);