feat: add pseudoinverse function based on SVD

mljs · Feb 28, 2017 · 3279a15 · 3279a15
1 parent 29ad468
commit 3279a15
Show file tree

Hide file tree

Showing 3 changed files with 86 additions and 12 deletions.
diff --git a/src/abstractMatrix.js b/src/abstractMatrix.js
@@ -3,6 +3,7 @@
 module.exports = abstractMatrix;
 
 var LuDecomposition = require('./dc/lu');
+var SvDecomposition = require('./dc/svd');
 var arrayUtils = require('ml-array-utils');
 var util = require('./util');
 var MatrixTransposeView = require('./views/transpose');
@@ -1545,6 +1546,32 @@ function abstractMatrix(superCtor) {
                 throw Error('Determinant can only be calculated for a square matrix.');
             }
         }
+
+        /**
+         * Returns inverse of a matrix if it exists or the pseudoinverse
+         * @param {number} threshold - threshold for taking inverse of singular values (default = 1e-15)
+         * @return {Matrix} the (pseudo)inverted matrix.
+         */
+        pseudoInverse(threshold) {
+            if (threshold === undefined) threshold = Number.EPSILON;
+            var svdSolution = new SvDecomposition(this, {autoTranspose: true});
+
+            var U = svdSolution.leftSingularVectors;
+            var V = svdSolution.rightSingularVectors;
+            var s = svdSolution.diagonal;
+
+            for (var i = 0; i < s.length; i++) {
+                if (Math.abs(s[i]) > threshold) {
+                    s[i] = 1.0 / s[i];
+                } else {
+                    s[i] = 0.0;
+                }
+            }
+
+            // convert list to diagonal
+            s = this.constructor[Symbol.species].diag(s);
+            return V.mmul(s.mmul(U.transposeView()));
+        }
     }
 
     Matrix.prototype.klass = 'Matrix';

diff --git a/src/dc/svd.js b/src/dc/svd.js
@@ -1,6 +1,6 @@
 'use strict';
 
-var Matrix = require('../matrix').Matrix;
+var Matrix = require('../matrix');
 var util = require('./util');
 var hypotenuse = util.hypotenuse;
 var getFilled2DArray = util.getFilled2DArray;
@@ -10,7 +10,7 @@ function SingularValueDecomposition(value, options) {
     if (!(this instanceof SingularValueDecomposition)) {
         return new SingularValueDecomposition(value, options);
     }
-    value = Matrix.checkMatrix(value);
+    value = Matrix.Matrix.checkMatrix(value);
 
     options = options || {};
 
@@ -418,26 +418,26 @@ SingularValueDecomposition.prototype = {
         return (Math.pow(2, -52) / 2) * Math.max(this.m, this.n) * this.s[0];
     },
     get leftSingularVectors() {
-        if (!Matrix.isMatrix(this.U)) {
-            this.U = new Matrix(this.U);
+        if (!Matrix.Matrix.isMatrix(this.U)) {
+            this.U = new Matrix.Matrix(this.U);
         }
         return this.U;
     },
     get rightSingularVectors() {
-        if (!Matrix.isMatrix(this.V)) {
-            this.V = new Matrix(this.V);
+        if (!Matrix.Matrix.isMatrix(this.V)) {
+            this.V = new Matrix.Matrix(this.V);
         }
         return this.V;
     },
     get diagonalMatrix() {
-        return Matrix.diag(this.s);
+        return Matrix.Matrix.diag(this.s);
     },
     solve: function (value) {
 
         var Y = value,
             e = this.threshold,
             scols = this.s.length,
-            Ls = Matrix.zeros(scols, scols),
+            Ls = Matrix.Matrix.zeros(scols, scols),
             i;
 
         for (i = 0; i < scols; i++) {
@@ -454,7 +454,7 @@ SingularValueDecomposition.prototype = {
         var VL = V.mmul(Ls),
             vrows = V.rows,
             urows = U.length,
-            VLU = Matrix.zeros(vrows, urows),
+            VLU = Matrix.Matrix.zeros(vrows, urows),
             j, k, sum;
 
         for (i = 0; i < vrows; i++) {
@@ -470,14 +470,14 @@ SingularValueDecomposition.prototype = {
         return VLU.mmul(Y);
     },
     solveForDiagonal: function (value) {
-        return this.solve(Matrix.diag(value));
+        return this.solve(Matrix.Matrix.diag(value));
     },
     inverse: function () {
         var V = this.V;
         var e = this.threshold,
             vrows = V.length,
             vcols = V[0].length,
-            X = new Matrix(vrows, this.s.length),
+            X = new Matrix.Matrix(vrows, this.s.length),
             i, j;
 
         for (i = 0; i < vrows; i++) {
@@ -494,7 +494,7 @@ SingularValueDecomposition.prototype = {
 
         var urows = U.length,
             ucols = U[0].length,
-            Y = new Matrix(vrows, urows),
+            Y = new Matrix.Matrix(vrows, urows),
             k, sum;
 
         for (i = 0; i < vrows; i++) {

diff --git a/test/matrix/utility.js b/test/matrix/utility.js
@@ -176,4 +176,51 @@ describe('utility methods', function () {
         matrix2 = new Matrix([[2, 1, 3], [7, 1, 1], [6, 2, 7]]);
         matrix.strassen3x3(matrix2).to2DArray().should.eql([[38, 8, 17], [33, 12, 36], [35, 12, 37]]);
     });
+
+    it('pseudoinverse', function () {
+        // Actual values calculated by the Numpy library
+
+        var matrix = new Matrix([[2, 4], [7, 1]]);
+        var result = matrix.pseudoInverse().to2DArray();
+
+        result[0][0].should.approximately(-0.03846153846153843, 1e-5);
+        result[0][1].should.approximately(0.15384615384615374, 1e-5);
+        result[1][0].should.approximately(0.2692307692307691, 1e-5);
+        result[1][1].should.approximately(-0.076923076923077, 1e-5);
+
+        matrix = new Matrix([[4, 7], [2, 6]]);
+        result = matrix.pseudoInverse().to2DArray();
+        result[0][0].should.approximately(0.6, 1e-5);
+        result[0][1].should.approximately(-0.7, 1e-5);
+        result[1][0].should.approximately(-0.2, 1e-5);
+        result[1][1].should.approximately(0.4, 1e-5);
+
+        // Example from http://reference.wolfram.com/language/ref/PseudoInverse.html
+        matrix = new Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]);
+        result = matrix.pseudoInverse().to2DArray();
+
+        result[0][0].should.approximately(-6.38888889e-01, 1e-5);
+        result[0][1].should.approximately(-1.66666667e-01, 1e-5);
+        result[0][2].should.approximately(3.05555556e-01, 1e-5);
+
+        result[1][0].should.approximately(-5.55555556e-02, 1e-5);
+        result[1][1].should.approximately(1.34337961e-16, 1e-5);
+        result[1][2].should.approximately(5.55555556e-02, 1e-5);
+
+        result[2][0].should.approximately(5.27777778e-01, 1e-5);
+        result[2][1].should.approximately(1.66666667e-01, 1e-5);
+        result[2][2].should.approximately(-1.94444444e-01, 1e-5);
+
+        // non-square matrix.
+        matrix = new Matrix([[1, 0, 1], [2, 4, 5]]);
+        result = matrix.pseudoInverse().to2DArray();
+
+        result[0][0].should.approximately(0.75609756, 1e-5);
+        result[0][1].should.approximately(-0.07317073, 1e-5);
+        result[1][0].should.approximately(-0.68292683, 1e-5);
+        result[1][1].should.approximately(0.19512195, 1e-5);
+        result[2][0].should.approximately(0.24390244, 1e-5);
+        result[2][1].should.approximately(0.07317073, 1e-5);
+
+    });
 });