From e19e5553fa206768a4bf65ccc226e58773ecbeff Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Wed, 4 Oct 2023 13:39:05 -0700 Subject: [PATCH 1/7] 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) From 6d9f759fa9d88c24b7511bb40f85ad3d6d6829df Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Thu, 5 Oct 2023 14:10:51 -0700 Subject: [PATCH 2/7] fix: Add test for precision in options arg of eigs And also a fix for a small bug that the new test uncovered. --- src/function/matrix/eigs/complexEigs.js | 2 +- test/unit-tests/function/matrix/eigs.test.js | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/function/matrix/eigs/complexEigs.js b/src/function/matrix/eigs/complexEigs.js index fc376577f1..dc29620a67 100644 --- a/src/function/matrix/eigs/complexEigs.js +++ b/src/function/matrix/eigs/complexEigs.js @@ -150,7 +150,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul } // return the diagonal row transformation matrix - return diag(Rdiag) + return findVectors ? diag(Rdiag) : null } /** diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js index bfdc88878e..c2ac4427e3 100644 --- a/test/unit-tests/function/matrix/eigs.test.js +++ b/test/unit-tests/function/matrix/eigs.test.js @@ -243,6 +243,9 @@ describe('eigs', function () { // only one. const poorm = eigs(matrix(difficult), 1e-14) assert.deepStrictEqual(poorm.values.size(), [3]) + // Make sure the precision argument can go in the options object + const stillbad = eigs(difficult, { precision: 1e-14, eigenvectors: false }) + assert.deepStrictEqual(stillbad.values, poor.values) }) it('diagonalizes matrix with bigNumber', function () { From c53a943c365b428013a4d0889325072305d09226 Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Thu, 5 Oct 2023 17:07:35 -0700 Subject: [PATCH 3/7] test: check eigs with matrix and options --- test/unit-tests/function/matrix/eigs.test.js | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js index c2ac4427e3..7cbd250679 100644 --- a/test/unit-tests/function/matrix/eigs.test.js +++ b/test/unit-tests/function/matrix/eigs.test.js @@ -30,13 +30,17 @@ describe('eigs', function () { vector => assert(Array.isArray(vector) && vector[0] instanceof Complex) ) - const realSymMatrix = eigs(matrix([[1, 0], [0, 1]])) + const id2 = matrix([[1, 0], [0, 1]]) + const realSymMatrix = eigs(id2) assert(realSymMatrix.values instanceof Matrix) assert.deepStrictEqual(size(realSymMatrix.values), matrix([2])) testEigenvectors(realSymMatrix, vector => { assert(vector instanceof Matrix) assert.deepStrictEqual(size(vector), matrix([2])) }) + // Check we get exact values in this trivial case with lower precision + const rough = eigs(id2, { precision: 1e-6 }) + assert.deepStrictEqual(realSymMatrix, rough) const genericMatrix = eigs(matrix([[0, 1], [-1, 0]])) assert(genericMatrix.values instanceof Matrix) From d68b28486a541649258da9c707094e3e9358972a Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Thu, 5 Oct 2023 17:13:47 -0700 Subject: [PATCH 4/7] refactor: remove dead code from complexEigs.js --- src/function/matrix/eigs/complexEigs.js | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/function/matrix/eigs/complexEigs.js b/src/function/matrix/eigs/complexEigs.js index dc29620a67..17d053f887 100644 --- a/src/function/matrix/eigs/complexEigs.js +++ b/src/function/matrix/eigs/complexEigs.js @@ -10,11 +10,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul * * @returns {{ values: number[], vectors: number[][] }} */ - function complexEigs (arr, N, prec, type, findVectors) { - if (findVectors === undefined) { - findVectors = true - } - + function complexEigs (arr, N, prec, type, findVectors = true) { // TODO check if any row/col are zero except the diagonal // make sure corresponding rows and columns have similar magnitude From cf60b8c6585192e4671adf9fcbdc5530d8223efd Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Sat, 14 Oct 2023 11:13:50 -0700 Subject: [PATCH 5/7] fix: add new signatures of eigs to typescript --- test/typescript-tests/testTypes.ts | 7 +++++++ types/index.d.ts | 10 ++++++++-- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/test/typescript-tests/testTypes.ts b/test/typescript-tests/testTypes.ts index 23a1a34fd9..c9ea0d5a3b 100644 --- a/test/typescript-tests/testTypes.ts +++ b/test/typescript-tests/testTypes.ts @@ -1291,6 +1291,13 @@ Matrices examples const eig = math.eigs(D) assert.ok(math.deepEqual(eig.values, [1, 1])) assert.deepStrictEqual(eig.eigenvectors, [{ value: 1, vector: [1, 0] }]) + const eigvv = math.eigs(D, { precision: 1e-6 }) + assert.ok(math.deepEqual(eigvv.values, [1, 1])) + assert.deepStrictEqual(eigvv.eigenvectors, [{ value: 1, vector: [1, 0] }]) + const eigv = math.eigs(D, { eigenvectors: false }) + assert.ok(math.deepEqual(eigv.values, [1, 1])) + //@ts-expect-error ...verify that eigenvectors not expected to be there + eigv.eigenvectors } // Fourier transform and inverse diff --git a/types/index.d.ts b/types/index.d.ts index a013d1c2d5..c91a465267 100644 --- a/types/index.d.ts +++ b/types/index.d.ts @@ -1763,7 +1763,10 @@ declare namespace math { */ eigs( x: MathCollection, - prec?: number | BigNumber + opts?: + | number + | BigNumber + | { precision?: number | BigNumber; eigenvectors?: true } ): { values: MathCollection eigenvectors: { @@ -1771,7 +1774,10 @@ declare namespace math { vector: MathCollection }[] } - + eigs( + x: MathCollection, + opts: { eigenvectors: false; precision?: number | BigNumber } + ): { values: MathCollection } /** * Compute the matrix exponential, expm(A) = e^A. The matrix must be * square. Not to be confused with exp(a), which performs element-wise From a5ebd0c212cf875b03471fbdf89b048783915043 Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Wed, 18 Oct 2023 11:36:28 -0700 Subject: [PATCH 6/7] test: ensure eigenvectors property not present with eigenvectors: false option --- src/function/matrix/eigs/complexEigs.js | 7 +++---- test/unit-tests/function/matrix/eigs.test.js | 10 ++++++---- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/function/matrix/eigs/complexEigs.js b/src/function/matrix/eigs/complexEigs.js index 17d053f887..26dcc0d38b 100644 --- a/src/function/matrix/eigs/complexEigs.js +++ b/src/function/matrix/eigs/complexEigs.js @@ -42,13 +42,12 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul // (So U = C^-1 arr C and the relationship between current arr // and original A is unchanged.) - let eigenvectors - if (findVectors) { - eigenvectors = findEigenvectors(arr, N, C, R, values, prec, type) + const eigenvectors = findEigenvectors(arr, N, C, R, values, prec, type) + return { values, eigenvectors } } - return { values, eigenvectors } + return { values } } /** diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js index 7cbd250679..764a6d90b9 100644 --- a/test/unit-tests/function/matrix/eigs.test.js +++ b/test/unit-tests/function/matrix/eigs.test.js @@ -107,10 +107,9 @@ describe('eigs', function () { approx.deepEqual(fullValues, [-0.9135495807127523, 2.26552473288741, 5.6502090685149735, -8.687249803623432] ) - assert.deepStrictEqual( - fullValues, - eigs(sym4, { eigenvectors: false }).values - ) + const justEigs = eigs(sym4, { eigenvectors: false }) + assert.deepStrictEqual(fullValues, justEigs.values) + assert.ok(!('eigenvectors' in justEigs)) }) it('calculates eigenvalues and eigenvectors for 5x5 matrix', function () { @@ -162,6 +161,7 @@ describe('eigs', function () { testEigenvectors(ans, (v, j) => approx.deepEqual(multiply(E[j], v), multiply(H, v)) ) + assert.ok(!('eigenvectors' in justvalues)) const Vcols = ans.eigenvectors.map(obj => obj.vector) const V = matrixFromColumns(...Vcols) const VtHV = multiply(transpose(V), H, V) @@ -250,6 +250,7 @@ describe('eigs', function () { // Make sure the precision argument can go in the options object const stillbad = eigs(difficult, { precision: 1e-14, eigenvectors: false }) assert.deepStrictEqual(stillbad.values, poor.values) + assert.ok(!('eigenvectors' in stillbad)) }) it('diagonalizes matrix with bigNumber', function () { @@ -268,6 +269,7 @@ describe('eigs', function () { const ans = eigs(H) const justvalues = eigs(H, { eigenvectors: false }) assert.deepStrictEqual(ans.values, justvalues.values) + assert.ok(!('eigenvectors' in justvalues)) const E = ans.values const Vcols = ans.eigenvectors.map(obj => obj.vector) const V = matrixFromColumns(...Vcols) From b4cfafebff8c11ffc11c978ee0e64e20529a0849 Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Wed, 18 Oct 2023 12:30:06 -0700 Subject: [PATCH 7/7] fix: correct balancing code in complexEigs --- src/function/matrix/eigs/complexEigs.js | 11 +++++------ test/unit-tests/function/matrix/eigs.test.js | 2 +- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/src/function/matrix/eigs/complexEigs.js b/src/function/matrix/eigs/complexEigs.js index 26dcc0d38b..8c0b2eb224 100644 --- a/src/function/matrix/eigs/complexEigs.js +++ b/src/function/matrix/eigs/complexEigs.js @@ -90,9 +90,8 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul for (let j = 0; j < N; j++) { if (i === j) continue - const c = abs(arr[i][j]) // should be real - colNorm = addScalar(colNorm, c) - rowNorm = addScalar(rowNorm, c) + colNorm = addScalar(colNorm, abs(arr[j][i])) + rowNorm = addScalar(rowNorm, abs(arr[i][j])) } if (!equal(colNorm, 0) && !equal(rowNorm, 0)) { @@ -131,13 +130,13 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul if (i === j) { continue } - arr[i][j] = multiplyScalar(arr[i][j], f) - arr[j][i] = multiplyScalar(arr[j][i], g) + arr[i][j] = multiplyScalar(arr[i][j], g) + arr[j][i] = multiplyScalar(arr[j][i], f) } // keep track of transformations if (findVectors) { - Rdiag[i] = multiplyScalar(Rdiag[i], f) + Rdiag[i] = multiplyScalar(Rdiag[i], g) } } } diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js index 764a6d90b9..9850224677 100644 --- a/test/unit-tests/function/matrix/eigs.test.js +++ b/test/unit-tests/function/matrix/eigs.test.js @@ -240,7 +240,7 @@ describe('eigs', function () { const difficult = [[2, 0, 0], [-1, -1, 9], [0, -1, 5]] const poor = eigs(difficult, 1e-14) assert.strictEqual(poor.values.length, 3) - approx.deepEqual(poor.values, [2, 2, 2], 6e-6) + approx.deepEqual(poor.values, [2, 2, 2], 7e-6) // Note the eigenvectors are junk, so we don't test them. The function // eigs thinks there are three of them, for example. Hopefully some // future iteration of mathjs will be able to discover there is really