From ab4e364138484349cf91178c4ce5d474ef57569c Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Wed, 4 Oct 2023 13:39:05 -0700 Subject: [PATCH] feat: Add option to eigs() to turn off eigenvector computation For large matrices, the eigenvector computation can be noticeably expensive and so it's worthwhile to have a way to turn it off if the eigenvectors will not be used. Resolves #2180. --- .../embeddedDocs/function/matrix/eigs.js | 5 +- src/function/matrix/eigs.js | 73 ++++++++++++------- src/function/matrix/eigs/realSymmetric.js | 65 ++++++++++------- test/unit-tests/function/matrix/eigs.test.js | 16 +++- 4 files changed, 102 insertions(+), 57 deletions(-) diff --git a/src/expression/embeddedDocs/function/matrix/eigs.js b/src/expression/embeddedDocs/function/matrix/eigs.js index ec389a198f..09171abf4c 100644 --- a/src/expression/embeddedDocs/function/matrix/eigs.js +++ b/src/expression/embeddedDocs/function/matrix/eigs.js @@ -4,9 +4,10 @@ export const eigsDocs = { syntax: [ 'eigs(x)' ], - description: 'Calculate the eigenvalues and eigenvectors of a real symmetric matrix', + description: 'Calculate the eigenvalues and optionally eigenvectors of a square matrix', examples: [ - 'eigs([[5, 2.3], [2.3, 1]])' + 'eigs([[5, 2.3], [2.3, 1]])', + 'eigs([[1, 2, 3], [4, 5, 6], [7, 8, 9]], { precision: 1e-6, eigenvectors: false }' ], seealso: [ 'inv' diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js index 50edd5cb70..140176126f 100644 --- a/src/function/matrix/eigs.js +++ b/src/function/matrix/eigs.js @@ -13,7 +13,7 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, const doComplexEigs = createComplexEigs({ config, addScalar, subtract, multiply, multiplyScalar, flatten, divideScalar, sqrt, abs, bignumber, diag, size, reshape, qr, inv, usolve, usolveAll, equal, complex, larger, smaller, matrixFromColumns, dot }) /** - * Compute eigenvalues and eigenvectors of a square matrix. + * Compute eigenvalues and optionally eigenvectors of a square matrix. * The eigenvalues are sorted by their absolute value, ascending, and * returned as a vector in the `values` property of the returned project. * An eigenvalue with algebraic multiplicity k will be listed k times, so @@ -31,9 +31,21 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, * in that case, however, you may still find useful information * in `err.values` and `err.vectors`. * + * Note that the 'precision' option does not directly specify the _accuracy_ + * of the returned eigenvalues. Rather, it determines how small an entry + * of the iterative approximations to an upper triangular matrix must be + * in order to be considered zero. The actual accuracy of the returned + * eigenvalues may be greater or less than the precision, depending on the + * conditioning of the matrix and how far apart or close the actual + * eigenvalues are. Note that currently, relatively simple, "traditional" + * methods of eigenvalue computation are being used; this is not a modern, + * high-precision eigenvalue computation. That said, it should typically + * produce fairly reasonable results. + * * Syntax: * * math.eigs(x, [prec]) + * math.eigs(x, {options}) * * Examples: * @@ -47,14 +59,17 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, * const UTxHxU = multiply(transpose(U), H, U) // diagonalizes H if possible * E[0] == UTxHxU[0][0] // returns true always * + * // Compute only approximate eigenvalues: + * const {values} = eigs(H, {eigenvectors: false, precision: 1e-6}) + * * See also: * * inv * * @param {Array | Matrix} x Matrix to be diagonalized * - * @param {number | BigNumber} [prec] Precision, default value: 1e-15 - * @return {{values: Array|Matrix, eigenvectors: Array}} Object containing an array of eigenvalues and an array of {value: number|BigNumber, vector: Array|Matrix} objects. + * @param {number | BigNumber | OptsObject} [opts] Object with keys `precision`, defaulting to config.epsilon, and `eigenvectors`, defaulting to true and specifying whether to compute eigenvectors. If just a number, specifies precision. + * @return {{values: Array|Matrix, eigenvectors?: Array}} Object containing an array of eigenvalues and an array of {value: number|BigNumber, vector: Array|Matrix} objects. The eigenvectors property is undefined if eigenvectors were not requested. * */ return typed('eigs', { @@ -67,38 +82,46 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, // is a roundabout way of doing type detection. Array: function (x) { return doEigs(matrix(x)) }, 'Array, number|BigNumber': function (x, prec) { - return doEigs(matrix(x), prec) + return doEigs(matrix(x), { precision: prec }) }, + 'Array, Object' (x, opts) { return doEigs(matrix(x), opts) }, Matrix: function (mat) { - return doEigs(mat, undefined, true) + return doEigs(mat, { matricize: true }) }, 'Matrix, number|BigNumber': function (mat, prec) { - return doEigs(mat, prec, true) + return doEigs(mat, { precision: prec, matricize: true }) + }, + 'Matrix, Object': function (mat, opts) { + const useOpts = { matricize: true } + Object.assign(useOpts, opts) + return doEigs(mat, useOpts) } }) - function doEigs (mat, prec, matricize = false) { - const result = computeValuesAndVectors(mat, prec) - if (matricize) { + function doEigs (mat, opts = {}) { + const computeVectors = 'eigenvectors' in opts ? opts.eigenvectors : true + const prec = opts.precision ?? config.epsilon + const result = computeValuesAndVectors(mat, prec, computeVectors) + if (opts.matricize) { result.values = matrix(result.values) - result.eigenvectors = result.eigenvectors.map(({ value, vector }) => - ({ value, vector: matrix(vector) })) - } - Object.defineProperty(result, 'vectors', { - enumerable: false, // to make sure that the eigenvectors can still be - // converted to string. - get: () => { - throw new Error('eigs(M).vectors replaced with eigs(M).eigenvectors') + if (computeVectors) { + result.eigenvectors = result.eigenvectors.map(({ value, vector }) => + ({ value, vector: matrix(vector) })) } - }) + } + if (computeVectors) { + Object.defineProperty(result, 'vectors', { + enumerable: false, // to make sure that the eigenvectors can still be + // converted to string. + get: () => { + throw new Error('eigs(M).vectors replaced with eigs(M).eigenvectors') + } + }) + } return result } - function computeValuesAndVectors (mat, prec) { - if (prec === undefined) { - prec = config.epsilon - } - + function computeValuesAndVectors (mat, prec, computeVectors) { const arr = mat.toArray() // NOTE: arr is guaranteed to be unaliased // and so safe to modify in place const asize = mat.size() @@ -114,12 +137,12 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, if (isSymmetric(arr, N, prec)) { const type = coerceTypes(mat, arr, N) // modifies arr by side effect - return doRealSymmetric(arr, N, prec, type) + return doRealSymmetric(arr, N, prec, type, computeVectors) } } const type = coerceTypes(mat, arr, N) // modifies arr by side effect - return doComplexEigs(arr, N, prec, type) + return doComplexEigs(arr, N, prec, type, computeVectors) } /** @return {boolean} */ diff --git a/src/function/matrix/eigs/realSymmetric.js b/src/function/matrix/eigs/realSymmetric.js index 7e24d406da..2090a9a343 100644 --- a/src/function/matrix/eigs/realSymmetric.js +++ b/src/function/matrix/eigs/realSymmetric.js @@ -7,28 +7,31 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c * @param {number} prec * @param {'number' | 'BigNumber'} type */ - function main (arr, N, prec = config.epsilon, type) { + function main (arr, N, prec = config.epsilon, type, computeVectors) { if (type === 'number') { - return diag(arr, prec) + return diag(arr, prec, computeVectors) } if (type === 'BigNumber') { - return diagBig(arr, prec) + return diagBig(arr, prec, computeVectors) } throw TypeError('Unsupported data type: ' + type) } // diagonalization implementation for number (efficient) - function diag (x, precision) { + function diag (x, precision, computeVectors) { const N = x.length const e0 = Math.abs(precision / N) let psi - let Sij = new Array(N) - // Sij is Identity Matrix - for (let i = 0; i < N; i++) { - Sij[i] = Array(N).fill(0) - Sij[i][i] = 1.0 + let Sij + if (computeVectors) { + Sij = new Array(N) + // Sij is Identity Matrix + for (let i = 0; i < N; i++) { + Sij[i] = Array(N).fill(0) + Sij[i][i] = 1.0 + } } // initial error let Vab = getAij(x) @@ -37,26 +40,29 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c const j = Vab[0][1] psi = getTheta(x[i][i], x[j][j], x[i][j]) x = x1(x, psi, i, j) - Sij = Sij1(Sij, psi, i, j) + if (computeVectors) Sij = Sij1(Sij, psi, i, j) Vab = getAij(x) } const Ei = Array(N).fill(0) // eigenvalues for (let i = 0; i < N; i++) { Ei[i] = x[i][i] } - return sorting(clone(Ei), clone(Sij)) + return sorting(clone(Ei), Sij, computeVectors) } // diagonalization implementation for bigNumber - function diagBig (x, precision) { + function diagBig (x, precision, computeVectors) { const N = x.length const e0 = abs(precision / N) let psi - let Sij = new Array(N) - // Sij is Identity Matrix - for (let i = 0; i < N; i++) { - Sij[i] = Array(N).fill(0) - Sij[i][i] = 1.0 + let Sij + if (computeVectors) { + Sij = new Array(N) + // Sij is Identity Matrix + for (let i = 0; i < N; i++) { + Sij[i] = Array(N).fill(0) + Sij[i][i] = 1.0 + } } // initial error let Vab = getAijBig(x) @@ -65,7 +71,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c const j = Vab[0][1] psi = getThetaBig(x[i][i], x[j][j], x[i][j]) x = x1Big(x, psi, i, j) - Sij = Sij1Big(Sij, psi, i, j) + if (computeVectors) Sij = Sij1Big(Sij, psi, i, j) Vab = getAijBig(x) } const Ei = Array(N).fill(0) // eigenvalues @@ -73,7 +79,7 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c Ei[i] = x[i][i] } // return [clone(Ei), clone(Sij)] - return sorting(clone(Ei), clone(Sij)) + return sorting(clone(Ei), Sij, computeVectors) } // get angle @@ -234,13 +240,15 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c } // sort results - function sorting (E, S) { + function sorting (E, S, computeVectors) { const N = E.length const values = Array(N) - const vecs = Array(N) - - for (let k = 0; k < N; k++) { - vecs[k] = Array(N) + let vecs + if (computeVectors) { + vecs = Array(N) + for (let k = 0; k < N; k++) { + vecs[k] = Array(N) + } } for (let i = 0; i < N; i++) { let minID = 0 @@ -252,11 +260,14 @@ export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, c } } values[i] = E.splice(minID, 1)[0] - for (let k = 0; k < N; k++) { - vecs[i][k] = S[k][minID] - S[k].splice(minID, 1) + if (computeVectors) { + for (let k = 0; k < N; k++) { + vecs[i][k] = S[k][minID] + S[k].splice(minID, 1) + } } } + if (!computeVectors) return { values } const eigenvectors = vecs.map((vector, i) => ({ value: values[i], vector })) return { values, eigenvectors } } diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js index f77862d651..bfdc88878e 100644 --- a/test/unit-tests/function/matrix/eigs.test.js +++ b/test/unit-tests/function/matrix/eigs.test.js @@ -94,12 +94,18 @@ describe('eigs', function () { [1.0, 1.0, 1.0]]).values, [0, 0, 3] ) - approx.deepEqual(eigs( + const sym4 = [[0.6163396801190624, -3.8571699139231796, 2.852995822026198, 4.1957619745869845], [-3.8571699139231796, 0.7047577966772156, 0.9122549659760404, 0.9232933211541949], [2.852995822026198, 0.9122549659760404, 1.6598316026960402, -1.2931270747054358], - [4.1957619745869845, 0.9232933211541949, -1.2931270747054358, -4.665994662426116]]).values, - [-0.9135495807127523, 2.26552473288741, 5.6502090685149735, -8.687249803623432] + [4.1957619745869845, 0.9232933211541949, -1.2931270747054358, -4.665994662426116]] + const fullValues = eigs(sym4).values + approx.deepEqual(fullValues, + [-0.9135495807127523, 2.26552473288741, 5.6502090685149735, -8.687249803623432] + ) + assert.deepStrictEqual( + fullValues, + eigs(sym4, { eigenvectors: false }).values ) }) @@ -147,6 +153,8 @@ describe('eigs', function () { [4.14, 4.27, 3.05, 2.24, 2.73, -4.47]] const ans = eigs(H) const E = ans.values + const justvalues = eigs(H, { eigenvectors: false }) + assert.deepStrictEqual(E, justvalues.values) testEigenvectors(ans, (v, j) => approx.deepEqual(multiply(E[j], v), multiply(H, v)) ) @@ -251,6 +259,8 @@ describe('eigs', function () { [4.24, -4.68, -3.33, 1.67, 2.80, 2.73], [4.14, 4.27, 3.05, 2.24, 2.73, -4.47]]) const ans = eigs(H) + const justvalues = eigs(H, { eigenvectors: false }) + assert.deepStrictEqual(ans.values, justvalues.values) const E = ans.values const Vcols = ans.eigenvectors.map(obj => obj.vector) const V = matrixFromColumns(...Vcols)