From 0256d8f9672b77cf8fb64c9f5719a78e0db1b030 Mon Sep 17 00:00:00 2001 From: Michal Grno Date: Mon, 13 Jul 2020 12:20:11 +0200 Subject: [PATCH 1/9] refactor solveValidation --- .../algebra/solver/utils/solveValidation.js | 164 +++++++++--------- 1 file changed, 78 insertions(+), 86 deletions(-) diff --git a/src/function/algebra/solver/utils/solveValidation.js b/src/function/algebra/solver/utils/solveValidation.js index 056a8f9db0..eb6dd66965 100644 --- a/src/function/algebra/solver/utils/solveValidation.js +++ b/src/function/algebra/solver/utils/solveValidation.js @@ -1,4 +1,4 @@ -import { isArray, isDenseMatrix, isMatrix } from '../../../../utils/is' +import { isArray, isMatrix, isDenseMatrix, isSparseMatrix } from '../../../../utils/is' import { arraySize } from '../../../../utils/array' import { format } from '../../../../utils/string' @@ -13,131 +13,123 @@ export function createSolveValidation ({ DenseMatrix }) { * @return {DenseMatrix} Dense column vector b */ return function solveValidation (m, b, copy) { - // matrix size - const size = m.size() - // validate matrix dimensions - if (size.length !== 2) { throw new RangeError('Matrix must be two dimensional (size: ' + format(size) + ')') } - // rows & columns - const rows = size[0] - const columns = size[1] - // validate rows & columns - if (rows !== columns) { throw new RangeError('Matrix must be square (size: ' + format(size) + ')') } - // vars - let data, i, bdata - // check b is matrix + const mSize = m.size() + + if (mSize.length !== 2) { + throw new RangeError('Matrix must be two dimensional (size: ' + format(mSize) + ')') + } + + const rows = mSize[0] + const columns = mSize[1] + + if (rows !== columns) { + throw new RangeError('Matrix must be square (size: ' + format(mSize) + ')') + } + + let data = [] + if (isMatrix(b)) { - // matrix size - const msize = b.size() - // vector - if (msize.length === 1) { - // check vector length - if (msize[0] !== rows) { throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') } - // create data array - data = [] - // matrix data (DenseMatrix) - bdata = b._data - // loop b data - for (i = 0; i < rows; i++) { - // row array + const bSize = b.size() + const bdata = b._data + + // 1-dim vector + if (bSize.length === 1) { + if (bSize[0] !== rows) { + throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') + } + + for (let i = 0; i < rows; i++) { data[i] = [bdata[i]] } - // return Dense Matrix + return new DenseMatrix({ data: data, size: [rows, 1], datatype: b._datatype }) } - // two dimensions - if (msize.length === 2) { - // array must be a column vector - if (msize[0] !== rows || msize[1] !== 1) { throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') } - // check matrix type + + // 2-dim column + if (bSize.length === 2) { + if (bSize[0] !== rows || bSize[1] !== 1) { + throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') + } + if (isDenseMatrix(b)) { - // check a copy is needed if (copy) { - // create data array data = [] - // matrix data (DenseMatrix) - bdata = b._data - // loop b data - for (i = 0; i < rows; i++) { - // row array + + for (let i = 0; i < rows; i++) { data[i] = [bdata[i][0]] } - // return Dense Matrix + return new DenseMatrix({ data: data, size: [rows, 1], datatype: b._datatype }) } - // b is already a column vector + return b } - // create data array - data = [] - for (i = 0; i < rows; i++) { data[i] = [0] } - // sparse matrix arrays - const values = b._values - const index = b._index - const ptr = b._ptr - // loop values in column 0 - for (let k1 = ptr[1], k = ptr[0]; k < k1; k++) { - // row - i = index[k] - // add to data - data[i][0] = values[k] + + if (isSparseMatrix(b)) { + for (let i = 0; i < rows; i++) { data[i] = [0] } + + const values = b._values + const index = b._index + const ptr = b._ptr + + for (let k1 = ptr[1], k = ptr[0]; k < k1; k++) { + const i = index[k] + data[i][0] = values[k] + } + + return new DenseMatrix({ + data: data, + size: [rows, 1], + datatype: b._datatype + }) } - // return Dense Matrix - return new DenseMatrix({ - data: data, - size: [rows, 1], - datatype: b._datatype - }) } - // throw error - throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') + + throw new RangeError('Dimension mismatch. The right side has to be either 1- or 2-dimensional vector.') } - // check b is array + if (isArray(b)) { - // size - const asize = arraySize(b) - // check matrix dimensions, vector - if (asize.length === 1) { - // check vector length - if (asize[0] !== rows) { throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') } - // create data array - data = [] - // loop b - for (i = 0; i < rows; i++) { - // row array + const bsize = arraySize(b) + + if (bsize.length === 1) { + if (bsize[0] !== rows) { + throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') + } + + for (let i = 0; i < rows; i++) { data[i] = [b[i]] } - // return Dense Matrix + return new DenseMatrix({ data: data, size: [rows, 1] }) } - if (asize.length === 2) { - // array must be a column vector - if (asize[0] !== rows || asize[1] !== 1) { throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') } - // create data array - data = [] - // loop b data - for (i = 0; i < rows; i++) { - // row array + + if (bsize.length === 2) { + if (bsize[0] !== rows || bsize[1] !== 1) { + throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') + } + + for (let i = 0; i < rows; i++) { data[i] = [b[i][0]] } - // return Dense Matrix + return new DenseMatrix({ data: data, size: [rows, 1] }) } - // throw error - throw new RangeError('Dimension mismatch. Matrix columns must match vector length.') + + throw new RangeError('Dimension mismatch. The right side has to be either 1- or 2-dimensional vector.') } } } From 6d2fef96dd71df059410df43f260c76f959daa74 Mon Sep 17 00:00:00 2001 From: Michal Grno Date: Mon, 13 Jul 2020 13:06:02 +0200 Subject: [PATCH 2/9] refactor usolve --- src/function/algebra/solver/usolve.js | 82 ++++++++++++++------------- 1 file changed, 42 insertions(+), 40 deletions(-) diff --git a/src/function/algebra/solver/usolve.js b/src/function/algebra/solver/usolve.js index 8767ba1753..fe0c4affaf 100644 --- a/src/function/algebra/solver/usolve.js +++ b/src/function/algebra/solver/usolve.js @@ -42,67 +42,64 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed return typed(name, { 'SparseMatrix, Array | Matrix': function (m, b) { - // process matrix return _sparseBackwardSubstitution(m, b) }, 'DenseMatrix, Array | Matrix': function (m, b) { - // process matrix return _denseBackwardSubstitution(m, b) }, 'Array, Array | Matrix': function (a, b) { - // create dense matrix from array const m = matrix(a) - // use matrix implementation const r = _denseBackwardSubstitution(m, b) - // result return r.valueOf() } }) function _denseBackwardSubstitution (m, b) { - // validate matrix and vector, return copy of column vector b + // make b into a column vector b = solveValidation(m, b, true) - // column vector data + const bdata = b._data - // rows & columns + const rows = m._size[0] const columns = m._size[1] + // result const x = [] - // arrays - const data = m._data - // backward solve m * x = b, loop columns (backwards) + + const mdata = m._data + // loop columns backwards for (let j = columns - 1; j >= 0; j--) { // b[j] const bj = bdata[j][0] || 0 // x[j] let xj - // backward substitution (outer product) avoids inner looping when bj === 0 + if (!equalScalar(bj, 0)) { - // value @ [j, j] - const vjj = data[j][j] - // check vjj + // value at [j, j] + const vjj = mdata[j][j] + if (equalScalar(vjj, 0)) { // system cannot be solved throw new Error('Linear system cannot be solved since matrix is singular') } - // calculate xj + xj = divideScalar(bj, vjj) + // loop rows for (let i = j - 1; i >= 0; i--) { // update copy of b - bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, data[i][j]))] + bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, mdata[i][j]))] } } else { - // zero value @ j + // zero value at j xj = 0 } // update x x[j] = [xj] } - // return column vector + return new DenseMatrix({ data: x, size: [rows, 1] @@ -110,37 +107,42 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed } function _sparseBackwardSubstitution (m, b) { - // validate matrix and vector, return copy of column vector b + // make b into a column vector b = solveValidation(m, b, true) - // column vector data + const bdata = b._data - // rows & columns + const rows = m._size[0] const columns = m._size[1] - // matrix arrays + const values = m._values const index = m._index const ptr = m._ptr - // vars + let i, k + // result const x = [] - // backward solve m * x = b, loop columns (backwards) + + // loop columns backwards for (let j = columns - 1; j >= 0; j--) { // b[j] const bj = bdata[j][0] || 0 - // backward substitution (outer product) avoids inner looping when bj === 0 + if (!equalScalar(bj, 0)) { - // value @ [j, j] + // value at [j, j] let vjj = 0 + // upper triangular matrix values & index (column j) - const jvalues = [] - const jindex = [] + const jValues = [] + const jIndices = [] + // first & last indeces in column - const f = ptr[j] - let l = ptr[j + 1] - // values in column, find value @ [j, j], loop backwards - for (k = l - 1; k >= f; k--) { + const firstIndex = ptr[j] + let lastIndex = ptr[j + 1] + + // values in column, find value at [j, j], loop backwards + for (k = lastIndex - 1; k >= firstIndex; k--) { // row i = index[k] // check row @@ -149,23 +151,23 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed vjj = values[k] } else if (i < j) { // store upper triangular - jvalues.push(values[k]) - jindex.push(i) + jValues.push(values[k]) + jIndices.push(i) } } - // at this point we must have a value @ [j, j] + // at this point we must have a value at [j, j] if (equalScalar(vjj, 0)) { - // system cannot be solved, there is no value @ [j, j] + // system cannot be solved, there is no value at [j, j] throw new Error('Linear system cannot be solved since matrix is singular') } // calculate xj const xj = divideScalar(bj, vjj) // loop upper triangular - for (k = 0, l = jindex.length; k < l; k++) { + for (k = 0, lastIndex = jIndices.length; k < lastIndex; k++) { // row - i = jindex[k] + i = jIndices[k] // update copy of b - bdata[i] = [subtract(bdata[i][0], multiplyScalar(xj, jvalues[k]))] + bdata[i] = [subtract(bdata[i][0], multiplyScalar(xj, jValues[k]))] } // update x x[j] = [xj] From 3eb7ad0058e4cb1d961fa08cd99d89875fc6d4f4 Mon Sep 17 00:00:00 2001 From: Michal Grno Date: Tue, 14 Jul 2020 20:30:18 +0200 Subject: [PATCH 3/9] usolve algorithm implemented (for square mat.) --- src/function/algebra/solver/usolve.js | 79 +++++++++++++++------------ 1 file changed, 45 insertions(+), 34 deletions(-) diff --git a/src/function/algebra/solver/usolve.js b/src/function/algebra/solver/usolve.js index fe0c4affaf..f07529c13e 100644 --- a/src/function/algebra/solver/usolve.js +++ b/src/function/algebra/solver/usolve.js @@ -56,54 +56,65 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed } }) - function _denseBackwardSubstitution (m, b) { - // make b into a column vector - b = solveValidation(m, b, true) + function _denseBackwardSubstitution (m, b_) { + // the algorithm is derived from + // https://www.overleaf.com/project/5e6c87c554a3190001a3fc93 - const bdata = b._data + // array of right-hand sides + const B = [solveValidation(m, b_, true)._data.map(e => e[0])] + const M = m._data const rows = m._size[0] const columns = m._size[1] - // result - const x = [] - - const mdata = m._data // loop columns backwards - for (let j = columns - 1; j >= 0; j--) { - // b[j] - const bj = bdata[j][0] || 0 - // x[j] - let xj + for (let i = columns - 1; i >= 0; i--) { + let L = B.length - if (!equalScalar(bj, 0)) { - // value at [j, j] - const vjj = mdata[j][j] + // loop right-hand sides + for (let k = 0; k < L; k++) { + const b = B[k] - if (equalScalar(vjj, 0)) { - // system cannot be solved - throw new Error('Linear system cannot be solved since matrix is singular') - } + if (!equalScalar(M[i][i], 0)) { + // non-singular row - xj = divideScalar(bj, vjj) + b[i] = divideScalar(b[i], M[i][i]) - // loop rows - for (let i = j - 1; i >= 0; i--) { - // update copy of b - bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, mdata[i][j]))] + for (let j = i - 1; j >= 0; j--) { + // b[j] -= b[i] * M[j,i] + b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i])) + } + + } else if (!equalScalar(b[i], 0)) { + // singular row, nonzero RHS + + if (k === 0) { + // There is no valid solution + throw new Error('Linear system cannot be solved since matrix is singular') + } else { + // This RHS is invalid but other solutions may still exist + B.splice(k, 1) + k -= 1 + L -= 1 + } + } else if (k === 0) { + // singular row, RHS is zero + + const bNew = [...b] + bNew[i] = 1 + + for (let j = i - 1; j >= 0; j--) { + bNew[j] = subtract(bNew[j], M[j][i]) + } + + B.push(bNew) } - } else { - // zero value at j - xj = 0 + } - // update x - x[j] = [xj] + } - return new DenseMatrix({ - data: x, - size: [rows, 1] - }) + return B } function _sparseBackwardSubstitution (m, b) { From cebad52c919a4e25f6d0eb3755ad02f67013e886 Mon Sep 17 00:00:00 2001 From: Michal Grno Date: Tue, 14 Jul 2020 22:26:38 +0200 Subject: [PATCH 4/9] added lsolve, consistent return type, fixed tests --- src/function/algebra/solver/lsolve.js | 110 ++++++++++-------- src/function/algebra/solver/usolve.js | 12 +- .../function/algebra/solver/lsolve.test.js | 20 ++-- .../function/algebra/solver/usolve.test.js | 20 ++-- 4 files changed, 85 insertions(+), 77 deletions(-) diff --git a/src/function/algebra/solver/lsolve.js b/src/function/algebra/solver/lsolve.js index d1e5b024dc..c1998f5e63 100644 --- a/src/function/algebra/solver/lsolve.js +++ b/src/function/algebra/solver/lsolve.js @@ -37,76 +37,84 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed * @param {Matrix, Array} L A N x N matrix or array (L) * @param {Matrix, Array} b A column vector with the b values * - * @return {DenseMatrix | Array} A column vector with the linear system solution (x) + * @return {DenseMatrix[] | Array[]} An array of affine-independent column vectors (x) that solve the linear system */ return typed(name, { 'SparseMatrix, Array | Matrix': function (m, b) { - // process matrix return _sparseForwardSubstitution(m, b) }, 'DenseMatrix, Array | Matrix': function (m, b) { - // process matrix return _denseForwardSubstitution(m, b) }, 'Array, Array | Matrix': function (a, b) { - // create dense matrix from array const m = matrix(a) - // use matrix implementation - const r = _denseForwardSubstitution(m, b) - // result - return r.valueOf() + const R = _denseForwardSubstitution(m, b) + return R.map(r => r.valueOf()) } }) - function _denseForwardSubstitution (m, b) { - // validate matrix and vector, return copy of column vector b - b = solveValidation(m, b, true) - // column vector data - const bdata = b._data - // rows & columns + function _denseForwardSubstitution (m, b_) { + // the algorithm is derived from + // https://www.overleaf.com/project/5e6c87c554a3190001a3fc93 + + // array of right-hand sides + const B = [solveValidation(m, b_, true)._data.map(e => e[0])] + + const M = m._data const rows = m._size[0] const columns = m._size[1] - // result - const x = [] - // data - const data = m._data - // forward solve m * x = b, loop columns - for (let j = 0; j < columns; j++) { - // b[j] - const bj = bdata[j][0] || 0 - // x[j] - let xj - // forward substitution (outer product) avoids inner looping when bj === 0 - if (!equalScalar(bj, 0)) { - // value @ [j, j] - const vjj = data[j][j] - // check vjj - if (equalScalar(vjj, 0)) { - // system cannot be solved - throw new Error('Linear system cannot be solved since matrix is singular') - } - // calculate xj - xj = divideScalar(bj, vjj) - // loop rows - for (let i = j + 1; i < rows; i++) { - // update copy of b - bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, data[i][j]))] + + // loop columns + for (let i = 0; i < columns; i++) { + let L = B.length + + // loop right-hand sides + for (let k = 0; k < L; k++) { + const b = B[k] + + if (!equalScalar(M[i][i], 0)) { + // non-singular row + + b[i] = divideScalar(b[i], M[i][i]) + + for (let j = i + 1; j < columns; j++) { + // b[j] -= b[i] * M[j,i] + b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i])) + } + + } else if (!equalScalar(b[i], 0)) { + // singular row, nonzero RHS + + if (k === 0) { + // There is no valid solution + throw new Error('Linear system cannot be solved since matrix is singular') + } else { + // This RHS is invalid but other solutions may still exist + B.splice(k, 1) + k -= 1 + L -= 1 + } + } else if (k === 0) { + // singular row, RHS is zero + + const bNew = [...b] + bNew[i] = 1 + + for (let j = i + 1; j < columns; j++) { + bNew[j] = subtract(bNew[j], M[j][i]) + } + + B.push(bNew) } - } else { - // zero @ j - xj = 0 + } - // update x - x[j] = [xj] + } - // return vector - return new DenseMatrix({ - data: x, - size: [rows, 1] - }) + + return B.map(x => new DenseMatrix({ data: x.map(e=>[e]), size: [rows, 1] })) } function _sparseForwardSubstitution (m, b) { @@ -174,9 +182,9 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed } } // return vector - return new DenseMatrix({ + return [new DenseMatrix({ data: x, size: [rows, 1] - }) + })] } }) diff --git a/src/function/algebra/solver/usolve.js b/src/function/algebra/solver/usolve.js index f07529c13e..5cd5b292f5 100644 --- a/src/function/algebra/solver/usolve.js +++ b/src/function/algebra/solver/usolve.js @@ -37,7 +37,7 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed * @param {Matrix, Array} U A N x N matrix or array (U) * @param {Matrix, Array} b A column vector with the b values * - * @return {DenseMatrix | Array} A column vector with the linear system solution (x) + * @return {DenseMatrix[] | Array[]} An array of affine-independent column vectors (x) that solve the linear system */ return typed(name, { @@ -51,8 +51,8 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed 'Array, Array | Matrix': function (a, b) { const m = matrix(a) - const r = _denseBackwardSubstitution(m, b) - return r.valueOf() + const R = _denseBackwardSubstitution(m, b) + return R.map(r => r.valueOf()) } }) @@ -114,7 +114,7 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed } - return B + return B.map(x => new DenseMatrix({ data: x.map(e=>[e]), size: [rows, 1] })) } function _sparseBackwardSubstitution (m, b) { @@ -188,9 +188,9 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed } } // return vector - return new DenseMatrix({ + return [new DenseMatrix({ data: x, size: [rows, 1] - }) + })] } }) diff --git a/test/unit-tests/function/algebra/solver/lsolve.test.js b/test/unit-tests/function/algebra/solver/lsolve.test.js index b1d9902398..402bcd6e20 100644 --- a/test/unit-tests/function/algebra/solver/lsolve.test.js +++ b/test/unit-tests/function/algebra/solver/lsolve.test.js @@ -17,7 +17,7 @@ describe('lsolve', function () { const x = math.lsolve(m, b) - approx.deepEqual(x, [[1], [1], [1], [1]]) + approx.deepEqual(x, [[[1], [1], [1], [1]]]) }) it('should solve linear system 4 x 4, array and column array', function () { @@ -36,7 +36,7 @@ describe('lsolve', function () { ] const x = math.lsolve(m, b) - approx.deepEqual(x, [[1], [1], [1], [1]]) + approx.deepEqual(x, [[[1], [1], [1], [1]]]) }) it('should solve linear system 4 x 4, matrices', function () { @@ -51,8 +51,8 @@ describe('lsolve', function () { const x = math.lsolve(m, b) - assert(x instanceof math.Matrix) - approx.deepEqual(x, math.matrix([[1], [1], [1], [1]])) + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) }) it('should solve linear system 4 x 4, sparse matrices', function () { @@ -67,8 +67,8 @@ describe('lsolve', function () { const x = math.lsolve(m, b) - assert(x instanceof math.Matrix) - approx.deepEqual(x, math.matrix([[1], [1], [1], [1]])) + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) }) it('should solve linear system 4 x 4, matrix and column matrix', function () { @@ -88,8 +88,8 @@ describe('lsolve', function () { const x = math.lsolve(m, b) - assert(x instanceof math.Matrix) - approx.deepEqual(x, math.matrix([[1], [1], [1], [1]])) + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) }) it('should solve linear system 4 x 4, sparse matrix and column matrix', function () { @@ -109,8 +109,8 @@ describe('lsolve', function () { const x = math.lsolve(m, b) - assert(x instanceof math.Matrix) - approx.deepEqual(x, math.matrix([[1], [1], [1], [1]])) + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) }) it('should throw exception when matrix is singular', function () { diff --git a/test/unit-tests/function/algebra/solver/usolve.test.js b/test/unit-tests/function/algebra/solver/usolve.test.js index b105da6d3a..3797553ed6 100644 --- a/test/unit-tests/function/algebra/solver/usolve.test.js +++ b/test/unit-tests/function/algebra/solver/usolve.test.js @@ -17,7 +17,7 @@ describe('usolve', function () { const x = math.usolve(m, b) - approx.deepEqual(x, [[-1], [-1], [-1], [4]]) + approx.deepEqual(x, [[[-1], [-1], [-1], [4]]]) }) it('should solve linear system 4 x 4, array and column array', function () { @@ -36,7 +36,7 @@ describe('usolve', function () { ] const x = math.usolve(m, b) - approx.deepEqual(x, [[-1], [-1], [-1], [4]]) + approx.deepEqual(x, [[[-1], [-1], [-1], [4]]]) }) it('should solve linear system 4 x 4, matrices', function () { @@ -51,8 +51,8 @@ describe('usolve', function () { const x = math.usolve(m, b) - assert(x instanceof math.Matrix) - approx.deepEqual(x, math.matrix([[-1], [-1], [-1], [4]])) + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) }) it('should solve linear system 4 x 4, sparse matrices', function () { @@ -67,8 +67,8 @@ describe('usolve', function () { const x = math.usolve(m, b) - assert(x instanceof math.Matrix) - approx.deepEqual(x, math.matrix([[-1], [-1], [-1], [4]])) + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) }) it('should solve linear system 4 x 4, matrix and column matrix', function () { @@ -88,8 +88,8 @@ describe('usolve', function () { const x = math.usolve(m, b) - assert(x instanceof math.Matrix) - approx.deepEqual(x, math.matrix([[-1], [-1], [-1], [4]])) + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) }) it('should solve linear system 4 x 4, sparse matrix and column matrix', function () { @@ -109,8 +109,8 @@ describe('usolve', function () { const x = math.usolve(m, b) - assert(x instanceof math.Matrix) - approx.deepEqual(x, math.matrix([[-1], [-1], [-1], [4]])) + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) }) it('should throw exception when matrix is singular', function () { From f3cb1a03b54830b14da7a384e29c5919e9567cc4 Mon Sep 17 00:00:00 2001 From: Michal Grno Date: Tue, 14 Jul 2020 23:32:48 +0200 Subject: [PATCH 5/9] fixed lusolve and its tests, fixed linting issues --- src/function/algebra/solver/lsolve.js | 7 +- src/function/algebra/solver/lusolve.js | 63 ++++++----- src/function/algebra/solver/usolve.js | 7 +- .../function/algebra/solver/lusolve.test.js | 100 +++++++++--------- 4 files changed, 85 insertions(+), 92 deletions(-) diff --git a/src/function/algebra/solver/lsolve.js b/src/function/algebra/solver/lsolve.js index c1998f5e63..398a039ffd 100644 --- a/src/function/algebra/solver/lsolve.js +++ b/src/function/algebra/solver/lsolve.js @@ -28,7 +28,7 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed * * const a = [[-2, 3], [2, 1]] * const b = [11, 9] - * const x = lsolve(a, b) // [[-5.5], [20]] + * const x = lsolve(a, b) // [ [[-5.5], [20]] ] * * See also: * @@ -84,7 +84,6 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed // b[j] -= b[i] * M[j,i] b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i])) } - } else if (!equalScalar(b[i], 0)) { // singular row, nonzero RHS @@ -109,12 +108,10 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed B.push(bNew) } - } - } - return B.map(x => new DenseMatrix({ data: x.map(e=>[e]), size: [rows, 1] })) + return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] })) } function _sparseForwardSubstitution (m, b) { diff --git a/src/function/algebra/solver/lusolve.js b/src/function/algebra/solver/lusolve.js index fe64bf00ae..a6429a71bb 100644 --- a/src/function/algebra/solver/lusolve.js +++ b/src/function/algebra/solver/lusolve.js @@ -19,25 +19,26 @@ export const createLusolve = /* #__PURE__ */ factory(name, dependencies, ({ type /** * Solves the linear system `A * x = b` where `A` is an [n x n] matrix and `b` is a [n] column vector. + * Returns an array of affine-independent vectors that solve the system. * * Syntax: * - * math.lusolve(A, b) // returns column vector with the solution to the linear system A * x = b - * math.lusolve(lup, b) // returns column vector with the solution to the linear system A * x = b, lup = math.lup(A) + * math.lusolve(A, b) // returns the solutions to the linear system A * x = b + * math.lusolve(lup, b) // returns the solutions to the linear system A * x = b, lup = math.lup(A) * * Examples: * * const m = [[1, 0, 0, 0], [0, 2, 0, 0], [0, 0, 3, 0], [0, 0, 0, 4]] * - * const x = math.lusolve(m, [-1, -1, -1, -1]) // x = [[-1], [-0.5], [-1/3], [-0.25]] + * const x = math.lusolve(m, [-1, -1, -1, -1]) // x = [ [[-1], [-0.5], [-1/3], [-0.25]] ] * * const f = math.lup(m) - * const x1 = math.lusolve(f, [-1, -1, -1, -1]) // x1 = [[-1], [-0.5], [-1/3], [-0.25]] - * const x2 = math.lusolve(f, [1, 2, 1, -1]) // x2 = [[1], [1], [1/3], [-0.25]] + * const x1 = math.lusolve(f, [-1, -1, -1, -1]) // x1 = [ [[-1], [-0.5], [-1/3], [-0.25]] ] + * const x2 = math.lusolve(f, [1, 2, 1, -1]) // x2 = [ [[1], [1], [1/3], [-0.25]] ] * * const a = [[-2, 3], [2, 1]] * const b = [11, 9] - * const x = math.lusolve(a, b) // [[2], [5]] + * const x = math.lusolve(a, b) // [ [[2], [5]] ] * * See also: * @@ -48,72 +49,70 @@ export const createLusolve = /* #__PURE__ */ factory(name, dependencies, ({ type * @param {number} [order] The Symbolic Ordering and Analysis order, see slu for details. Matrix must be a SparseMatrix * @param {Number} [threshold] Partial pivoting threshold (1 for partial pivoting), see slu for details. Matrix must be a SparseMatrix. * - * @return {DenseMatrix | Array} Column vector with the solution to the linear system A * x = b + * @return {DenseMatrix[] | Array[]} Column vector with the solution to the linear system A * x = b */ return typed(name, { 'Array, Array | Matrix': function (a, b) { - // convert a to matrix a = matrix(a) - // matrix lup decomposition const d = lup(a) - // solve - const x = _lusolve(d.L, d.U, d.p, null, b) - // convert result to array - return x.valueOf() + const xs = _lusolve(d.L, d.U, d.p, null, b) + return xs.map(x => x.valueOf()) }, 'DenseMatrix, Array | Matrix': function (a, b) { - // matrix lup decomposition const d = lup(a) - // solve return _lusolve(d.L, d.U, d.p, null, b) }, 'SparseMatrix, Array | Matrix': function (a, b) { - // matrix lup decomposition const d = lup(a) - // solve return _lusolve(d.L, d.U, d.p, null, b) }, 'SparseMatrix, Array | Matrix, number, number': function (a, b, order, threshold) { - // matrix lu decomposition const d = slu(a, order, threshold) - // solve return _lusolve(d.L, d.U, d.p, d.q, b) }, 'Object, Array | Matrix': function (d, b) { - // solve return _lusolve(d.L, d.U, d.p, d.q, b) } }) function _toMatrix (a) { - // check it is a matrix if (isMatrix(a)) { return a } - // check array if (isArray(a)) { return matrix(a) } - // throw throw new TypeError('Invalid Matrix LU decomposition') } function _lusolve (l, u, p, q, b) { - // verify L, U, P + // verify decomposition l = _toMatrix(l) u = _toMatrix(u) - // validate matrix and vector - b = solveValidation(l, b, false) + // apply row permutations if needed (b is a DenseMatrix) - if (p) { b._data = csIpvec(p, b._data) } + if (p) { + b = solveValidation(l, b, false) + b._data = csIpvec(p, b._data) + } + // use forward substitution to resolve L * y = b - const y = lsolve(l, b) + const ys = lsolve(l, b) + const xs = [] + // use backward substitution to resolve U * x = y - const x = usolve(u, y) + for (const y of ys) { + xs.push(...usolve(u, y)) + } + // apply column permutations if needed (x is a DenseMatrix) - if (q) { x._data = csIpvec(q, x._data) } - // return solution - return x + if (q) { + for (const x of xs) { + x._data = csIpvec(q, x._data) + } + } + + return xs } }) diff --git a/src/function/algebra/solver/usolve.js b/src/function/algebra/solver/usolve.js index 5cd5b292f5..e6295413d8 100644 --- a/src/function/algebra/solver/usolve.js +++ b/src/function/algebra/solver/usolve.js @@ -28,7 +28,7 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed * * const a = [[-2, 3], [2, 1]] * const b = [11, 9] - * const x = usolve(a, b) // [[8], [9]] + * const x = usolve(a, b) // [ [[8], [9]] ] * * See also: * @@ -84,7 +84,6 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed // b[j] -= b[i] * M[j,i] b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i])) } - } else if (!equalScalar(b[i], 0)) { // singular row, nonzero RHS @@ -109,12 +108,10 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed B.push(bNew) } - } - } - return B.map(x => new DenseMatrix({ data: x.map(e=>[e]), size: [rows, 1] })) + return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] })) } function _sparseBackwardSubstitution (m, b) { diff --git a/test/unit-tests/function/algebra/solver/lusolve.test.js b/test/unit-tests/function/algebra/solver/lusolve.test.js index b76fa4cdd2..525439e351 100644 --- a/test/unit-tests/function/algebra/solver/lusolve.test.js +++ b/test/unit-tests/function/algebra/solver/lusolve.test.js @@ -15,9 +15,9 @@ describe('lusolve', function () { ] const b = [-1, -1, -1, -1] - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - approx.deepEqual(x, [[-1], [-0.5], [-1 / 3], [-0.25]]) + approx.deepEqual(xs, [[[-1], [-0.5], [-1 / 3], [-0.25]]]) }) it('should solve linear system 4 x 4, array and column array', function () { @@ -34,9 +34,9 @@ describe('lusolve', function () { [-1], [-1] ] - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - approx.deepEqual(x, [[-1], [-0.5], [-1 / 3], [-0.25]]) + approx.deepEqual(xs, [[[-1], [-0.5], [-1 / 3], [-0.25]]]) }) it('should solve linear system 4 x 4, matrices', function () { @@ -49,10 +49,10 @@ describe('lusolve', function () { ]) const b = math.matrix([-1, -1, -1, -1]) - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - assert(x instanceof math.Matrix) - approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) + assert(xs[0] instanceof math.Matrix) + approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) }) it('should solve linear system 4 x 4, sparse matrices', function () { @@ -65,10 +65,10 @@ describe('lusolve', function () { ], 'sparse') const b = math.matrix([[-1], [-1], [-1], [-1]], 'sparse') - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - assert(x instanceof math.Matrix) - approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) + assert(xs[0] instanceof math.Matrix) + approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) }) it('should solve linear system 4 x 4, matrix and column matrix', function () { @@ -86,10 +86,10 @@ describe('lusolve', function () { [-1] ]) - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - assert(x instanceof math.Matrix) - approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) + assert(xs[0] instanceof math.Matrix) + approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) }) it('should solve linear system 4 x 4, sparse matrix and column matrix', function () { @@ -107,10 +107,10 @@ describe('lusolve', function () { [-1] ], 'sparse') - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - assert(x instanceof math.Matrix) - approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) + assert(xs[0] instanceof math.Matrix) + approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) }) it('should solve linear system 4 x 4, LUP decomposition (array)', function () { @@ -123,11 +123,11 @@ describe('lusolve', function () { ] const lup = math.lup(m) - const x = math.lusolve(lup, [-1, -1, -1, -1]) - approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) + const xs = math.lusolve(lup, [-1, -1, -1, -1]) + approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) - const y = math.lusolve(lup, [1, 2, 1, -1]) - approx.deepEqual(y, math.matrix([[1], [1], [1 / 3], [-0.25]])) + const ys = math.lusolve(lup, [1, 2, 1, -1]) + approx.deepEqual(ys, [math.matrix([[1], [1], [1 / 3], [-0.25]])]) }) it('should solve linear system 4 x 4, LUP decomposition (matrix)', function () { @@ -140,11 +140,11 @@ describe('lusolve', function () { ]) const lup = math.lup(m) - const x = math.lusolve(lup, [-1, -1, -1, -1]) - approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) + const xs = math.lusolve(lup, [-1, -1, -1, -1]) + approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) - const y = math.lusolve(lup, [1, 2, 1, -1]) - approx.deepEqual(y, math.matrix([[1], [1], [1 / 3], [-0.25]])) + const ys = math.lusolve(lup, [1, 2, 1, -1]) + approx.deepEqual(ys, [math.matrix([[1], [1], [1 / 3], [-0.25]])]) }) it('should solve linear system 4 x 4, LUP decomposition (sparse matrix)', function () { @@ -157,11 +157,11 @@ describe('lusolve', function () { ], 'sparse') const lup = math.lup(m) - const x = math.lusolve(lup, [-1, -1, -1, -1]) - approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) + const xs = math.lusolve(lup, [-1, -1, -1, -1]) + approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) - const y = math.lusolve(lup, [1, 2, 1, -1]) - approx.deepEqual(y, math.matrix([[1], [1], [1 / 3], [-0.25]])) + const ys = math.lusolve(lup, [1, 2, 1, -1]) + approx.deepEqual(ys, [math.matrix([[1], [1], [1 / 3], [-0.25]])]) }) it('should solve linear system 3 x 3, no permutations, arrays', function () { @@ -173,9 +173,9 @@ describe('lusolve', function () { ] const b = [-2, 4, 2] - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - approx.deepEqual(x, [[-5 / 3], [7 / 3], [-1]]) + approx.deepEqual(xs, [[[-5 / 3], [7 / 3], [-1]]]) }) it('should solve linear system 3 x 3, no permutations, matrix', function () { @@ -187,9 +187,9 @@ describe('lusolve', function () { ]) const b = [-2, 4, 2] - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - approx.deepEqual(x, math.matrix([[-5 / 3], [7 / 3], [-1]])) + approx.deepEqual(xs, [math.matrix([[-5 / 3], [7 / 3], [-1]])]) }) it('should solve linear system 3 x 3, no permutations, sparse matrix', function () { @@ -201,9 +201,9 @@ describe('lusolve', function () { ], 'sparse') const b = [-2, 4, 2] - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - approx.deepEqual(x, math.matrix([[-5 / 3], [7 / 3], [-1]])) + approx.deepEqual(xs, [math.matrix([[-5 / 3], [7 / 3], [-1]])]) }) it('should solve linear system 3 x 3, permutations, arrays', function () { @@ -215,9 +215,9 @@ describe('lusolve', function () { ] const b = [4, -2, 2] - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - approx.deepEqual(x, [[-5 / 3], [7 / 3], [-1]]) + approx.deepEqual(xs, [[[-5 / 3], [7 / 3], [-1]]]) }) it('should solve linear system 4 x 4, permutations, matrix - Issue 437', function () { @@ -231,9 +231,9 @@ describe('lusolve', function () { const b = [0.1, 0.2, 0.15, 0.1] - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - approx.deepEqual(x, math.matrix([[0.025], [-0.075], [0], [0.2]])) + approx.deepEqual(xs, [math.matrix([[0.025], [-0.075], [0], [0.2]])]) }) it('should solve linear system 4 x 4, permutations, sparse - Issue 437', function () { @@ -247,9 +247,9 @@ describe('lusolve', function () { const b = [0.1, 0.2, 0.15, 0.1] - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - approx.deepEqual(x, math.matrix([[0.025], [-0.075], [0], [0.2]])) + approx.deepEqual(xs, [math.matrix([[0.025], [-0.075], [0], [0.2]])]) }) it('should solve linear system 3 x 3, permutations, sparse matrix', function () { @@ -261,9 +261,9 @@ describe('lusolve', function () { ], 'sparse') const b = [4, -2, 2] - const x = math.lusolve(m, b) + const xs = math.lusolve(m, b) - approx.deepEqual(x, math.matrix([[-5 / 3], [7 / 3], [-1]])) + approx.deepEqual(xs, [math.matrix([[-5 / 3], [7 / 3], [-1]])]) }) it('should solve linear system 4 x 4, natural ordering (order=0), partial pivoting, sparse matrix', function () { @@ -277,9 +277,9 @@ describe('lusolve', function () { const b = [1.000000, 1.250000, 1.500000, 1.750000] - const x = math.lusolve(m, b, 0, 1) + const xs = math.lusolve(m, b, 0, 1) - approx.deepEqual(x, math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])) + approx.deepEqual(xs, [math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])]) }) it('should solve linear system 4 x 4, amd(A+A\') (order=1), partial pivoting, sparse matrix', function () { @@ -293,9 +293,9 @@ describe('lusolve', function () { const b = [1.000000, 1.250000, 1.500000, 1.750000] - const x = math.lusolve(m, b, 1, 1) + const xs = math.lusolve(m, b, 1, 1) - approx.deepEqual(x, math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])) + approx.deepEqual(xs, [math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])]) }) it('should solve linear system 4 x 4, amd(A\'*A) (order=2), partial pivoting, sparse matrix', function () { @@ -309,9 +309,9 @@ describe('lusolve', function () { const b = [1.000000, 1.250000, 1.500000, 1.750000] - const x = math.lusolve(m, b, 2, 1) + const xs = math.lusolve(m, b, 2, 1) - approx.deepEqual(x, math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])) + approx.deepEqual(xs, [math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])]) }) it('should solve linear system 4 x 4, amd(A\'*A) (order=3), partial pivoting, sparse matrix', function () { @@ -325,9 +325,9 @@ describe('lusolve', function () { const b = [1.000000, 1.250000, 1.500000, 1.750000] - const x = math.lusolve(m, b, 3, 1) + const xs = math.lusolve(m, b, 3, 1) - approx.deepEqual(x, math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])) + approx.deepEqual(xs, [math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])]) }) it('should throw exception when matrix is singular', function () { From 97d9925936690938c797318396923f612039d68a Mon Sep 17 00:00:00 2001 From: Michal Grno Date: Wed, 15 Jul 2020 14:18:42 +0200 Subject: [PATCH 6/9] added tests for usolve&lsolve, try-catch in lusolve --- src/function/algebra/solver/lusolve.js | 7 ++++++- .../function/algebra/solver/lsolve.test.js | 17 +++++++++++++++++ .../function/algebra/solver/usolve.test.js | 17 +++++++++++++++++ 3 files changed, 40 insertions(+), 1 deletion(-) diff --git a/src/function/algebra/solver/lusolve.js b/src/function/algebra/solver/lusolve.js index a6429a71bb..93cfea5d97 100644 --- a/src/function/algebra/solver/lusolve.js +++ b/src/function/algebra/solver/lusolve.js @@ -103,7 +103,12 @@ export const createLusolve = /* #__PURE__ */ factory(name, dependencies, ({ type // use backward substitution to resolve U * x = y for (const y of ys) { - xs.push(...usolve(u, y)) + try { + xs.push(...usolve(u, y)) + } catch(e) { + // no solution found + // not a big deal + } } // apply column permutations if needed (x is a DenseMatrix) diff --git a/test/unit-tests/function/algebra/solver/lsolve.test.js b/test/unit-tests/function/algebra/solver/lsolve.test.js index 402bcd6e20..77f66e3f6d 100644 --- a/test/unit-tests/function/algebra/solver/lsolve.test.js +++ b/test/unit-tests/function/algebra/solver/lsolve.test.js @@ -118,4 +118,21 @@ describe('lsolve', function () { assert.throws(function () { math.lsolve(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) assert.throws(function () { math.lsolve(math.matrix([[1, 1], [0, 0]], 'sparse'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) }) + + it('should solve systems with singular matrices', function () { + approx.deepEqual( + math.lsolve([[2, 0, 0], [1, 0, 0], [-1, 1, 1]], [4, 2, 1]), + [[[2], [0], [3]], [[2], [1], [2]]] + ) + + approx.deepEqual( + math.lsolve([[0, 0, 0], [1, 1, 0], [2, 1, 0]], [0, 2, 2]), + [[[0], [2], [0]], [[0], [2], [1]]] + ) + + approx.deepEqual( + math.lsolve([[0, 0, 0], [1, 1, 0], [1, 1, 0]], [0, 2, 2]), + [[[0], [2], [0]], [[1], [1], [0]], [[0], [2], [1]]] + ) + }) }) diff --git a/test/unit-tests/function/algebra/solver/usolve.test.js b/test/unit-tests/function/algebra/solver/usolve.test.js index 3797553ed6..3f579acc80 100644 --- a/test/unit-tests/function/algebra/solver/usolve.test.js +++ b/test/unit-tests/function/algebra/solver/usolve.test.js @@ -118,4 +118,21 @@ describe('usolve', function () { assert.throws(function () { math.usolve(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) assert.throws(function () { math.usolve(math.matrix([[1, 1], [0, 0]], 'sparse'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) }) + + it('should solve systems with singular matrices', function () { + approx.deepEqual( + math.usolve([[1, 1, -1], [0, 0, 1], [0, 0, 2]], [1, 2, 4]), + [[[3], [0], [2]], [[2], [1], [2]]] + ) + + approx.deepEqual( + math.usolve([[0, 1, 2], [0, 1, 1], [0, 0, 0]], [2, 2, 0]), + [[[0], [2], [0]], [[1], [2], [0]]] + ) + + approx.deepEqual( + math.usolve([[0, 1, 1], [0, 1, 1], [0, 0, 0]], [2, 2, 0]), + [[[0], [2], [0]], [[0], [1], [1]], [[1], [2], [0]]] + ) + }) }) From b4ec1d638394ca207b4068b916e9229f3d6f8017 Mon Sep 17 00:00:00 2001 From: Michal Grno Date: Sun, 19 Jul 2020 23:07:17 +0200 Subject: [PATCH 7/9] put changes into separate files (u-/lsolveAll), revert changes to u-, l- and lusolve --- src/factoriesAny.js | 2 + src/function/algebra/solver/lsolve.js | 149 +++++++------- src/function/algebra/solver/lsolveAll.js | 187 +++++++++++++++++ src/function/algebra/solver/lusolve.js | 42 ++-- src/function/algebra/solver/usolve.js | 125 +++++------- src/function/algebra/solver/usolveAll.js | 193 ++++++++++++++++++ .../function/algebra/solver/lsolve.test.js | 37 +--- .../function/algebra/solver/lsolveAll.test.js | 138 +++++++++++++ .../function/algebra/solver/lusolve.test.js | 100 ++++----- .../function/algebra/solver/usolve.test.js | 37 +--- .../function/algebra/solver/usolveAll.test.js | 138 +++++++++++++ 11 files changed, 867 insertions(+), 281 deletions(-) create mode 100644 src/function/algebra/solver/lsolveAll.js create mode 100644 src/function/algebra/solver/usolveAll.js create mode 100644 test/unit-tests/function/algebra/solver/lsolveAll.test.js create mode 100644 test/unit-tests/function/algebra/solver/usolveAll.test.js diff --git a/src/factoriesAny.js b/src/factoriesAny.js index de60e761aa..a87428fd74 100644 --- a/src/factoriesAny.js +++ b/src/factoriesAny.js @@ -103,6 +103,8 @@ export { createDotPow } from './function/arithmetic/dotPow' export { createDotDivide } from './function/arithmetic/dotDivide' export { createLsolve } from './function/algebra/solver/lsolve' export { createUsolve } from './function/algebra/solver/usolve' +export { createLsolveAll } from './function/algebra/solver/lsolveAll' +export { createUsolveAll } from './function/algebra/solver/usolveAll' export { createLeftShift } from './function/bitwise/leftShift' export { createRightArithShift } from './function/bitwise/rightArithShift' export { createRightLogShift } from './function/bitwise/rightLogShift' diff --git a/src/function/algebra/solver/lsolve.js b/src/function/algebra/solver/lsolve.js index 398a039ffd..64ba9bc331 100644 --- a/src/function/algebra/solver/lsolve.js +++ b/src/function/algebra/solver/lsolve.js @@ -28,7 +28,7 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed * * const a = [[-2, 3], [2, 1]] * const b = [11, 9] - * const x = lsolve(a, b) // [ [[-5.5], [20]] ] + * const x = lsolve(a, b) // [[-5.5], [20]] * * See also: * @@ -37,7 +37,7 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed * @param {Matrix, Array} L A N x N matrix or array (L) * @param {Matrix, Array} b A column vector with the b values * - * @return {DenseMatrix[] | Array[]} An array of affine-independent column vectors (x) that solve the linear system + * @return {DenseMatrix | Array} A column vector with the linear system solution (x) */ return typed(name, { @@ -51,137 +51,126 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed 'Array, Array | Matrix': function (a, b) { const m = matrix(a) - const R = _denseForwardSubstitution(m, b) - return R.map(r => r.valueOf()) + const r = _denseForwardSubstitution(m, b) + return r.valueOf() } }) - function _denseForwardSubstitution (m, b_) { - // the algorithm is derived from - // https://www.overleaf.com/project/5e6c87c554a3190001a3fc93 - - // array of right-hand sides - const B = [solveValidation(m, b_, true)._data.map(e => e[0])] + function _denseForwardSubstitution (m, b) { + // validate matrix and vector, return copy of column vector b + b = solveValidation(m, b, true) + const bdata = b._data - const M = m._data const rows = m._size[0] const columns = m._size[1] - // loop columns - for (let i = 0; i < columns; i++) { - let L = B.length + // result + const x = [] - // loop right-hand sides - for (let k = 0; k < L; k++) { - const b = B[k] + const mdata = m._data - if (!equalScalar(M[i][i], 0)) { - // non-singular row + // loop columns + for (let j = 0; j < columns; j++) { + const bj = bdata[j][0] || 0 + let xj - b[i] = divideScalar(b[i], M[i][i]) + if (!equalScalar(bj, 0)) { + // non-degenerate row, find solution - for (let j = i + 1; j < columns; j++) { - // b[j] -= b[i] * M[j,i] - b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i])) - } - } else if (!equalScalar(b[i], 0)) { - // singular row, nonzero RHS - - if (k === 0) { - // There is no valid solution - throw new Error('Linear system cannot be solved since matrix is singular') - } else { - // This RHS is invalid but other solutions may still exist - B.splice(k, 1) - k -= 1 - L -= 1 - } - } else if (k === 0) { - // singular row, RHS is zero + const vjj = mdata[j][j] - const bNew = [...b] - bNew[i] = 1 + if (equalScalar(vjj, 0)) { + throw new Error('Linear system cannot be solved since matrix is singular') + } - for (let j = i + 1; j < columns; j++) { - bNew[j] = subtract(bNew[j], M[j][i]) - } + xj = divideScalar(bj, vjj) - B.push(bNew) + // loop rows + for (let i = j + 1; i < rows; i++) { + bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, mdata[i][j]))] } + } else { + // degenerate row, we can choose any value + xj = 0 } + + x[j] = [xj] } - return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] })) + return new DenseMatrix({ + data: x, + size: [rows, 1] + }) } function _sparseForwardSubstitution (m, b) { // validate matrix and vector, return copy of column vector b b = solveValidation(m, b, true) - // column vector data + const bdata = b._data - // rows & columns + const rows = m._size[0] const columns = m._size[1] - // matrix arrays + const values = m._values const index = m._index const ptr = m._ptr - // vars - let i, k + // result const x = [] - // forward solve m * x = b, loop columns + + // loop columns for (let j = 0; j < columns; j++) { - // b[j] const bj = bdata[j][0] || 0 - // forward substitution (outer product) avoids inner looping when bj === 0 + if (!equalScalar(bj, 0)) { - // value @ [j, j] + // non-degenerate row, find solution + let vjj = 0 - // lower triangular matrix values & index (column j) - const jvalues = [] - const jindex = [] - // last index in column - let l = ptr[j + 1] - // values in column, find value @ [j, j] - for (k = ptr[j]; k < l; k++) { - // row - i = index[k] + // matrix values & indices (column j) + const jValues = [] + const jIndices = [] + + // first and last index in the column + const firstIndex = ptr[j] + const lastIndex = ptr[j + 1] + + // values in column, find value at [j, j] + for (let k = firstIndex; k < lastIndex; k++) { + const i = index[k] + // check row (rows are not sorted!) if (i === j) { - // update vjj vjj = values[k] } else if (i > j) { // store lower triangular - jvalues.push(values[k]) - jindex.push(i) + jValues.push(values[k]) + jIndices.push(i) } } - // at this point we must have a value @ [j, j] + + // at this point we must have a value in vjj if (equalScalar(vjj, 0)) { - // system cannot be solved, there is no value @ [j, j] throw new Error('Linear system cannot be solved since matrix is singular') } - // calculate xj + const xj = divideScalar(bj, vjj) - // loop lower triangular - for (k = 0, l = jindex.length; k < l; k++) { - // row - i = jindex[k] - // update copy of b - bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, jvalues[k]))] + + for (let k = 0, l = jIndices.length; k < l; k++) { + const i = jIndices[k] + bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, jValues[k]))] } - // update x + x[j] = [xj] } else { - // update x + // degenerate row, we can choose any value x[j] = [0] } } - // return vector - return [new DenseMatrix({ + + return new DenseMatrix({ data: x, size: [rows, 1] - })] + }) } }) diff --git a/src/function/algebra/solver/lsolveAll.js b/src/function/algebra/solver/lsolveAll.js new file mode 100644 index 0000000000..9be4a7231b --- /dev/null +++ b/src/function/algebra/solver/lsolveAll.js @@ -0,0 +1,187 @@ +import { factory } from '../../../utils/factory' +import { createSolveValidation } from './utils/solveValidation' + +const name = 'lsolveAll' +const dependencies = [ + 'typed', + 'matrix', + 'divideScalar', + 'multiplyScalar', + 'subtract', + 'equalScalar', + 'DenseMatrix' +] + +export const createLsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, divideScalar, multiplyScalar, subtract, equalScalar, DenseMatrix }) => { + const solveValidation = createSolveValidation({ DenseMatrix }) + + /** + * Solves the linear equation system by forwards substitution. Matrix must be a lower triangular matrix. + * + * `L * x = b` + * + * Syntax: + * + * math.lsolve(L, b) + * + * Examples: + * + * const a = [[-2, 3], [2, 1]] + * const b = [11, 9] + * const x = lsolve(a, b) // [ [[-5.5], [20]] ] + * + * See also: + * + * lup, slu, usolve, lusolve + * + * @param {Matrix, Array} L A N x N matrix or array (L) + * @param {Matrix, Array} b A column vector with the b values + * + * @return {DenseMatrix[] | Array[]} An array of affine-independent column vectors (x) that solve the linear system + */ + return typed(name, { + + 'SparseMatrix, Array | Matrix': function (m, b) { + return _sparseForwardSubstitution(m, b) + }, + + 'DenseMatrix, Array | Matrix': function (m, b) { + return _denseForwardSubstitution(m, b) + }, + + 'Array, Array | Matrix': function (a, b) { + const m = matrix(a) + const R = _denseForwardSubstitution(m, b) + return R.map(r => r.valueOf()) + } + }) + + function _denseForwardSubstitution (m, b_) { + // the algorithm is derived from + // https://www.overleaf.com/project/5e6c87c554a3190001a3fc93 + + // array of right-hand sides + const B = [solveValidation(m, b_, true)._data.map(e => e[0])] + + const M = m._data + const rows = m._size[0] + const columns = m._size[1] + + // loop columns + for (let i = 0; i < columns; i++) { + let L = B.length + + // loop right-hand sides + for (let k = 0; k < L; k++) { + const b = B[k] + + if (!equalScalar(M[i][i], 0)) { + // non-singular row + + b[i] = divideScalar(b[i], M[i][i]) + + for (let j = i + 1; j < columns; j++) { + // b[j] -= b[i] * M[j,i] + b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i])) + } + } else if (!equalScalar(b[i], 0)) { + // singular row, nonzero RHS + + if (k === 0) { + // There is no valid solution + throw new Error('Linear system cannot be solved since matrix is singular') + } else { + // This RHS is invalid but other solutions may still exist + B.splice(k, 1) + k -= 1 + L -= 1 + } + } else if (k === 0) { + // singular row, RHS is zero + + const bNew = [...b] + bNew[i] = 1 + + for (let j = i + 1; j < columns; j++) { + bNew[j] = subtract(bNew[j], M[j][i]) + } + + B.push(bNew) + } + } + } + + return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] })) + } + + function _sparseForwardSubstitution (m, b) { + // validate matrix and vector, return copy of column vector b + b = solveValidation(m, b, true) + // column vector data + const bdata = b._data + // rows & columns + const rows = m._size[0] + const columns = m._size[1] + // matrix arrays + const values = m._values + const index = m._index + const ptr = m._ptr + // vars + let i, k + // result + const x = [] + // forward solve m * x = b, loop columns + for (let j = 0; j < columns; j++) { + // b[j] + const bj = bdata[j][0] || 0 + // forward substitution (outer product) avoids inner looping when bj === 0 + if (!equalScalar(bj, 0)) { + // value @ [j, j] + let vjj = 0 + // lower triangular matrix values & index (column j) + const jvalues = [] + const jindex = [] + // last index in column + let l = ptr[j + 1] + // values in column, find value @ [j, j] + for (k = ptr[j]; k < l; k++) { + // row + i = index[k] + // check row (rows are not sorted!) + if (i === j) { + // update vjj + vjj = values[k] + } else if (i > j) { + // store lower triangular + jvalues.push(values[k]) + jindex.push(i) + } + } + // at this point we must have a value @ [j, j] + if (equalScalar(vjj, 0)) { + // system cannot be solved, there is no value @ [j, j] + throw new Error('Linear system cannot be solved since matrix is singular') + } + // calculate xj + const xj = divideScalar(bj, vjj) + // loop lower triangular + for (k = 0, l = jindex.length; k < l; k++) { + // row + i = jindex[k] + // update copy of b + bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, jvalues[k]))] + } + // update x + x[j] = [xj] + } else { + // update x + x[j] = [0] + } + } + // return vector + return [new DenseMatrix({ + data: x, + size: [rows, 1] + })] + } +}) diff --git a/src/function/algebra/solver/lusolve.js b/src/function/algebra/solver/lusolve.js index 93cfea5d97..8bebb86051 100644 --- a/src/function/algebra/solver/lusolve.js +++ b/src/function/algebra/solver/lusolve.js @@ -19,26 +19,25 @@ export const createLusolve = /* #__PURE__ */ factory(name, dependencies, ({ type /** * Solves the linear system `A * x = b` where `A` is an [n x n] matrix and `b` is a [n] column vector. - * Returns an array of affine-independent vectors that solve the system. * * Syntax: * - * math.lusolve(A, b) // returns the solutions to the linear system A * x = b - * math.lusolve(lup, b) // returns the solutions to the linear system A * x = b, lup = math.lup(A) + * math.lusolve(A, b) // returns column vector with the solution to the linear system A * x = b + * math.lusolve(lup, b) // returns column vector with the solution to the linear system A * x = b, lup = math.lup(A) * * Examples: * * const m = [[1, 0, 0, 0], [0, 2, 0, 0], [0, 0, 3, 0], [0, 0, 0, 4]] * - * const x = math.lusolve(m, [-1, -1, -1, -1]) // x = [ [[-1], [-0.5], [-1/3], [-0.25]] ] + * const x = math.lusolve(m, [-1, -1, -1, -1]) // x = [[-1], [-0.5], [-1/3], [-0.25]] * * const f = math.lup(m) - * const x1 = math.lusolve(f, [-1, -1, -1, -1]) // x1 = [ [[-1], [-0.5], [-1/3], [-0.25]] ] - * const x2 = math.lusolve(f, [1, 2, 1, -1]) // x2 = [ [[1], [1], [1/3], [-0.25]] ] + * const x1 = math.lusolve(f, [-1, -1, -1, -1]) // x1 = [[-1], [-0.5], [-1/3], [-0.25]] + * const x2 = math.lusolve(f, [1, 2, 1, -1]) // x2 = [[1], [1], [1/3], [-0.25]] * * const a = [[-2, 3], [2, 1]] * const b = [11, 9] - * const x = math.lusolve(a, b) // [ [[2], [5]] ] + * const x = math.lusolve(a, b) // [[2], [5]] * * See also: * @@ -49,15 +48,15 @@ export const createLusolve = /* #__PURE__ */ factory(name, dependencies, ({ type * @param {number} [order] The Symbolic Ordering and Analysis order, see slu for details. Matrix must be a SparseMatrix * @param {Number} [threshold] Partial pivoting threshold (1 for partial pivoting), see slu for details. Matrix must be a SparseMatrix. * - * @return {DenseMatrix[] | Array[]} Column vector with the solution to the linear system A * x = b + * @return {DenseMatrix | Array} Column vector with the solution to the linear system A * x = b */ return typed(name, { 'Array, Array | Matrix': function (a, b) { a = matrix(a) const d = lup(a) - const xs = _lusolve(d.L, d.U, d.p, null, b) - return xs.map(x => x.valueOf()) + const x = _lusolve(d.L, d.U, d.p, null, b) + return x.valueOf() }, 'DenseMatrix, Array | Matrix': function (a, b) { @@ -93,31 +92,18 @@ export const createLusolve = /* #__PURE__ */ factory(name, dependencies, ({ type // apply row permutations if needed (b is a DenseMatrix) if (p) { - b = solveValidation(l, b, false) + b = solveValidation(l, b, true) b._data = csIpvec(p, b._data) } // use forward substitution to resolve L * y = b - const ys = lsolve(l, b) - const xs = [] - + const y = lsolve(l, b) // use backward substitution to resolve U * x = y - for (const y of ys) { - try { - xs.push(...usolve(u, y)) - } catch(e) { - // no solution found - // not a big deal - } - } + const x = usolve(u, y) // apply column permutations if needed (x is a DenseMatrix) - if (q) { - for (const x of xs) { - x._data = csIpvec(q, x._data) - } - } + if (q) { x._data = csIpvec(q, x._data) } - return xs + return x } }) diff --git a/src/function/algebra/solver/usolve.js b/src/function/algebra/solver/usolve.js index e6295413d8..aea9e6cc0d 100644 --- a/src/function/algebra/solver/usolve.js +++ b/src/function/algebra/solver/usolve.js @@ -28,7 +28,7 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed * * const a = [[-2, 3], [2, 1]] * const b = [11, 9] - * const x = usolve(a, b) // [ [[8], [9]] ] + * const x = usolve(a, b) // [[8], [9]] * * See also: * @@ -37,7 +37,7 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed * @param {Matrix, Array} U A N x N matrix or array (U) * @param {Matrix, Array} b A column vector with the b values * - * @return {DenseMatrix[] | Array[]} An array of affine-independent column vectors (x) that solve the linear system + * @return {DenseMatrix | Array} A column vector with the linear system solution (x) */ return typed(name, { @@ -51,67 +51,59 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed 'Array, Array | Matrix': function (a, b) { const m = matrix(a) - const R = _denseBackwardSubstitution(m, b) - return R.map(r => r.valueOf()) + const r = _denseBackwardSubstitution(m, b) + return r.valueOf() } }) - function _denseBackwardSubstitution (m, b_) { - // the algorithm is derived from - // https://www.overleaf.com/project/5e6c87c554a3190001a3fc93 + function _denseBackwardSubstitution (m, b) { + // make b into a column vector + b = solveValidation(m, b, true) - // array of right-hand sides - const B = [solveValidation(m, b_, true)._data.map(e => e[0])] + const bdata = b._data - const M = m._data const rows = m._size[0] const columns = m._size[1] - // loop columns backwards - for (let i = columns - 1; i >= 0; i--) { - let L = B.length - - // loop right-hand sides - for (let k = 0; k < L; k++) { - const b = B[k] - - if (!equalScalar(M[i][i], 0)) { - // non-singular row + // result + const x = [] - b[i] = divideScalar(b[i], M[i][i]) + const mdata = m._data + // loop columns backwards + for (let j = columns - 1; j >= 0; j--) { + // b[j] + const bj = bdata[j][0] || 0 + // x[j] + let xj - for (let j = i - 1; j >= 0; j--) { - // b[j] -= b[i] * M[j,i] - b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i])) - } - } else if (!equalScalar(b[i], 0)) { - // singular row, nonzero RHS - - if (k === 0) { - // There is no valid solution - throw new Error('Linear system cannot be solved since matrix is singular') - } else { - // This RHS is invalid but other solutions may still exist - B.splice(k, 1) - k -= 1 - L -= 1 - } - } else if (k === 0) { - // singular row, RHS is zero + if (!equalScalar(bj, 0)) { + // value at [j, j] + const vjj = mdata[j][j] - const bNew = [...b] - bNew[i] = 1 + if (equalScalar(vjj, 0)) { + // system cannot be solved + throw new Error('Linear system cannot be solved since matrix is singular') + } - for (let j = i - 1; j >= 0; j--) { - bNew[j] = subtract(bNew[j], M[j][i]) - } + xj = divideScalar(bj, vjj) - B.push(bNew) + // loop rows + for (let i = j - 1; i >= 0; i--) { + // update copy of b + bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, mdata[i][j]))] } + } else { + // zero value at j + xj = 0 } + // update x + x[j] = [xj] } - return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] })) + return new DenseMatrix({ + data: x, + size: [rows, 1] + }) } function _sparseBackwardSubstitution (m, b) { @@ -127,18 +119,16 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed const index = m._index const ptr = m._ptr - let i, k - // result const x = [] // loop columns backwards for (let j = columns - 1; j >= 0; j--) { - // b[j] const bj = bdata[j][0] || 0 if (!equalScalar(bj, 0)) { - // value at [j, j] + // non-degenerate row, find solution + let vjj = 0 // upper triangular matrix values & index (column j) @@ -147,15 +137,14 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed // first & last indeces in column const firstIndex = ptr[j] - let lastIndex = ptr[j + 1] + const lastIndex = ptr[j + 1] // values in column, find value at [j, j], loop backwards - for (k = lastIndex - 1; k >= firstIndex; k--) { - // row - i = index[k] - // check row + for (let k = lastIndex - 1; k >= firstIndex; k--) { + const i = index[k] + + // check row (rows are not sorted!) if (i === j) { - // update vjj vjj = values[k] } else if (i < j) { // store upper triangular @@ -163,31 +152,29 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed jIndices.push(i) } } - // at this point we must have a value at [j, j] + + // at this point we must have a value in vjj if (equalScalar(vjj, 0)) { - // system cannot be solved, there is no value at [j, j] throw new Error('Linear system cannot be solved since matrix is singular') } - // calculate xj + const xj = divideScalar(bj, vjj) - // loop upper triangular - for (k = 0, lastIndex = jIndices.length; k < lastIndex; k++) { - // row - i = jIndices[k] - // update copy of b + + for (let k = 0, lastIndex = jIndices.length; k < lastIndex; k++) { + const i = jIndices[k] bdata[i] = [subtract(bdata[i][0], multiplyScalar(xj, jValues[k]))] } - // update x + x[j] = [xj] } else { - // update x + // degenerate row, we can choose any value x[j] = [0] } } - // return vector - return [new DenseMatrix({ + + return new DenseMatrix({ data: x, size: [rows, 1] - })] + }) } }) diff --git a/src/function/algebra/solver/usolveAll.js b/src/function/algebra/solver/usolveAll.js new file mode 100644 index 0000000000..c07fafe723 --- /dev/null +++ b/src/function/algebra/solver/usolveAll.js @@ -0,0 +1,193 @@ +import { factory } from '../../../utils/factory' +import { createSolveValidation } from './utils/solveValidation' + +const name = 'usolveAll' +const dependencies = [ + 'typed', + 'matrix', + 'divideScalar', + 'multiplyScalar', + 'subtract', + 'equalScalar', + 'DenseMatrix' +] + +export const createUsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, divideScalar, multiplyScalar, subtract, equalScalar, DenseMatrix }) => { + const solveValidation = createSolveValidation({ DenseMatrix }) + + /** + * Solves the linear equation system by backward substitution. Matrix must be an upper triangular matrix. + * + * `U * x = b` + * + * Syntax: + * + * math.usolve(U, b) + * + * Examples: + * + * const a = [[-2, 3], [2, 1]] + * const b = [11, 9] + * const x = usolve(a, b) // [ [[8], [9]] ] + * + * See also: + * + * lup, slu, usolve, lusolve + * + * @param {Matrix, Array} U A N x N matrix or array (U) + * @param {Matrix, Array} b A column vector with the b values + * + * @return {DenseMatrix[] | Array[]} An array of affine-independent column vectors (x) that solve the linear system + */ + return typed(name, { + + 'SparseMatrix, Array | Matrix': function (m, b) { + return _sparseBackwardSubstitution(m, b) + }, + + 'DenseMatrix, Array | Matrix': function (m, b) { + return _denseBackwardSubstitution(m, b) + }, + + 'Array, Array | Matrix': function (a, b) { + const m = matrix(a) + const R = _denseBackwardSubstitution(m, b) + return R.map(r => r.valueOf()) + } + }) + + function _denseBackwardSubstitution (m, b_) { + // the algorithm is derived from + // https://www.overleaf.com/project/5e6c87c554a3190001a3fc93 + + // array of right-hand sides + const B = [solveValidation(m, b_, true)._data.map(e => e[0])] + + const M = m._data + const rows = m._size[0] + const columns = m._size[1] + + // loop columns backwards + for (let i = columns - 1; i >= 0; i--) { + let L = B.length + + // loop right-hand sides + for (let k = 0; k < L; k++) { + const b = B[k] + + if (!equalScalar(M[i][i], 0)) { + // non-singular row + + b[i] = divideScalar(b[i], M[i][i]) + + for (let j = i - 1; j >= 0; j--) { + // b[j] -= b[i] * M[j,i] + b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i])) + } + } else if (!equalScalar(b[i], 0)) { + // singular row, nonzero RHS + + if (k === 0) { + // There is no valid solution + throw new Error('Linear system cannot be solved since matrix is singular') + } else { + // This RHS is invalid but other solutions may still exist + B.splice(k, 1) + k -= 1 + L -= 1 + } + } else if (k === 0) { + // singular row, RHS is zero + + const bNew = [...b] + bNew[i] = 1 + + for (let j = i - 1; j >= 0; j--) { + bNew[j] = subtract(bNew[j], M[j][i]) + } + + B.push(bNew) + } + } + } + + return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] })) + } + + function _sparseBackwardSubstitution (m, b) { + // make b into a column vector + b = solveValidation(m, b, true) + + const bdata = b._data + + const rows = m._size[0] + const columns = m._size[1] + + const values = m._values + const index = m._index + const ptr = m._ptr + + let i, k + + // result + const x = [] + + // loop columns backwards + for (let j = columns - 1; j >= 0; j--) { + // b[j] + const bj = bdata[j][0] || 0 + + if (!equalScalar(bj, 0)) { + // value at [j, j] + let vjj = 0 + + // upper triangular matrix values & index (column j) + const jValues = [] + const jIndices = [] + + // first & last indeces in column + const firstIndex = ptr[j] + let lastIndex = ptr[j + 1] + + // values in column, find value at [j, j], loop backwards + for (k = lastIndex - 1; k >= firstIndex; k--) { + // row + i = index[k] + // check row + if (i === j) { + // update vjj + vjj = values[k] + } else if (i < j) { + // store upper triangular + jValues.push(values[k]) + jIndices.push(i) + } + } + // at this point we must have a value at [j, j] + if (equalScalar(vjj, 0)) { + // system cannot be solved, there is no value at [j, j] + throw new Error('Linear system cannot be solved since matrix is singular') + } + // calculate xj + const xj = divideScalar(bj, vjj) + // loop upper triangular + for (k = 0, lastIndex = jIndices.length; k < lastIndex; k++) { + // row + i = jIndices[k] + // update copy of b + bdata[i] = [subtract(bdata[i][0], multiplyScalar(xj, jValues[k]))] + } + // update x + x[j] = [xj] + } else { + // update x + x[j] = [0] + } + } + // return vector + return [new DenseMatrix({ + data: x, + size: [rows, 1] + })] + } +}) diff --git a/test/unit-tests/function/algebra/solver/lsolve.test.js b/test/unit-tests/function/algebra/solver/lsolve.test.js index 77f66e3f6d..b1d9902398 100644 --- a/test/unit-tests/function/algebra/solver/lsolve.test.js +++ b/test/unit-tests/function/algebra/solver/lsolve.test.js @@ -17,7 +17,7 @@ describe('lsolve', function () { const x = math.lsolve(m, b) - approx.deepEqual(x, [[[1], [1], [1], [1]]]) + approx.deepEqual(x, [[1], [1], [1], [1]]) }) it('should solve linear system 4 x 4, array and column array', function () { @@ -36,7 +36,7 @@ describe('lsolve', function () { ] const x = math.lsolve(m, b) - approx.deepEqual(x, [[[1], [1], [1], [1]]]) + approx.deepEqual(x, [[1], [1], [1], [1]]) }) it('should solve linear system 4 x 4, matrices', function () { @@ -51,8 +51,8 @@ describe('lsolve', function () { const x = math.lsolve(m, b) - assert(x[0] instanceof math.Matrix) - approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) + assert(x instanceof math.Matrix) + approx.deepEqual(x, math.matrix([[1], [1], [1], [1]])) }) it('should solve linear system 4 x 4, sparse matrices', function () { @@ -67,8 +67,8 @@ describe('lsolve', function () { const x = math.lsolve(m, b) - assert(x[0] instanceof math.Matrix) - approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) + assert(x instanceof math.Matrix) + approx.deepEqual(x, math.matrix([[1], [1], [1], [1]])) }) it('should solve linear system 4 x 4, matrix and column matrix', function () { @@ -88,8 +88,8 @@ describe('lsolve', function () { const x = math.lsolve(m, b) - assert(x[0] instanceof math.Matrix) - approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) + assert(x instanceof math.Matrix) + approx.deepEqual(x, math.matrix([[1], [1], [1], [1]])) }) it('should solve linear system 4 x 4, sparse matrix and column matrix', function () { @@ -109,8 +109,8 @@ describe('lsolve', function () { const x = math.lsolve(m, b) - assert(x[0] instanceof math.Matrix) - approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) + assert(x instanceof math.Matrix) + approx.deepEqual(x, math.matrix([[1], [1], [1], [1]])) }) it('should throw exception when matrix is singular', function () { @@ -118,21 +118,4 @@ describe('lsolve', function () { assert.throws(function () { math.lsolve(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) assert.throws(function () { math.lsolve(math.matrix([[1, 1], [0, 0]], 'sparse'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) }) - - it('should solve systems with singular matrices', function () { - approx.deepEqual( - math.lsolve([[2, 0, 0], [1, 0, 0], [-1, 1, 1]], [4, 2, 1]), - [[[2], [0], [3]], [[2], [1], [2]]] - ) - - approx.deepEqual( - math.lsolve([[0, 0, 0], [1, 1, 0], [2, 1, 0]], [0, 2, 2]), - [[[0], [2], [0]], [[0], [2], [1]]] - ) - - approx.deepEqual( - math.lsolve([[0, 0, 0], [1, 1, 0], [1, 1, 0]], [0, 2, 2]), - [[[0], [2], [0]], [[1], [1], [0]], [[0], [2], [1]]] - ) - }) }) diff --git a/test/unit-tests/function/algebra/solver/lsolveAll.test.js b/test/unit-tests/function/algebra/solver/lsolveAll.test.js new file mode 100644 index 0000000000..0415a037d6 --- /dev/null +++ b/test/unit-tests/function/algebra/solver/lsolveAll.test.js @@ -0,0 +1,138 @@ +// test lsolveAll +import assert from 'assert' + +import approx from '../../../../../tools/approx' +import math from '../../../../../src/bundleAny' + +describe('lsolveAll', function () { + it('should solve linear system 4 x 4, arrays', function () { + const m = + [ + [1, 0, 0, 0], + [1, 1, 0, 0], + [1, 1, 1, 0], + [1, 1, 1, 1] + ] + const b = [1, 2, 3, 4] + + const x = math.lsolveAll(m, b) + + approx.deepEqual(x, [[[1], [1], [1], [1]]]) + }) + + it('should solve linear system 4 x 4, array and column array', function () { + const m = + [ + [1, 0, 0, 0], + [1, 1, 0, 0], + [1, 1, 1, 0], + [1, 1, 1, 1] + ] + const b = [ + [1], + [2], + [3], + [4] + ] + const x = math.lsolveAll(m, b) + + approx.deepEqual(x, [[[1], [1], [1], [1]]]) + }) + + it('should solve linear system 4 x 4, matrices', function () { + const m = math.matrix( + [ + [1, 0, 0, 0], + [1, 1, 0, 0], + [1, 1, 1, 0], + [1, 1, 1, 1] + ]) + const b = math.matrix([1, 2, 3, 4]) + + const x = math.lsolveAll(m, b) + + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) + }) + + it('should solve linear system 4 x 4, sparse matrices', function () { + const m = math.sparse( + [ + [1, 0, 0, 0], + [1, 1, 0, 0], + [1, 1, 1, 0], + [1, 1, 1, 1] + ]) + const b = math.matrix([[1], [2], [3], [4]], 'sparse') + + const x = math.lsolveAll(m, b) + + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) + }) + + it('should solve linear system 4 x 4, matrix and column matrix', function () { + const m = math.matrix( + [ + [1, 0, 0, 0], + [1, 1, 0, 0], + [1, 1, 1, 0], + [1, 1, 1, 1] + ]) + const b = math.matrix([ + [1], + [2], + [3], + [4] + ]) + + const x = math.lsolveAll(m, b) + + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) + }) + + it('should solve linear system 4 x 4, sparse matrix and column matrix', function () { + const m = math.matrix( + [ + [1, 0, 0, 0], + [1, 1, 0, 0], + [1, 1, 1, 0], + [1, 1, 1, 1] + ], 'sparse') + const b = math.matrix([ + [1], + [2], + [3], + [4] + ], 'sparse') + + const x = math.lsolveAll(m, b) + + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) + }) + + it('should throw exception when matrix is singular', function () { + assert.throws(function () { math.lsolveAll([[1, 1], [0, 0]], [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) + assert.throws(function () { math.lsolveAll(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) + assert.throws(function () { math.lsolveAll(math.matrix([[1, 1], [0, 0]], 'sparse'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) + }) + + it('should solve systems with singular matrices', function () { + approx.deepEqual( + math.lsolveAll([[2, 0, 0], [1, 0, 0], [-1, 1, 1]], [4, 2, 1]), + [[[2], [0], [3]], [[2], [1], [2]]] + ) + + approx.deepEqual( + math.lsolveAll([[0, 0, 0], [1, 1, 0], [2, 1, 0]], [0, 2, 2]), + [[[0], [2], [0]], [[0], [2], [1]]] + ) + + approx.deepEqual( + math.lsolveAll([[0, 0, 0], [1, 1, 0], [1, 1, 0]], [0, 2, 2]), + [[[0], [2], [0]], [[1], [1], [0]], [[0], [2], [1]]] + ) + }) +}) diff --git a/test/unit-tests/function/algebra/solver/lusolve.test.js b/test/unit-tests/function/algebra/solver/lusolve.test.js index 525439e351..b76fa4cdd2 100644 --- a/test/unit-tests/function/algebra/solver/lusolve.test.js +++ b/test/unit-tests/function/algebra/solver/lusolve.test.js @@ -15,9 +15,9 @@ describe('lusolve', function () { ] const b = [-1, -1, -1, -1] - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - approx.deepEqual(xs, [[[-1], [-0.5], [-1 / 3], [-0.25]]]) + approx.deepEqual(x, [[-1], [-0.5], [-1 / 3], [-0.25]]) }) it('should solve linear system 4 x 4, array and column array', function () { @@ -34,9 +34,9 @@ describe('lusolve', function () { [-1], [-1] ] - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - approx.deepEqual(xs, [[[-1], [-0.5], [-1 / 3], [-0.25]]]) + approx.deepEqual(x, [[-1], [-0.5], [-1 / 3], [-0.25]]) }) it('should solve linear system 4 x 4, matrices', function () { @@ -49,10 +49,10 @@ describe('lusolve', function () { ]) const b = math.matrix([-1, -1, -1, -1]) - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - assert(xs[0] instanceof math.Matrix) - approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) + assert(x instanceof math.Matrix) + approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) }) it('should solve linear system 4 x 4, sparse matrices', function () { @@ -65,10 +65,10 @@ describe('lusolve', function () { ], 'sparse') const b = math.matrix([[-1], [-1], [-1], [-1]], 'sparse') - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - assert(xs[0] instanceof math.Matrix) - approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) + assert(x instanceof math.Matrix) + approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) }) it('should solve linear system 4 x 4, matrix and column matrix', function () { @@ -86,10 +86,10 @@ describe('lusolve', function () { [-1] ]) - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - assert(xs[0] instanceof math.Matrix) - approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) + assert(x instanceof math.Matrix) + approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) }) it('should solve linear system 4 x 4, sparse matrix and column matrix', function () { @@ -107,10 +107,10 @@ describe('lusolve', function () { [-1] ], 'sparse') - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - assert(xs[0] instanceof math.Matrix) - approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) + assert(x instanceof math.Matrix) + approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) }) it('should solve linear system 4 x 4, LUP decomposition (array)', function () { @@ -123,11 +123,11 @@ describe('lusolve', function () { ] const lup = math.lup(m) - const xs = math.lusolve(lup, [-1, -1, -1, -1]) - approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) + const x = math.lusolve(lup, [-1, -1, -1, -1]) + approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) - const ys = math.lusolve(lup, [1, 2, 1, -1]) - approx.deepEqual(ys, [math.matrix([[1], [1], [1 / 3], [-0.25]])]) + const y = math.lusolve(lup, [1, 2, 1, -1]) + approx.deepEqual(y, math.matrix([[1], [1], [1 / 3], [-0.25]])) }) it('should solve linear system 4 x 4, LUP decomposition (matrix)', function () { @@ -140,11 +140,11 @@ describe('lusolve', function () { ]) const lup = math.lup(m) - const xs = math.lusolve(lup, [-1, -1, -1, -1]) - approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) + const x = math.lusolve(lup, [-1, -1, -1, -1]) + approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) - const ys = math.lusolve(lup, [1, 2, 1, -1]) - approx.deepEqual(ys, [math.matrix([[1], [1], [1 / 3], [-0.25]])]) + const y = math.lusolve(lup, [1, 2, 1, -1]) + approx.deepEqual(y, math.matrix([[1], [1], [1 / 3], [-0.25]])) }) it('should solve linear system 4 x 4, LUP decomposition (sparse matrix)', function () { @@ -157,11 +157,11 @@ describe('lusolve', function () { ], 'sparse') const lup = math.lup(m) - const xs = math.lusolve(lup, [-1, -1, -1, -1]) - approx.deepEqual(xs, [math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])]) + const x = math.lusolve(lup, [-1, -1, -1, -1]) + approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1 / 3], [-0.25]])) - const ys = math.lusolve(lup, [1, 2, 1, -1]) - approx.deepEqual(ys, [math.matrix([[1], [1], [1 / 3], [-0.25]])]) + const y = math.lusolve(lup, [1, 2, 1, -1]) + approx.deepEqual(y, math.matrix([[1], [1], [1 / 3], [-0.25]])) }) it('should solve linear system 3 x 3, no permutations, arrays', function () { @@ -173,9 +173,9 @@ describe('lusolve', function () { ] const b = [-2, 4, 2] - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - approx.deepEqual(xs, [[[-5 / 3], [7 / 3], [-1]]]) + approx.deepEqual(x, [[-5 / 3], [7 / 3], [-1]]) }) it('should solve linear system 3 x 3, no permutations, matrix', function () { @@ -187,9 +187,9 @@ describe('lusolve', function () { ]) const b = [-2, 4, 2] - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - approx.deepEqual(xs, [math.matrix([[-5 / 3], [7 / 3], [-1]])]) + approx.deepEqual(x, math.matrix([[-5 / 3], [7 / 3], [-1]])) }) it('should solve linear system 3 x 3, no permutations, sparse matrix', function () { @@ -201,9 +201,9 @@ describe('lusolve', function () { ], 'sparse') const b = [-2, 4, 2] - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - approx.deepEqual(xs, [math.matrix([[-5 / 3], [7 / 3], [-1]])]) + approx.deepEqual(x, math.matrix([[-5 / 3], [7 / 3], [-1]])) }) it('should solve linear system 3 x 3, permutations, arrays', function () { @@ -215,9 +215,9 @@ describe('lusolve', function () { ] const b = [4, -2, 2] - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - approx.deepEqual(xs, [[[-5 / 3], [7 / 3], [-1]]]) + approx.deepEqual(x, [[-5 / 3], [7 / 3], [-1]]) }) it('should solve linear system 4 x 4, permutations, matrix - Issue 437', function () { @@ -231,9 +231,9 @@ describe('lusolve', function () { const b = [0.1, 0.2, 0.15, 0.1] - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - approx.deepEqual(xs, [math.matrix([[0.025], [-0.075], [0], [0.2]])]) + approx.deepEqual(x, math.matrix([[0.025], [-0.075], [0], [0.2]])) }) it('should solve linear system 4 x 4, permutations, sparse - Issue 437', function () { @@ -247,9 +247,9 @@ describe('lusolve', function () { const b = [0.1, 0.2, 0.15, 0.1] - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - approx.deepEqual(xs, [math.matrix([[0.025], [-0.075], [0], [0.2]])]) + approx.deepEqual(x, math.matrix([[0.025], [-0.075], [0], [0.2]])) }) it('should solve linear system 3 x 3, permutations, sparse matrix', function () { @@ -261,9 +261,9 @@ describe('lusolve', function () { ], 'sparse') const b = [4, -2, 2] - const xs = math.lusolve(m, b) + const x = math.lusolve(m, b) - approx.deepEqual(xs, [math.matrix([[-5 / 3], [7 / 3], [-1]])]) + approx.deepEqual(x, math.matrix([[-5 / 3], [7 / 3], [-1]])) }) it('should solve linear system 4 x 4, natural ordering (order=0), partial pivoting, sparse matrix', function () { @@ -277,9 +277,9 @@ describe('lusolve', function () { const b = [1.000000, 1.250000, 1.500000, 1.750000] - const xs = math.lusolve(m, b, 0, 1) + const x = math.lusolve(m, b, 0, 1) - approx.deepEqual(xs, [math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])]) + approx.deepEqual(x, math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])) }) it('should solve linear system 4 x 4, amd(A+A\') (order=1), partial pivoting, sparse matrix', function () { @@ -293,9 +293,9 @@ describe('lusolve', function () { const b = [1.000000, 1.250000, 1.500000, 1.750000] - const xs = math.lusolve(m, b, 1, 1) + const x = math.lusolve(m, b, 1, 1) - approx.deepEqual(xs, [math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])]) + approx.deepEqual(x, math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])) }) it('should solve linear system 4 x 4, amd(A\'*A) (order=2), partial pivoting, sparse matrix', function () { @@ -309,9 +309,9 @@ describe('lusolve', function () { const b = [1.000000, 1.250000, 1.500000, 1.750000] - const xs = math.lusolve(m, b, 2, 1) + const x = math.lusolve(m, b, 2, 1) - approx.deepEqual(xs, [math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])]) + approx.deepEqual(x, math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])) }) it('should solve linear system 4 x 4, amd(A\'*A) (order=3), partial pivoting, sparse matrix', function () { @@ -325,9 +325,9 @@ describe('lusolve', function () { const b = [1.000000, 1.250000, 1.500000, 1.750000] - const xs = math.lusolve(m, b, 3, 1) + const x = math.lusolve(m, b, 3, 1) - approx.deepEqual(xs, [math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])]) + approx.deepEqual(x, math.matrix([[-0.186372], [-0.131621], [0.574586], [2.454950]])) }) it('should throw exception when matrix is singular', function () { diff --git a/test/unit-tests/function/algebra/solver/usolve.test.js b/test/unit-tests/function/algebra/solver/usolve.test.js index 3f579acc80..b105da6d3a 100644 --- a/test/unit-tests/function/algebra/solver/usolve.test.js +++ b/test/unit-tests/function/algebra/solver/usolve.test.js @@ -17,7 +17,7 @@ describe('usolve', function () { const x = math.usolve(m, b) - approx.deepEqual(x, [[[-1], [-1], [-1], [4]]]) + approx.deepEqual(x, [[-1], [-1], [-1], [4]]) }) it('should solve linear system 4 x 4, array and column array', function () { @@ -36,7 +36,7 @@ describe('usolve', function () { ] const x = math.usolve(m, b) - approx.deepEqual(x, [[[-1], [-1], [-1], [4]]]) + approx.deepEqual(x, [[-1], [-1], [-1], [4]]) }) it('should solve linear system 4 x 4, matrices', function () { @@ -51,8 +51,8 @@ describe('usolve', function () { const x = math.usolve(m, b) - assert(x[0] instanceof math.Matrix) - approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) + assert(x instanceof math.Matrix) + approx.deepEqual(x, math.matrix([[-1], [-1], [-1], [4]])) }) it('should solve linear system 4 x 4, sparse matrices', function () { @@ -67,8 +67,8 @@ describe('usolve', function () { const x = math.usolve(m, b) - assert(x[0] instanceof math.Matrix) - approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) + assert(x instanceof math.Matrix) + approx.deepEqual(x, math.matrix([[-1], [-1], [-1], [4]])) }) it('should solve linear system 4 x 4, matrix and column matrix', function () { @@ -88,8 +88,8 @@ describe('usolve', function () { const x = math.usolve(m, b) - assert(x[0] instanceof math.Matrix) - approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) + assert(x instanceof math.Matrix) + approx.deepEqual(x, math.matrix([[-1], [-1], [-1], [4]])) }) it('should solve linear system 4 x 4, sparse matrix and column matrix', function () { @@ -109,8 +109,8 @@ describe('usolve', function () { const x = math.usolve(m, b) - assert(x[0] instanceof math.Matrix) - approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) + assert(x instanceof math.Matrix) + approx.deepEqual(x, math.matrix([[-1], [-1], [-1], [4]])) }) it('should throw exception when matrix is singular', function () { @@ -118,21 +118,4 @@ describe('usolve', function () { assert.throws(function () { math.usolve(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) assert.throws(function () { math.usolve(math.matrix([[1, 1], [0, 0]], 'sparse'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) }) - - it('should solve systems with singular matrices', function () { - approx.deepEqual( - math.usolve([[1, 1, -1], [0, 0, 1], [0, 0, 2]], [1, 2, 4]), - [[[3], [0], [2]], [[2], [1], [2]]] - ) - - approx.deepEqual( - math.usolve([[0, 1, 2], [0, 1, 1], [0, 0, 0]], [2, 2, 0]), - [[[0], [2], [0]], [[1], [2], [0]]] - ) - - approx.deepEqual( - math.usolve([[0, 1, 1], [0, 1, 1], [0, 0, 0]], [2, 2, 0]), - [[[0], [2], [0]], [[0], [1], [1]], [[1], [2], [0]]] - ) - }) }) diff --git a/test/unit-tests/function/algebra/solver/usolveAll.test.js b/test/unit-tests/function/algebra/solver/usolveAll.test.js new file mode 100644 index 0000000000..a80a58e0d5 --- /dev/null +++ b/test/unit-tests/function/algebra/solver/usolveAll.test.js @@ -0,0 +1,138 @@ +// test usolveAll +import assert from 'assert' + +import approx from '../../../../../tools/approx' +import math from '../../../../../src/bundleAny' + +describe('usolveAll', function () { + it('should solve linear system 4 x 4, arrays', function () { + const m = + [ + [1, 1, 1, 1], + [0, 1, 1, 1], + [0, 0, 1, 1], + [0, 0, 0, 1] + ] + const b = [1, 2, 3, 4] + + const x = math.usolveAll(m, b) + + approx.deepEqual(x, [[[-1], [-1], [-1], [4]]]) + }) + + it('should solve linear system 4 x 4, array and column array', function () { + const m = + [ + [1, 1, 1, 1], + [0, 1, 1, 1], + [0, 0, 1, 1], + [0, 0, 0, 1] + ] + const b = [ + [1], + [2], + [3], + [4] + ] + const x = math.usolveAll(m, b) + + approx.deepEqual(x, [[[-1], [-1], [-1], [4]]]) + }) + + it('should solve linear system 4 x 4, matrices', function () { + const m = math.matrix( + [ + [1, 1, 1, 1], + [0, 1, 1, 1], + [0, 0, 1, 1], + [0, 0, 0, 1] + ]) + const b = math.matrix([1, 2, 3, 4]) + + const x = math.usolveAll(m, b) + + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) + }) + + it('should solve linear system 4 x 4, sparse matrices', function () { + const m = math.sparse( + [ + [1, 1, 1, 1], + [0, 1, 1, 1], + [0, 0, 1, 1], + [0, 0, 0, 1] + ]) + const b = math.matrix([[1], [2], [3], [4]], 'sparse') + + const x = math.usolveAll(m, b) + + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) + }) + + it('should solve linear system 4 x 4, matrix and column matrix', function () { + const m = math.matrix( + [ + [1, 1, 1, 1], + [0, 1, 1, 1], + [0, 0, 1, 1], + [0, 0, 0, 1] + ]) + const b = math.matrix([ + [1], + [2], + [3], + [4] + ]) + + const x = math.usolveAll(m, b) + + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) + }) + + it('should solve linear system 4 x 4, sparse matrix and column matrix', function () { + const m = math.matrix( + [ + [1, 1, 1, 1], + [0, 1, 1, 1], + [0, 0, 1, 1], + [0, 0, 0, 1] + ], 'sparse') + const b = math.matrix([ + [1], + [2], + [3], + [4] + ], 'sparse') + + const x = math.usolveAll(m, b) + + assert(x[0] instanceof math.Matrix) + approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) + }) + + it('should throw exception when matrix is singular', function () { + assert.throws(function () { math.usolveAll([[1, 1], [0, 0]], [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) + assert.throws(function () { math.usolveAll(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) + assert.throws(function () { math.usolveAll(math.matrix([[1, 1], [0, 0]], 'sparse'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) + }) + + it('should solve systems with singular matrices', function () { + approx.deepEqual( + math.usolveAll([[1, 1, -1], [0, 0, 1], [0, 0, 2]], [1, 2, 4]), + [[[3], [0], [2]], [[2], [1], [2]]] + ) + + approx.deepEqual( + math.usolveAll([[0, 1, 2], [0, 1, 1], [0, 0, 0]], [2, 2, 0]), + [[[0], [2], [0]], [[1], [2], [0]]] + ) + + approx.deepEqual( + math.usolveAll([[0, 1, 1], [0, 1, 1], [0, 0, 0]], [2, 2, 0]), + [[[0], [2], [0]], [[0], [1], [1]], [[1], [2], [0]]] + ) + }) +}) From a42b1104017f226b75397bca23a2c5484a8469a0 Mon Sep 17 00:00:00 2001 From: Michal Grno Date: Tue, 25 Aug 2020 20:42:04 +0200 Subject: [PATCH 8/9] made *solveAll return [] for non-solvable, implemented sparse algorithms --- src/function/algebra/solver/lsolveAll.js | 128 ++++++++++-------- src/function/algebra/solver/usolveAll.js | 116 ++++++++-------- .../function/algebra/solver/lsolveAll.test.js | 27 +++- .../function/algebra/solver/usolveAll.test.js | 27 +++- 4 files changed, 174 insertions(+), 124 deletions(-) diff --git a/src/function/algebra/solver/lsolveAll.js b/src/function/algebra/solver/lsolveAll.js index 9be4a7231b..8c9e395601 100644 --- a/src/function/algebra/solver/lsolveAll.js +++ b/src/function/algebra/solver/lsolveAll.js @@ -89,7 +89,7 @@ export const createLsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ ty if (k === 0) { // There is no valid solution - throw new Error('Linear system cannot be solved since matrix is singular') + return [] } else { // This RHS is invalid but other solutions may still exist B.splice(k, 1) @@ -114,74 +114,84 @@ export const createLsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ ty return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] })) } - function _sparseForwardSubstitution (m, b) { - // validate matrix and vector, return copy of column vector b - b = solveValidation(m, b, true) - // column vector data - const bdata = b._data - // rows & columns + function _sparseForwardSubstitution (m, b_) { + // array of right-hand sides + const B = [solveValidation(m, b_, true)._data.map(e => e[0])] + const rows = m._size[0] const columns = m._size[1] - // matrix arrays + const values = m._values const index = m._index const ptr = m._ptr - // vars - let i, k - // result - const x = [] - // forward solve m * x = b, loop columns - for (let j = 0; j < columns; j++) { - // b[j] - const bj = bdata[j][0] || 0 - // forward substitution (outer product) avoids inner looping when bj === 0 - if (!equalScalar(bj, 0)) { - // value @ [j, j] - let vjj = 0 - // lower triangular matrix values & index (column j) - const jvalues = [] - const jindex = [] - // last index in column - let l = ptr[j + 1] - // values in column, find value @ [j, j] - for (k = ptr[j]; k < l; k++) { - // row - i = index[k] - // check row (rows are not sorted!) - if (i === j) { - // update vjj - vjj = values[k] - } else if (i > j) { + + // loop columns + for (let i = 0; i < columns; i++) { + let L = B.length + + // loop right-hand sides + for (let k = 0; k < L; k++) { + const b = B[k] + + // values & indices (column i) + const iValues = [] + const iIndices = [] + + // first & last indeces in column + const firstIndex = ptr[i] + const lastIndex = ptr[i + 1] + + // find the value at [i, i] + let Mii = 0 + for (let j = firstIndex; j < lastIndex; j++) { + const J = index[j] + // check row + if (J === i) { + Mii = values[j] + } else if (J > i) { // store lower triangular - jvalues.push(values[k]) - jindex.push(i) + iValues.push(values[j]) + iIndices.push(J) } } - // at this point we must have a value @ [j, j] - if (equalScalar(vjj, 0)) { - // system cannot be solved, there is no value @ [j, j] - throw new Error('Linear system cannot be solved since matrix is singular') - } - // calculate xj - const xj = divideScalar(bj, vjj) - // loop lower triangular - for (k = 0, l = jindex.length; k < l; k++) { - // row - i = jindex[k] - // update copy of b - bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, jvalues[k]))] + + if (!equalScalar(Mii, 0)) { + // non-singular row + + b[i] = divideScalar(b[i], Mii) + + for (let j = 0, lastIndex = iIndices.length; j < lastIndex; j++) { + const J = iIndices[j] + b[J] = subtract(b[J], multiplyScalar(b[i], iValues[j])) + } + } else if (!equalScalar(b[i], 0)) { + // singular row, nonzero RHS + + if (k === 0) { + // There is no valid solution + return [] + } else { + // This RHS is invalid but other solutions may still exist + B.splice(k, 1) + k -= 1 + L -= 1 + } + } else if (k === 0) { + // singular row, RHS is zero + + const bNew = [...b] + bNew[i] = 1 + + for (let j = 0, lastIndex = iIndices.length; j < lastIndex; j++) { + const J = iIndices[j] + bNew[J] = subtract(bNew[J], iValues[j]) + } + + B.push(bNew) } - // update x - x[j] = [xj] - } else { - // update x - x[j] = [0] } } - // return vector - return [new DenseMatrix({ - data: x, - size: [rows, 1] - })] + + return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] })) } }) diff --git a/src/function/algebra/solver/usolveAll.js b/src/function/algebra/solver/usolveAll.js index c07fafe723..0b7e56d03b 100644 --- a/src/function/algebra/solver/usolveAll.js +++ b/src/function/algebra/solver/usolveAll.js @@ -89,7 +89,7 @@ export const createUsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ ty if (k === 0) { // There is no valid solution - throw new Error('Linear system cannot be solved since matrix is singular') + return [] } else { // This RHS is invalid but other solutions may still exist B.splice(k, 1) @@ -114,11 +114,9 @@ export const createUsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ ty return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] })) } - function _sparseBackwardSubstitution (m, b) { - // make b into a column vector - b = solveValidation(m, b, true) - - const bdata = b._data + function _sparseBackwardSubstitution (m, b_) { + // array of right-hand sides + const B = [solveValidation(m, b_, true)._data.map(e => e[0])] const rows = m._size[0] const columns = m._size[1] @@ -127,67 +125,75 @@ export const createUsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ ty const index = m._index const ptr = m._ptr - let i, k - - // result - const x = [] - // loop columns backwards - for (let j = columns - 1; j >= 0; j--) { - // b[j] - const bj = bdata[j][0] || 0 + for (let i = columns - 1; i >= 0; i--) { + let L = B.length - if (!equalScalar(bj, 0)) { - // value at [j, j] - let vjj = 0 + // loop right-hand sides + for (let k = 0; k < L; k++) { + const b = B[k] - // upper triangular matrix values & index (column j) - const jValues = [] - const jIndices = [] + // values & indices (column i) + const iValues = [] + const iIndices = [] // first & last indeces in column - const firstIndex = ptr[j] - let lastIndex = ptr[j + 1] + const firstIndex = ptr[i] + const lastIndex = ptr[i + 1] - // values in column, find value at [j, j], loop backwards - for (k = lastIndex - 1; k >= firstIndex; k--) { - // row - i = index[k] + // find the value at [i, i] + let Mii = 0 + for (let j = lastIndex - 1; j >= firstIndex; j--) { + const J = index[j] // check row - if (i === j) { - // update vjj - vjj = values[k] - } else if (i < j) { + if (J === i) { + Mii = values[j] + } else if (J < i) { // store upper triangular - jValues.push(values[k]) - jIndices.push(i) + iValues.push(values[j]) + iIndices.push(J) } } - // at this point we must have a value at [j, j] - if (equalScalar(vjj, 0)) { - // system cannot be solved, there is no value at [j, j] - throw new Error('Linear system cannot be solved since matrix is singular') - } - // calculate xj - const xj = divideScalar(bj, vjj) - // loop upper triangular - for (k = 0, lastIndex = jIndices.length; k < lastIndex; k++) { - // row - i = jIndices[k] - // update copy of b - bdata[i] = [subtract(bdata[i][0], multiplyScalar(xj, jValues[k]))] + + if (!equalScalar(Mii, 0)) { + // non-singular row + + b[i] = divideScalar(b[i], Mii) + + // loop upper triangular + for (let j = 0, lastIndex = iIndices.length; j < lastIndex; j++) { + const J = iIndices[j] + b[J] = subtract(b[J], multiplyScalar(b[i], iValues[j])) + } + } else if (!equalScalar(b[i], 0)) { + // singular row, nonzero RHS + + if (k === 0) { + // There is no valid solution + return [] + } else { + // This RHS is invalid but other solutions may still exist + B.splice(k, 1) + k -= 1 + L -= 1 + } + } else if (k === 0) { + // singular row, RHS is zero + + const bNew = [...b] + bNew[i] = 1 + + // loop upper triangular + for (let j = 0, lastIndex = iIndices.length; j < lastIndex; j++) { + const J = iIndices[j] + bNew[J] = subtract(bNew[J], iValues[j]) + } + + B.push(bNew) } - // update x - x[j] = [xj] - } else { - // update x - x[j] = [0] } } - // return vector - return [new DenseMatrix({ - data: x, - size: [rows, 1] - })] + + return B.map(x => new DenseMatrix({ data: x.map(e => [e]), size: [rows, 1] })) } }) diff --git a/test/unit-tests/function/algebra/solver/lsolveAll.test.js b/test/unit-tests/function/algebra/solver/lsolveAll.test.js index 0415a037d6..6df922ba84 100644 --- a/test/unit-tests/function/algebra/solver/lsolveAll.test.js +++ b/test/unit-tests/function/algebra/solver/lsolveAll.test.js @@ -113,13 +113,13 @@ describe('lsolveAll', function () { approx.deepEqual(x, [math.matrix([[1], [1], [1], [1]])]) }) - it('should throw exception when matrix is singular', function () { - assert.throws(function () { math.lsolveAll([[1, 1], [0, 0]], [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) - assert.throws(function () { math.lsolveAll(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) - assert.throws(function () { math.lsolveAll(math.matrix([[1, 1], [0, 0]], 'sparse'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) + it('should return an empty array when there is no solution', function () { + assert.deepStrictEqual([], math.lsolveAll([[1, 1], [0, 0]], [1, 1])) + assert.deepStrictEqual([], math.lsolveAll(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1])) + assert.deepStrictEqual([], math.lsolveAll(math.matrix([[1, 1], [0, 0]], 'sparse'), [1, 1])) }) - it('should solve systems with singular matrices', function () { + it('should solve systems with singular dense matrices', function () { approx.deepEqual( math.lsolveAll([[2, 0, 0], [1, 0, 0], [-1, 1, 1]], [4, 2, 1]), [[[2], [0], [3]], [[2], [1], [2]]] @@ -135,4 +135,21 @@ describe('lsolveAll', function () { [[[0], [2], [0]], [[1], [1], [0]], [[0], [2], [1]]] ) }) + + it('should solve systems with singular sparse matrices', function () { + approx.deepEqual( + math.lsolveAll(math.matrix([[2, 0, 0], [1, 0, 0], [-1, 1, 1]], 'sparse'), [4, 2, 1]), + [math.matrix([[2], [0], [3]]), math.matrix([[2], [1], [2]])] + ) + + approx.deepEqual( + math.lsolveAll(math.matrix([[0, 0, 0], [1, 1, 0], [2, 1, 0]], 'sparse'), [0, 2, 2]), + [math.matrix([[0], [2], [0]]), math.matrix([[0], [2], [1]])] + ) + + approx.deepEqual( + math.lsolveAll(math.matrix([[0, 0, 0], [1, 1, 0], [1, 1, 0]], 'sparse'), [0, 2, 2]), + [math.matrix([[0], [2], [0]]), math.matrix([[1], [1], [0]]), math.matrix([[0], [2], [1]])] + ) + }) }) diff --git a/test/unit-tests/function/algebra/solver/usolveAll.test.js b/test/unit-tests/function/algebra/solver/usolveAll.test.js index a80a58e0d5..f7e3453146 100644 --- a/test/unit-tests/function/algebra/solver/usolveAll.test.js +++ b/test/unit-tests/function/algebra/solver/usolveAll.test.js @@ -113,13 +113,13 @@ describe('usolveAll', function () { approx.deepEqual(x, [math.matrix([[-1], [-1], [-1], [4]])]) }) - it('should throw exception when matrix is singular', function () { - assert.throws(function () { math.usolveAll([[1, 1], [0, 0]], [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) - assert.throws(function () { math.usolveAll(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) - assert.throws(function () { math.usolveAll(math.matrix([[1, 1], [0, 0]], 'sparse'), [1, 1]) }, /Error: Linear system cannot be solved since matrix is singular/) + it('should return an empty array when there is no solution', function () { + assert.deepStrictEqual([], math.usolveAll([[1, 1], [0, 0]], [1, 1])) + assert.deepStrictEqual([], math.usolveAll(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1])) + assert.deepStrictEqual([], math.usolveAll(math.matrix([[1, 1], [0, 0]], 'sparse'), [1, 1])) }) - it('should solve systems with singular matrices', function () { + it('should solve systems with singular dense matrices', function () { approx.deepEqual( math.usolveAll([[1, 1, -1], [0, 0, 1], [0, 0, 2]], [1, 2, 4]), [[[3], [0], [2]], [[2], [1], [2]]] @@ -135,4 +135,21 @@ describe('usolveAll', function () { [[[0], [2], [0]], [[0], [1], [1]], [[1], [2], [0]]] ) }) + + it('should solve systems with singular sparse matrices', function () { + approx.deepEqual( + math.usolveAll(math.matrix([[1, 1, -1], [0, 0, 1], [0, 0, 2]], 'sparse'), [1, 2, 4]), + [math.matrix([[3], [0], [2]]), math.matrix([[2], [1], [2]])] + ) + + approx.deepEqual( + math.usolveAll(math.matrix([[0, 1, 2], [0, 1, 1], [0, 0, 0]], 'sparse'), [2, 2, 0]), + [math.matrix([[0], [2], [0]]), math.matrix([[1], [2], [0]])] + ) + + approx.deepEqual( + math.usolveAll(math.matrix([[0, 1, 1], [0, 1, 1], [0, 0, 0]], 'sparse'), [2, 2, 0]), + [math.matrix([[0], [2], [0]]), math.matrix([[0], [1], [1]]), math.matrix([[1], [2], [0]])] + ) + }) }) From c90cbbd07417f34790194c881c8c89a1fc2ea58b Mon Sep 17 00:00:00 2001 From: Michal Grno Date: Tue, 25 Aug 2020 21:01:05 +0200 Subject: [PATCH 9/9] improved documentation for *solve(All) --- src/expression/embeddedDocs/embeddedDocs.js | 4 ++++ .../embeddedDocs/function/algebra/lsolve.js | 4 ++-- .../embeddedDocs/function/algebra/lsolveAll.js | 17 +++++++++++++++++ .../embeddedDocs/function/algebra/usolve.js | 4 ++-- .../embeddedDocs/function/algebra/usolveAll.js | 15 +++++++++++++++ src/function/algebra/solver/lsolve.js | 4 ++-- src/function/algebra/solver/lsolveAll.js | 4 ++-- src/function/algebra/solver/usolve.js | 4 ++-- src/function/algebra/solver/usolveAll.js | 4 ++-- 9 files changed, 48 insertions(+), 12 deletions(-) create mode 100644 src/expression/embeddedDocs/function/algebra/lsolveAll.js create mode 100644 src/expression/embeddedDocs/function/algebra/usolveAll.js diff --git a/src/expression/embeddedDocs/embeddedDocs.js b/src/expression/embeddedDocs/embeddedDocs.js index 5c8bfb2267..840555c970 100644 --- a/src/expression/embeddedDocs/embeddedDocs.js +++ b/src/expression/embeddedDocs/embeddedDocs.js @@ -167,11 +167,13 @@ import { addDocs } from './function/arithmetic/add' import { absDocs } from './function/arithmetic/abs' import { qrDocs } from './function/algebra/qr' import { usolveDocs } from './function/algebra/usolve' +import { usolveAllDocs } from './function/algebra/usolveAll' import { sluDocs } from './function/algebra/slu' import { rationalizeDocs } from './function/algebra/rationalize' import { simplifyDocs } from './function/algebra/simplify' import { lupDocs } from './function/algebra/lup' import { lsolveDocs } from './function/algebra/lsolve' +import { lsolveAllDocs } from './function/algebra/lsolveAll' import { derivativeDocs } from './function/algebra/derivative' import { versionDocs } from './constants/version' import { trueDocs } from './constants/true' @@ -310,12 +312,14 @@ export const embeddedDocs = { // functions - algebra derivative: derivativeDocs, lsolve: lsolveDocs, + lsolveAll: lsolveAllDocs, lup: lupDocs, lusolve: lusolveDocs, simplify: simplifyDocs, rationalize: rationalizeDocs, slu: sluDocs, usolve: usolveDocs, + usolveAll: usolveAllDocs, qr: qrDocs, // functions - arithmetic diff --git a/src/expression/embeddedDocs/function/algebra/lsolve.js b/src/expression/embeddedDocs/function/algebra/lsolve.js index 0cf27bc9d6..7fc0474e6b 100644 --- a/src/expression/embeddedDocs/function/algebra/lsolve.js +++ b/src/expression/embeddedDocs/function/algebra/lsolve.js @@ -5,13 +5,13 @@ export const lsolveDocs = { 'x=lsolve(L, b)' ], description: - 'Solves the linear system L * x = b where L is an [n x n] lower triangular matrix and b is a [n] column vector.', + 'Finds one solution of the linear system L * x = b where L is an [n x n] lower triangular matrix and b is a [n] column vector.', examples: [ 'a = [-2, 3; 2, 1]', 'b = [11, 9]', 'x = lsolve(a, b)' ], seealso: [ - 'lup', 'lusolve', 'usolve', 'matrix', 'sparse' + 'lsolveAll', 'lup', 'lusolve', 'usolve', 'matrix', 'sparse' ] } diff --git a/src/expression/embeddedDocs/function/algebra/lsolveAll.js b/src/expression/embeddedDocs/function/algebra/lsolveAll.js new file mode 100644 index 0000000000..e5ff44cf24 --- /dev/null +++ b/src/expression/embeddedDocs/function/algebra/lsolveAll.js @@ -0,0 +1,17 @@ +export const lsolveAllDocs = { + name: 'lsolveAll', + category: 'Algebra', + syntax: [ + 'x=lsolveAll(L, b)' + ], + description: + 'Finds all solutions of the linear system L * x = b where L is an [n x n] lower triangular matrix and b is a [n] column vector.', + examples: [ + 'a = [-2, 3; 2, 1]', + 'b = [11, 9]', + 'x = lsolve(a, b)' + ], + seealso: [ + 'lsolve', 'lup', 'lusolve', 'usolve', 'matrix', 'sparse' + ] +} diff --git a/src/expression/embeddedDocs/function/algebra/usolve.js b/src/expression/embeddedDocs/function/algebra/usolve.js index 7c0903a1de..9a22974a37 100644 --- a/src/expression/embeddedDocs/function/algebra/usolve.js +++ b/src/expression/embeddedDocs/function/algebra/usolve.js @@ -5,11 +5,11 @@ export const usolveDocs = { 'x=usolve(U, b)' ], description: - 'Solves the linear system U * x = b where U is an [n x n] upper triangular matrix and b is a [n] column vector.', + 'Finds one solution of the linear system U * x = b where U is an [n x n] upper triangular matrix and b is a [n] column vector.', examples: [ 'x=usolve(sparse([1, 1, 1, 1; 0, 1, 1, 1; 0, 0, 1, 1; 0, 0, 0, 1]), [1; 2; 3; 4])' ], seealso: [ - 'lup', 'lusolve', 'lsolve', 'matrix', 'sparse' + 'usolveAll', 'lup', 'lusolve', 'lsolve', 'matrix', 'sparse' ] } diff --git a/src/expression/embeddedDocs/function/algebra/usolveAll.js b/src/expression/embeddedDocs/function/algebra/usolveAll.js new file mode 100644 index 0000000000..c2131c1be1 --- /dev/null +++ b/src/expression/embeddedDocs/function/algebra/usolveAll.js @@ -0,0 +1,15 @@ +export const usolveAllDocs = { + name: 'usolveAll', + category: 'Algebra', + syntax: [ + 'x=usolve(U, b)' + ], + description: + 'Finds all solutions of the linear system U * x = b where U is an [n x n] upper triangular matrix and b is a [n] column vector.', + examples: [ + 'x=usolve(sparse([1, 1, 1, 1; 0, 1, 1, 1; 0, 0, 1, 1; 0, 0, 0, 1]), [1; 2; 3; 4])' + ], + seealso: [ + 'usolve', 'lup', 'lusolve', 'lsolve', 'matrix', 'sparse' + ] +} diff --git a/src/function/algebra/solver/lsolve.js b/src/function/algebra/solver/lsolve.js index 64ba9bc331..eb1623930f 100644 --- a/src/function/algebra/solver/lsolve.js +++ b/src/function/algebra/solver/lsolve.js @@ -16,7 +16,7 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed const solveValidation = createSolveValidation({ DenseMatrix }) /** - * Solves the linear equation system by forwards substitution. Matrix must be a lower triangular matrix. + * Finds one solution of a linear equation system by forwards substitution. Matrix must be a lower triangular matrix. Throws an error if there's no solution. * * `L * x = b` * @@ -32,7 +32,7 @@ export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed * * See also: * - * lup, slu, usolve, lusolve + * lsolveAll, lup, slu, usolve, lusolve * * @param {Matrix, Array} L A N x N matrix or array (L) * @param {Matrix, Array} b A column vector with the b values diff --git a/src/function/algebra/solver/lsolveAll.js b/src/function/algebra/solver/lsolveAll.js index 8c9e395601..7a4c825785 100644 --- a/src/function/algebra/solver/lsolveAll.js +++ b/src/function/algebra/solver/lsolveAll.js @@ -16,7 +16,7 @@ export const createLsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ ty const solveValidation = createSolveValidation({ DenseMatrix }) /** - * Solves the linear equation system by forwards substitution. Matrix must be a lower triangular matrix. + * Finds all solutions of a linear equation system by forwards substitution. Matrix must be a lower triangular matrix. * * `L * x = b` * @@ -32,7 +32,7 @@ export const createLsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ ty * * See also: * - * lup, slu, usolve, lusolve + * lsolve, lup, slu, usolve, lusolve * * @param {Matrix, Array} L A N x N matrix or array (L) * @param {Matrix, Array} b A column vector with the b values diff --git a/src/function/algebra/solver/usolve.js b/src/function/algebra/solver/usolve.js index aea9e6cc0d..87bc9c93b2 100644 --- a/src/function/algebra/solver/usolve.js +++ b/src/function/algebra/solver/usolve.js @@ -16,7 +16,7 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed const solveValidation = createSolveValidation({ DenseMatrix }) /** - * Solves the linear equation system by backward substitution. Matrix must be an upper triangular matrix. + * Finds one solution of a linear equation system by backward substitution. Matrix must be an upper triangular matrix. Throws an error if there's no solution. * * `U * x = b` * @@ -32,7 +32,7 @@ export const createUsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed * * See also: * - * lup, slu, usolve, lusolve + * usolveAll, lup, slu, usolve, lusolve * * @param {Matrix, Array} U A N x N matrix or array (U) * @param {Matrix, Array} b A column vector with the b values diff --git a/src/function/algebra/solver/usolveAll.js b/src/function/algebra/solver/usolveAll.js index 0b7e56d03b..9e1b9822fd 100644 --- a/src/function/algebra/solver/usolveAll.js +++ b/src/function/algebra/solver/usolveAll.js @@ -16,7 +16,7 @@ export const createUsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ ty const solveValidation = createSolveValidation({ DenseMatrix }) /** - * Solves the linear equation system by backward substitution. Matrix must be an upper triangular matrix. + * Finds all solutions of a linear equation system by backward substitution. Matrix must be an upper triangular matrix. * * `U * x = b` * @@ -32,7 +32,7 @@ export const createUsolveAll = /* #__PURE__ */ factory(name, dependencies, ({ ty * * See also: * - * lup, slu, usolve, lusolve + * usolve, lup, slu, usolve, lusolve * * @param {Matrix, Array} U A N x N matrix or array (U) * @param {Matrix, Array} b A column vector with the b values