From cbf5df0f7002be2bb39ab7901ff7740c292e285e Mon Sep 17 00:00:00 2001 From: rjbaucells Date: Fri, 24 Apr 2015 00:07:06 -0400 Subject: [PATCH] cs_chol() --- lib/function/algebra/sparse/sparse_chol.js | 157 ++++++++++++++++++ lib/function/algebra/sparse/sparse_cumsum.js | 38 +++++ lib/function/algebra/sparse/sparse_ereach.js | 72 ++++++++ lib/function/algebra/sparse/sparse_symperm.js | 97 +++++++++++ 4 files changed, 364 insertions(+) create mode 100644 lib/function/algebra/sparse/sparse_chol.js create mode 100644 lib/function/algebra/sparse/sparse_cumsum.js create mode 100644 lib/function/algebra/sparse/sparse_ereach.js create mode 100644 lib/function/algebra/sparse/sparse_symperm.js diff --git a/lib/function/algebra/sparse/sparse_chol.js b/lib/function/algebra/sparse/sparse_chol.js new file mode 100644 index 000000000..68d6657f5 --- /dev/null +++ b/lib/function/algebra/sparse/sparse_chol.js @@ -0,0 +1,157 @@ +'use strict'; + +function factory (type, config, load) { + + var divideScalar = load(require('../../arithmetic/divideScalar')); + var sqrt = load(require('../../arithmetic/sqrt')); + var subtract = load(require('../../arithmetic/subtract')); + var multiply = load(require('../../arithmetic/multiply')); + var im = load(require('../../complex/im')); + var re = load(require('../../complex/re')); + var conj = load(require('../../complex/conj')); + var equal = load(require('../../relational/equal')); + var smallerEq = load(require('../../relational/smallerEq')); + + var sparse_symperm = load(require('./sparse_symperm')); + var sparse_ereach = load(require('./sparse_ereach')); + + var CcsMatrix = type.CcsMatrix; + + /** + * Computes the Cholesky factorization of matrix A. It computes L and P so + * L * L' = P * A * P' + * + * @param {Matrix} m The A Matrix to factorize, only upper triangular part used + * @param {Object} s The symbolic analysis from sparse_schol() + * + * @return {Number} The numeric Cholesky factorization of A or null + * + * Reference: http://faculty.cse.tamu.edu/davis/publications.html + */ + var sparse_chol = function (m, s) { + // validate input + if (!m) + return null; + // m arrays + var size = m._size; + // columns + var n = size[1]; + // symbolic analysis result + var parent = s.parent; + var cp = s.cp; + var pinv = s.pinv; + // nonzero elements (estimate) + var lnz = cp[n]; + // L arrays + var lvalues = new Array(lnz); + var lindex = new Array(lnz); + var lptr = new Array(n + 1); + // L + var L = new CcsMatrix({ + values: lvalues, + index: lindex, + ptr: lptr, + size:[n, n] + }); + // vars + var c = new Array(2 * n); + var x = new Array(n); + // compute C = P * A * P' + var cm = pinv ? sparse_symperm (m, pinv, 1) : m; + // C matrix arrays + var cvalues = cm._values; + var cindex = cm._index; + var cptr = cm._ptr; + // vars + var k, p; + // initialize variables + for (k = 0; k < n; k++) + lptr[k] = c[k] = cp[k]; + // compute L(k,:) for L*L' = C + for (k = 0; k < n; k++) { + // nonzero pattern of L(k,:) + var top = sparse_ereach(cm, k, parent, c); + // x (0:k) is now zero + x[k] = 0; + // x = full(triu(C(:,k))) + for (p = cptr[k]; p < cptr[k+1]; p++) { + if (cindex[p] <= k) + x[cindex[p]] = cvalues[p]; + } + // d = C(k,k) + var d = x[k]; + // clear x for k+1st iteration + x[k] = 0; + // solve L(0:k-1,0:k-1) * x = C(:,k) + for (; top < n; top++) { + // s[top..n-1] is pattern of L(k,:) + var i = s[top]; + // L(k,i) = x (i) / L(i,i) + var lki = divideScalar(x[i], lvalues[lptr[i]]); + // clear x for k+1st iteration + x[i] = 0; + for (p = lptr[i] + 1; p < c[i]; p++) { + // row + var r = lindex[p]; + // update x[r] + x[r] = subtract(x[r], multiply(lvalues[p], lki)); + } + // d = d - L(k,i)*L(k,i) + d = subtract(d, multiply(lki, conj(lki))); + p = c[i]++; + // store L(k,i) in column i + lindex[p] = k; + lvalues[p] = conj(lki); + } + // compute L(k,k) + if (smallerEq(re(d), 0) || !equal(im(d), 0)) { + // not pos def + return null; + } + p = c[k]++; + // store L(k,k) = sqrt(d) in column k + lindex[p] = k; + lvalues[p] = sqrt(d); + } + // finalize L + lptr[n] = cp[n]; + // P matrix + var P; + // check we need to calculate P + if (pinv) { + // P arrays + var pvalues = []; + var pindex = []; + var pptr = new Array(n + 1); + // create P matrix + for (p = 0; p < n; p++) { + // initialize ptr (one value per column) + pptr[p] = p; + // index (apply permutation vector) + pindex.push(pinv[p]); + // value 1 + pvalues.push(1); + } + // update ptr + pptr[n] = n; + // P + P = new CcsMatrix({ + values: pvalues, + index: pindex, + ptr: pptr, + size: [n, n] + }); + } + // return L & P + return { + L: L, + P: P + }; + }; + + return sparse_chol; +} + +exports.name = 'sparse_chol'; +exports.path = 'sparse'; +exports.factory = factory; diff --git a/lib/function/algebra/sparse/sparse_cumsum.js b/lib/function/algebra/sparse/sparse_cumsum.js new file mode 100644 index 000000000..1e45627d0 --- /dev/null +++ b/lib/function/algebra/sparse/sparse_cumsum.js @@ -0,0 +1,38 @@ +'use strict'; + +function factory () { + + /** + * It sets the p[i] equal to the sum of c[0] through c[i-1]. + * + * @param {Array} ptr The CCS matrix ptr array + * @param {Array} c The CCS matrix ptr array + * @param {Number} n The number of columns + * + * Reference: http://faculty.cse.tamu.edu/davis/publications.html + */ + var sparse_cumsum = function (ptr, c, n) { + // variables + var i; + var nz = 0; + + for (i = 0; i < n; i++) { + // initialize ptr @ i + ptr[i] = nz; + // increment number of nonzeros + nz += c[i]; + // also copy p[0..n-1] back into c[0..n-1] + c[i] = ptr[i]; + } + // finalize ptr + ptr[n] = nz; + // return sum (c [0..n-1]) + return nz; + }; + + return sparse_cumsum; +} + +exports.name = 'sparse_cumsum'; +exports.path = 'sparse'; +exports.factory = factory; diff --git a/lib/function/algebra/sparse/sparse_ereach.js b/lib/function/algebra/sparse/sparse_ereach.js new file mode 100644 index 000000000..aa78a9321 --- /dev/null +++ b/lib/function/algebra/sparse/sparse_ereach.js @@ -0,0 +1,72 @@ +'use strict'; + +function factory (type, config, load) { + + var sparse_marked = load(require('./sparse_marked')); + var sparse_mark = load(require('./sparse_mark')); + + /** + * Find nonzero pattern of Cholesky L(k,1:k-1) using etree and triu(A(:,k)) + * + * @param {Matrix} a The A matrix + * @param {Number} k The kth column in A + * @param {Array} parent The parent vector from the symbolic analysis result + * @param {Array} w The nonzero pattern xi[top] .. xi[n - 1], an array of size = 2 * n + * The first n entries is the nonzero pattern, the last n entries is the stack + * + * @return {Number} The index for the nonzero pattern + * + * Reference: http://faculty.cse.tamu.edu/davis/publications.html + */ + var sparse_ereach = function (a, k, parent, w) { + // a arrays + var aindex = a._index; + var aptr = a._ptr; + var asize = a._size; + // columns + var n = asize[1]; + // initialize top + var top = n; + // vars + var p, p0, p1, len; + // mark node k as visited + sparse_mark(w, k); + // loop values & index for column k + for (p0 = aptr[k], p1 = aptr[k + 1], p = p0; p < p1; p++) { + // A(i,k) is nonzero + var i = aindex[p]; + // only use upper triangular part of A + if (i > k) + continue; + // traverse up etree + for (len = 0; !sparse_marked(w, i); i = parent[i]) { + // L(k,i) is nonzero, last n entries in w + w[n + len++] = i; + // mark i as visited + sparse_mark(w, i); + } + while (len > 0) { + // decrement top & len + --top; + --len; + // push path onto stack, last n entries in w + w[n + top] = w[n + len]; + } + } + // unmark all nodes + for (p = top; p < n; p++) { + // use stack value, last n entries in w + sparse_mark(w, w[n + p]); + } + // unmark node k + sparse_mark(w, k); + // s[top..n-1] contains pattern of L(k,:) + return top; + }; + + return sparse_ereach; +} + +exports.name = 'sparse_ereach'; +exports.path = 'sparse'; +exports.factory = factory; diff --git a/lib/function/algebra/sparse/sparse_symperm.js b/lib/function/algebra/sparse/sparse_symperm.js new file mode 100644 index 000000000..bd9eb7a0c --- /dev/null +++ b/lib/function/algebra/sparse/sparse_symperm.js @@ -0,0 +1,97 @@ +'use strict'; + +function factory (type, config, load) { + + var sparse_cumsum = load(require('./sparse_cumsum')); + var conj = load(require('../../complex/conj')); + + var CcsMatrix = type.CcsMatrix; + + /** + * Computes the symmetric permutation of matrix A accessing only + * the upper triangular part of A. + * + * C = P * A * P' + * + * @param {Matrix} a The A matrix + * @param {Array} pinv The inverse of permutation vector + * @param {boolean} values Process matrix values (true) + * + * @return {Matrix} The C matrix, C = P * A * P' + * + * Reference: http://faculty.cse.tamu.edu/davis/publications.html + */ + var sparse_symperm = function (a, pinv, values) { + // A matrix arrays + var avalues = a._values; + var aindex = a._index; + var aptr = a._ptr; + var asize = a._size; + // columns + var n = asize[1]; + // number of nonzero elements in C + var nz = aptr[n]; + // C matrix arrays + var cvalues = values && avalues ? new Array(nz) : null; + var cindex = new Array(nz); + var cptr = new Array(n + 1); + // variables + var i, i2, j, j2, p, p0, p1; + // create workspace vector + var w = new Array(n); + // count entries in each column of C + for (j = 0; j < n; j++) { + // column j of A is column j2 of C + j2 = pinv ? pinv[j] : j; + // loop values in column j + for (p0 = aptr[j], p1 = aptr[j + 1], p = p0; p < p1; p++) { + // row + i = aindex[p]; + // skip lower triangular part of A + if (i > j) + continue; + // row i of A is row i2 of C + i2 = pinv ? pinv[i] : i; + // column count of C + w[Math.max(i2, j2)]++; + } + } + // compute column pointers of C + sparse_cumsum(cptr, w, n); + // loop columns + for (j = 0; j < n; j++) { + // column j of A is column j2 of C + j2 = pinv ? pinv[j] : j; + // loop values in column j + for (p0 = aptr[j], p1 = aptr[j + 1], p = p0; p < p1; p++) { + // row + i = aindex[p]; + // skip lower triangular part of A + if (i > j) + continue; + // row i of A is row i2 of C + i2 = pinv ? pinv[i] : i; + // C index for column j2 + var q = w[Math.max(i2, j2)]++; + // update C index for entry q + cindex[q] = Math.min(i2, j2); + // check we need to process values + if (cvalues) + cvalues[q] = (i2 <= j2) ? avalues[p] : conj(avalues[p]); + } + } + // return C matrix + return new CcsMatrix({ + values: cvalues, + index: cindex, + ptr: cptr, + size: [n, n] + }); + }; + + return sparse_symperm; +} + +exports.name = 'sparse_symperm'; +exports.path = 'sparse'; +exports.factory = factory;