mirror of
https://github.com/josdejong/mathjs.git
synced 2026-01-18 14:59:29 +00:00
cs_chol()
This commit is contained in:
rjbaucells
parent
4197743500
commit
cbf5df0f70
157
lib/function/algebra/sparse/sparse_chol.js
Normal file
157
lib/function/algebra/sparse/sparse_chol.js
Normal file
@ -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;
|
||||
38
lib/function/algebra/sparse/sparse_cumsum.js
Normal file
38
lib/function/algebra/sparse/sparse_cumsum.js
Normal file
@ -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;
|
||||
72
lib/function/algebra/sparse/sparse_ereach.js
Normal file
72
lib/function/algebra/sparse/sparse_ereach.js
Normal file
@ -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;
|
||||
97
lib/function/algebra/sparse/sparse_symperm.js
Normal file
97
lib/function/algebra/sparse/sparse_symperm.js
Normal file
@ -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;
|
||||
Loading…
x
Reference in New Issue
Block a user