/**
* @classdesc A library for inverting matrices and pseudo inverse
* @class matrix/invert
* @hideconstructor
*/
var U = require('./mUtils');
var matrix = require('../matrix');
//Inverts non singular matrix
//Source: http://blog.acipo.com/matrix-inversion-in-javascript/
function matrix_invert(M){
// I use Guassian Elimination to calculate the inverse:
// (1) 'augment' the matrix (left) by the identity (on the right)
// (2) Turn the matrix on the left into the identity by elemetry row ops
// (3) The matrix on the right is the inverse (was the identity matrix)
// There are 3 elemtary row ops: (I combine b and c in my code)
// (a) Swap 2 rows
// (b) Multiply a row by a scalar
// (c) Add 2 rows
//if the matrix isn't square: exit (error)
if(M.length !== M[0].length){
throw 'Matrix must be square to invert';
}
//create the identity matrix (I), and a copy (C) of the original
var i=0, ii=0, j=0, dim=M.length, e=0, t=0;
var I = [], C = [];
for(i=0; i<dim; i+=1){
// Create the row
I[I.length]=[];
C[C.length]=[];
for(j=0; j<dim; j+=1){
//if we're on the diagonal, put a 1 (for identity)
if(i==j){ I[i][j] = 1; }
else{ I[i][j] = 0; }
// Also, make the copy of the original
C[i][j] = M[i][j];
}
}
// Perform elementary row operations
for(i=0; i<dim; i+=1){
// get the element e on the diagonal
e = C[i][i];
// if we have a 0 on the diagonal (we'll need to swap with a lower row)
if(e==0){
//look through every row below the i'th row
for(ii=i+1; ii<dim; ii+=1){
//if the ii'th row has a non-0 in the i'th col
if(C[ii][i] != 0){
//it would make the diagonal have a non-0 so swap it
for(j=0; j<dim; j++){
e = C[i][j]; //temp store i'th row
C[i][j] = C[ii][j];//replace i'th row by ii'th
C[ii][j] = e; //repace ii'th by temp
e = I[i][j]; //temp store i'th row
I[i][j] = I[ii][j];//replace i'th row by ii'th
I[ii][j] = e; //repace ii'th by temp
}
//don't bother checking other rows since we've swapped
break;
}
}
//get the new diagonal
e = C[i][i];
//if it's still 0, not invertable (error)
if(e==0){
console.log('Matrix is singular and cannot be inverted. Attempting SVD');
return pinv(M);
}
}
// Scale this row down by e (so we have a 1 on the diagonal)
for(j=0; j<dim; j++){
C[i][j] = C[i][j]/e; //apply to original matrix
I[i][j] = I[i][j]/e; //apply to identity
}
// Subtract this row (scaled appropriately for each row) from ALL of
// the other rows so that there will be 0's in this column in the
// rows above and below this one
for(ii=0; ii<dim; ii++){
// Only apply to other rows (we want a 1 on the diagonal)
if(ii==i){continue;}
// We want to change this element to 0
e = C[ii][i];
// Subtract (the row above(or below) scaled by e) from (the
// current row) but start at the i'th column and assume all the
// stuff left of diagonal is 0 (which it should be if we made this
// algorithm correctly)
for(j=0; j<dim; j++){
C[ii][j] -= e*C[i][j]; //apply to original matrix
I[ii][j] -= e*I[i][j]; //apply to identity
}
}
}
//we've done all operations, C should be the identity
//matrix I should be the inverse:
return I;
}
/**
* Calculates the inverse of a matrix via Gaussian elimination
* @function invert
* @param {matrix} self - Adds a scalar/array/matrix to self
* @returns {matrix}
* @memberof matrix/invert
* @example
* // Invert a non-singular matrix
* var M = require('./src/matrix_lib');
* var A = new M([[1,2,3],[4,5,6],[7,8,5]]); // non singular matrix
* console.log(A.invert().print());
* // result:
* //-1.917 1.167 -0.250
* //1.833 -1.333 0.500
* //-0.250 0.500 -0.250
*/
matrix.prototype.invert = function()
{
var M = matrix_invert(this.value);
return matrix.make(M);
};
matrix.invert = function(A)
{
return new matrix(matrix_invert(matrix.make(A).value));
}
var localEps = 2.220446049250313e-16;
function pinv(A) {
var z = svd(A), foo = z.S[0];
var U = z.U, S = z.S, V = z.V;
var m = A.length, n = A[0].length, tol = Math.max(m,n)*localEps*foo,M = S.length;
var i,Sinv = new Array(M);
for(i=M-1;i!==-1;i--) { if(S[i]>tol) Sinv[i] = 1/S[i]; else Sinv[i] = 0; }
//return numeric.dot(numeric.dot(V,numeric.diag(Sinv)),numeric.transpose(U))
var SinvM = [];
for(var i = 0; i < Sinv.length; i++)
{
SinvM[i] = [];
for(var j = 0; j < Sinv.length; j++)
{
var sinvVal = (i==j) ? Sinv[i] : 0;
SinvM[i][j] = sinvVal;
}
}
SinvM = matrix.make(SinvM);
V = matrix.make(V);
U = matrix.make(U);
return V.multiply(SinvM).multiply(U.transpose());
}
matrix.prototype.pinv = function()
{
var M = pinv(this.value);
return matrix.make(M);
};
/**
* Calculates the pseudoinverse of a matrix via SVD
* @function pinv
* @param {matrix} self - Adds a scalar/array/matrix to self
* @returns {matrix}
* @memberof matrix/invert
* @example
* // Invert a singular matrix
* var M = require('./src/matrix_lib');
* var A = new M([[1,2,3],[4,5,6],[7,8,9]]); // Singular matrix
* console.log(A.pinv().print());
* // result:
* //-0.639 -0.167 0.306
* //-0.056 -0.000 0.056
* //0.528 0.167 -0.194
*/
matrix.pinv = function(A)
{
return new matrix(pinv(matrix.make(A).value));
}
function rep(s,v,k) {
if(typeof k === "undefined") { k=0; }
var n = s[k], ret = Array(n), i;
if(k === s.length-1) {
for(i=n-2;i>=0;i-=2) { ret[i+1] = v; ret[i] = v; }
if(i===-1) { ret[0] = v; }
return ret;
}
for(i=n-1;i>=0;i--) { ret[i] = rep(s,v,k+1); }
return ret;
}
//Source:Numerics library
//Taken from https://github.com/sloisel/numeric/blob/master/src/svd.js.
//Left as is except for small modificaitons.
//Shanti Rao sent me this routine by private email. I had to modify it
//slightly to work on Arrays instead of using a Matrix object.
//It is apparently translated from http://stitchpanorama.sourceforge.net/Python/svd.py
function svd(A) {
var temp;
//Compute the thin SVD from G. H. Golub and C. Reinsch, Numer. Math. 14, 403-420 (1970)
//var prec= numeric.epsilon; //Math.pow(2,-52) // assumes double prec
var prec = 2.220446049250313e-16;
var tolerance= 1.e-64/prec;
var itmax= 50;
var c=0;
var i=0;
var j=0;
var k=0;
var l=0;
//var u= numeric.clone(A);
var u = U.matrix_copy(A);
var m= u.length;
var n= u[0].length;
if (m < n) throw "Need more rows than columns"
var e = new Array(n);
var q = new Array(n);
for (i=0; i<n; i++) e[i] = q[i] = 0.0;
var v = rep([n,n],0);
// v.zero();
function pythag(a,b)
{
a = Math.abs(a)
b = Math.abs(b)
if (a > b)
return a*Math.sqrt(1.0+(b*b/a/a))
else if (b == 0.0)
return a
return b*Math.sqrt(1.0+(a*a/b/b))
}
//Householder's reduction to bidiagonal form
var f= 0.0;
var g= 0.0;
var h= 0.0;
var x= 0.0;
var y= 0.0;
var z= 0.0;
var s= 0.0;
for (i=0; i < n; i++)
{
e[i]= g;
s= 0.0;
l= i+1;
for (j=i; j < m; j++)
s += (u[j][i]*u[j][i]);
if (s <= tolerance)
g= 0.0;
else
{
f= u[i][i];
g= Math.sqrt(s);
if (f >= 0.0) g= -g;
h= f*g-s
u[i][i]=f-g;
for (j=l; j < n; j++)
{
s= 0.0
for (k=i; k < m; k++)
s += u[k][i]*u[k][j]
f= s/h
for (k=i; k < m; k++)
u[k][j]+=f*u[k][i]
}
}
q[i]= g
s= 0.0
for (j=l; j < n; j++)
s= s + u[i][j]*u[i][j]
if (s <= tolerance)
g= 0.0
else
{
f= u[i][i+1]
g= Math.sqrt(s)
if (f >= 0.0) g= -g
h= f*g - s
u[i][i+1] = f-g;
for (j=l; j < n; j++) e[j]= u[i][j]/h
for (j=l; j < m; j++)
{
s=0.0
for (k=l; k < n; k++)
s += (u[j][k]*u[i][k])
for (k=l; k < n; k++)
u[j][k]+=s*e[k]
}
}
y= Math.abs(q[i])+Math.abs(e[i])
if (y>x)
x=y
}
// accumulation of right hand gtransformations
for (i=n-1; i != -1; i+= -1)
{
if (g != 0.0)
{
h= g*u[i][i+1]
for (j=l; j < n; j++)
v[j][i]=u[i][j]/h
for (j=l; j < n; j++)
{
s=0.0
for (k=l; k < n; k++)
s += u[i][k]*v[k][j]
for (k=l; k < n; k++)
v[k][j]+=(s*v[k][i])
}
}
for (j=l; j < n; j++)
{
v[i][j] = 0;
v[j][i] = 0;
}
v[i][i] = 1;
g= e[i]
l= i
}
// accumulation of left hand transformations
for (i=n-1; i != -1; i+= -1)
{
l= i+1
g= q[i]
for (j=l; j < n; j++)
u[i][j] = 0;
if (g != 0.0)
{
h= u[i][i]*g
for (j=l; j < n; j++)
{
s=0.0
for (k=l; k < m; k++) s += u[k][i]*u[k][j];
f= s/h
for (k=i; k < m; k++) u[k][j]+=f*u[k][i];
}
for (j=i; j < m; j++) u[j][i] = u[j][i]/g;
}
else
for (j=i; j < m; j++) u[j][i] = 0;
u[i][i] += 1;
}
// diagonalization of the bidiagonal form
prec= prec*x
for (k=n-1; k != -1; k+= -1)
{
for (var iteration=0; iteration < itmax; iteration++)
{ // test f splitting
var test_convergence = false
for (l=k; l != -1; l+= -1)
{
if (Math.abs(e[l]) <= prec)
{ test_convergence= true
break
}
if (Math.abs(q[l-1]) <= prec)
break
}
if (!test_convergence)
{ // cancellation of e[l] if l>0
c= 0.0
s= 1.0
var l1= l-1
for (i =l; i<k+1; i++)
{
f= s*e[i]
e[i]= c*e[i]
if (Math.abs(f) <= prec)
break
g= q[i]
h= pythag(f,g)
q[i]= h
c= g/h
s= -f/h
for (j=0; j < m; j++)
{
y= u[j][l1]
z= u[j][i]
u[j][l1] = y*c+(z*s)
u[j][i] = -y*s+(z*c)
}
}
}
// test f convergence
z= q[k]
if (l== k)
{ //convergence
if (z<0.0)
{ //q[k] is made non-negative
q[k]= -z
for (j=0; j < n; j++)
v[j][k] = -v[j][k]
}
break //break out of iteration loop and move on to next k value
}
if (iteration >= itmax-1)
throw 'Error: no convergence.'
// shift from bottom 2x2 minor
x= q[l]
y= q[k-1]
g= e[k-1]
h= e[k]
f= ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y)
g= pythag(f,1.0)
if (f < 0.0)
f= ((x-z)*(x+z)+h*(y/(f-g)-h))/x
else
f= ((x-z)*(x+z)+h*(y/(f+g)-h))/x
// next QR transformation
c= 1.0
s= 1.0
for (i=l+1; i< k+1; i++)
{
g= e[i]
y= q[i]
h= s*g
g= c*g
z= pythag(f,h)
e[i-1]= z
c= f/z
s= h/z
f= x*c+g*s
g= -x*s+g*c
h= y*s
y= y*c
for (j=0; j < n; j++)
{
x= v[j][i-1]
z= v[j][i]
v[j][i-1] = x*c+z*s
v[j][i] = -x*s+z*c
}
z= pythag(f,h)
q[i-1]= z
c= f/z
s= h/z
f= c*g+s*y
x= -s*g+c*y
for (j=0; j < m; j++)
{
y= u[j][i-1]
z= u[j][i]
u[j][i-1] = y*c+z*s
u[j][i] = -y*s+z*c
}
}
e[l]= 0.0
e[k]= f
q[k]= x
}
}
//vt= transpose(v)
//return (u,q,vt)
for (i=0;i<q.length; i++)
if (q[i] < prec) q[i] = 0
//sort eigenvalues
for (i=0; i< n; i++)
{
//writeln(q)
for (j=i-1; j >= 0; j--)
{
if (q[j] < q[i])
{
// writeln(i,'-',j)
c = q[j]
q[j] = q[i]
q[i] = c
for(k=0;k<u.length;k++) { temp = u[k][i]; u[k][i] = u[k][j]; u[k][j] = temp; }
for(k=0;k<v.length;k++) { temp = v[k][i]; v[k][i] = v[k][j]; v[k][j] = temp; }
// u.swapCols(i,j)
// v.swapCols(i,j)
i = j
}
}
}
return {U:u,S:q,V:v}
};