var PCA = (function () {
/**
* The first step is to subtract the mean and center data
*
* @param {Array} matrix - data in an mXn matrix format
* @returns
*/
function computeDeviationMatrix(matrix) {
var unit = unitSquareMatrix(matrix.length);
return subtract(matrix, scale(multiply(unit, matrix), 1 / matrix.length));
}
/**
* Computes variance from deviation
*
* @param {Array} deviation - data minus mean as calculated from computeDeviationMatrix
* @returns
*/
function computeDeviationScores(deviation) {
var devSumOfSquares = multiply(transpose(deviation), deviation);
return devSumOfSquares;
}
/**
* Calculates the var covar square matrix using either population or sample
*
* @param {Array} devSumOfSquares
* @param {boolean} sample - true/false whether data is from sample or not
* @returns
*/
function computeVarianceCovariance(devSumOfSquares, sample) {
var varianceCovariance;
if (sample)
varianceCovariance = scale(devSumOfSquares, 1 / (devSumOfSquares.length - 1));
else
varianceCovariance = scale(devSumOfSquares, 1 / (devSumOfSquares.length));
return varianceCovariance;
}
/**
* Matrix is the deviation sum of squares as computed earlier
*
* @param {Array} matrix - output of computeDeviationScores
* @returns
*/
function computeSVD(matrix) {
var result = svd(matrix);
var eigenvectors = result.U;
var eigenvalues = result.S;
var results = eigenvalues.map(function (value, i) {
var obj = {};
obj.eigenvalue = value;
obj.vector = eigenvectors.map(function (vector, j) {
return -1 * vector[i]; //HACK prevent completely negative vectors
});
return obj;
});
return results;
}
/**
* Get reduced dataset after removing some dimensions
*
* @param {Array} data - initial matrix started out with
* @param {rest} vectors - eigenvectors selected as part of process
* @returns
*/
function computeAdjustedData(data, ...vectorObjs) {
//FIXME no need to transpose vectors since they're already in row normal form
var vectors = vectorObjs.map(function(v){return v.vector});
var matrixMinusMean = computeDeviationMatrix(data);
var adjustedData = multiply(vectors, transpose(matrixMinusMean));
var unit = unitSquareMatrix(data.length);
var avgData = scale(multiply(unit, data), -1 / data.length); //NOTE get the averages to add back
var formattedAdjustData = formatData(adjustedData, 2);
return {
adjustedData: adjustedData,
formattedAdjustedData: formattedAdjustData,
avgData: avgData,
selectedVectors: vectors
};
}
/**
* Get original data set from reduced data set (decompress)
* @param {*} adjustedData = formatted or unformatted adjusted data
* @param {*} vectors = selectedVectors
* @param {*} avgData = avgData
*/
function computeOriginalData(adjustedData, vectors, avgData) {
var originalWithoutMean = transpose(multiply(transpose(vectors), adjustedData));
var originalWithMean = subtract(originalWithoutMean, avgData);
var formattedData = formatData(originalWithMean, 2);
return {
originalData: originalWithMean,
formattedOriginalData: formattedData
}
}
/**
* Get percentage explained, or loss
* @param {*} vectors
* @param {*} selected
*/
function computePercentageExplained(vectors, ...selected) {
var total = vectors.map(function (v) {
return v.eigenvalue
}).reduce(function (a, b) {
return a + b;
});
var explained = selected.map(function (v) {
return v.eigenvalue
}).reduce(function (a, b) {
return a + b;
});
return (explained / total);
}
function getEigenVectors(data) {
return computeSVD(computeVarianceCovariance(computeDeviationScores(computeDeviationMatrix(data)), false));
}
function analyseTopResult(data) {
var eigenVectors = getEigenVectors(data);
var sorted = eigenVectors.sort(function (a, b) {
return b.eigenvalue - a.eigenvalue;
});
console.log('Sorted Vectors', sorted);
var selected = sorted[0].vector;
return computeAdjustedData(data, selected);
}
function formatData(data, precision) {
var TEN = Math.pow(10, precision || 2);
return data.map(function (d, i) {
return d.map(function (n) {
return Math.round(n * TEN) / TEN;
})
})
}
/**
* Multiplies AxB, where A and B are matrices of nXm and mXn dimensions
* @param {} a
* @param {*} b
*/
function multiply(a, b) {
if (!a[0] || !b[0] || !a.length || !b.length) {
throw new Error('Both A and B should be matrices');
}
if (b.length !== a[0].length) {
throw new Error('Columns in A should be the same as the number of rows in B');
}
var product = [];
for (var i = 0; i < a.length; i++) {
product[i] = []; //initialize a new row
for (var j = 0; j < b[0].length; j++) {
for (var k = 0; k < a[0].length; k++) {
(product[i])[j] = !!(product[i])[j] ? (product[i])[j] + (a[i])[k] * (b[k])[j] : (a[i])[k] * (b[k])[j];
}
}
}
return product;
}
/**
* Utility function to subtract matrix b from a
*
* @param {any} a
* @param {any} b
* @returns
*/
function subtract(a, b) {
if (!(a.length === b.length && a[0].length === b[0].length))
throw new Error('Both A and B should have the same dimensions');
var result = [];
for (var i = 0; i < a.length; i++) {
result[i] = [];
for (var j = 0; j < b[0].length; j++) {
(result[i])[j] = (a[i])[j] - (b[i])[j];
}
}
return result;
}
/**
* Multiplies a matrix into a factor
*
* @param {any} matrix
* @param {any} factor
* @returns
*/
function scale(matrix, factor) {
var result = [];
for (var i = 0; i < matrix.length; i++) {
result[i] = [];
for (var j = 0; j < matrix[0].length; j++) {
(result[i])[j] = (matrix[i])[j] * factor;
}
}
return result;
}
/**
* Generates a unit square matrix
* @param {*} rows = number of rows to fill
*/
function unitSquareMatrix(rows) {
var result = [];
for (var i = 0; i < rows; i++) {
result[i] = [];
for (var j = 0; j < rows; j++) {
(result[i])[j] = 1;
}
}
return result;
}
/**
* Transposes a matrix, converts rows to columns
* @param {*} matrix
*/
function transpose(matrix) {
var operated = clone(matrix);
return operated[0].map(function (m, c) {
return matrix.map(function (r) {
return r[c];
});
});
}
/**
* Deep Clones a matrix
* @param {*} arr
*/
function clone(arr) {
var string = JSON.stringify(arr);
var result = JSON.parse(string);
return result;
}
/**
* Compute the thin SVD from G. H. Golub and C. Reinsch, Numer. Math. 14, 403-420 (1970)
* From the Numeric JS Implementation Copyright (C) 2011 by Sébastien Loisel
* The C implementation from which this has been taken may be found here: http://www.public.iastate.edu/~dicook/JSS/paper/code/svd.c
* @param {*} A = m*n matrix
*/
function svd(A) {
var temp;
var prec = Math.pow(2, -52) // assumes double prec
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 = clone(A);
var m = u.length;
var n = u[0].length;
if (m < n) throw "Need more rows than columns"
var e = new Array(n); //vector1
var q = new Array(n); //vector2
for (i = 0; i < n; i++) e[i] = q[i] = 0.0;
var v = rep([n, n], 0);
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))
}
//rep function
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;
}
//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; //vector
s = 0.0; //sum
l = i + 1; //stays 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 transformations
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 //u is array, v is square of columns
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
}
}
for (i = 0; i < q.length; i++)
if (q[i] < prec) q[i] = 0
//sort eigenvalues
for (i = 0; i < n; i++) {
for (j = i - 1; j >= 0; j--) {
if (q[j] < q[i]) {
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;
}
i = j
}
}
}
return {
U: u,
S: q,
V: v
}
}
return {
computeDeviationScores: computeDeviationScores,
computeDeviationMatrix: computeDeviationMatrix,
computeSVD: computeSVD,
computePercentageExplained: computePercentageExplained,
computeOriginalData: computeOriginalData,
computeVarianceCovariance: computeVarianceCovariance,
computeAdjustedData: computeAdjustedData,
getEigenVectors: getEigenVectors,
analyseTopResult: analyseTopResult,
transpose: transpose,
multiply: multiply,
clone: clone,
scale: scale
}
})();
if(typeof module !== 'undefined')
module.exports = PCA;