Maybe this will help.
Check out Array Math.svd(Matrix A) which computes the singular value decomposition of A (square or rectangular). Here is a recap from PTVF Numerical Recipes book:
[u, w, v] = Math.svd(a). u, w, and v form a decomposition of a. u and v are orthogonal, w is a vector that represents a diagonal matrix, a = u.diag(w).transpose(v). Since u and v are orthogonal, inverse(a) is easy to compute inverse(a) = v.diag(1/w).transpose(u). If a is singular or near singular some wi values may be zero or near zero, so 1/wi may be undefined. Solution, amazingly, is to simply replace 1/wi by zero when wi is zero or near zero. This replacement is controlled by a tolerance (i.e. like 1e-10).
Here is code that returns the so-called pseudo inverse. Tested to solve over-determined systems (least squares solutions to over-determined linear systems) on ~30 x 100k matricies with good performance.
Mike
// Gives the pseudo inverse, singular values smaller than threshold will be dropped.
function pseudoInverse(matrix, threshold) {
var svd = Math.svd(matrix);
var u = svd[0];
var w = svd[1];
var v = svd[2];
var iw = new Matrix(w.length, w.length);
for (var row = 0; row != w.length; ++row) {
for (var col = 0; col != w.length; ++col) {
iw.at(row, col, row == col ? Math.abs(w.at(row)) < threshold ? 0 : 1 / w.at(row) : 0);
}
}
return v.mul(iw).mul(u.transpose());
};
