From 7c77bc29d5ea86d165601c50f7d170bf6188df4b Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Wed, 4 Oct 2023 11:28:25 -0700 Subject: [PATCH] feat: Use simple shifts in QR eigenvalue iterations that improve convergence Although we might want to use better shifts in the future, we might just use a library instead. But for now I think this: Resolves #2178. Also responds to the review feedback provided in PR #3037. --- src/function/matrix/eigs.js | 57 ++++++++------- src/function/matrix/eigs/complexEigs.js | 48 +++++++------ test/typescript-tests/testTypes.ts | 2 +- test/unit-tests/function/matrix/eigs.test.js | 75 ++++++++++---------- 4 files changed, 95 insertions(+), 87 deletions(-) diff --git a/src/function/matrix/eigs.js b/src/function/matrix/eigs.js index 472d544719..50edd5cb70 100644 --- a/src/function/matrix/eigs.js +++ b/src/function/matrix/eigs.js @@ -65,25 +65,42 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, // streamlined. It is done because the Matrix object carries some // type information about its entries, and so constructing the matrix // is a roundabout way of doing type detection. - Array: function (x) { return computeValuesAndVectors(matrix(x)) }, + Array: function (x) { return doEigs(matrix(x)) }, 'Array, number|BigNumber': function (x, prec) { - return computeValuesAndVectors(matrix(x), prec) + return doEigs(matrix(x), prec) }, Matrix: function (mat) { - return computeValuesAndVectors(mat, undefined, true) + return doEigs(mat, undefined, true) }, 'Matrix, number|BigNumber': function (mat, prec) { - return computeValuesAndVectors(mat, prec, true) + return doEigs(mat, prec, true) } }) - function computeValuesAndVectors (mat, prec, matricize = false) { + function doEigs (mat, prec, matricize = false) { + const result = computeValuesAndVectors(mat, prec) + if (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') + } + }) + return result + } + + function computeValuesAndVectors (mat, prec) { if (prec === undefined) { prec = config.epsilon } - let result - const arr = mat.toArray() + const arr = mat.toArray() // NOTE: arr is guaranteed to be unaliased + // and so safe to modify in place const asize = mat.size() if (asize.length !== 2 || asize[0] !== asize[1]) { @@ -93,30 +110,16 @@ export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, const N = asize[0] if (isReal(arr, N, prec)) { - coerceReal(arr, N) + coerceReal(arr, N) // modifies arr by side effect if (isSymmetric(arr, N, prec)) { - const type = coerceTypes(mat, arr, N) - result = doRealSymmetric(arr, N, prec, type) + const type = coerceTypes(mat, arr, N) // modifies arr by side effect + return doRealSymmetric(arr, N, prec, type) } } - if (!result) { - const type = coerceTypes(mat, arr, N) - result = doComplexEigs(arr, N, prec, type) - } - if (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') - } - }) - return result + + const type = coerceTypes(mat, arr, N) // modifies arr by side effect + return doComplexEigs(arr, N, prec, type) } /** @return {boolean} */ diff --git a/src/function/matrix/eigs/complexEigs.js b/src/function/matrix/eigs/complexEigs.js index 5bfa6944ef..fc376577f1 100644 --- a/src/function/matrix/eigs/complexEigs.js +++ b/src/function/matrix/eigs/complexEigs.js @@ -23,9 +23,9 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul const R = balance(arr, N, prec, type, findVectors) // R is the row transformation matrix - // arr = A' = R A R⁻¹, A is the original matrix + // arr = A' = R A R^-1, A is the original matrix // (if findVectors is false, R is undefined) - // (And so to return to original matrix: A = R⁻¹ arr R) + // (And so to return to original matrix: A = R^-1 arr R) // TODO if magnitudes of elements vary over many orders, // move greatest elements to the top left corner @@ -35,7 +35,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul // updates the transformation matrix R with new row operationsq // MODIFIES arr by side effect! reduceToHessenberg(arr, N, prec, type, findVectors, R) - // still true that original A = R⁻¹ arr R) + // still true that original A = R^-1 arr R) // find eigenvalues const { values, C } = iterateUntilTriangular(arr, N, prec, type, findVectors) @@ -43,7 +43,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul // values is the list of eigenvalues, C is the column // transformation matrix that transforms arr, the hessenberg // matrix, to upper triangular - // (So U = C⁻¹ arr C and the relationship between current arr + // (So U = C^-1 arr C and the relationship between current arr // and original A is unchanged.) let eigenvectors @@ -254,7 +254,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul // The Francis Algorithm // The core idea of this algorithm is that doing successive - // A' = Q⁺AQ transformations will eventually converge to block- + // A' = QtAQ transformations will eventually converge to block- // upper-triangular with diagonal blocks either 1x1 or 2x2. // The Q here is the one from the QR decomposition, A = QR. // Since the eigenvalues of a block-upper-triangular matrix are @@ -276,7 +276,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul // N×N matrix describing the overall transformation done during the QR algorithm let Qtotal = findVectors ? diag(Array(N).fill(one)) : undefined - // n×n matrix describing the QR transformations done since last convergence + // nxn matrix describing the QR transformations done since last convergence let Qpartial = findVectors ? diag(Array(n).fill(one)) : undefined // last eigenvalue converged before this many steps @@ -289,7 +289,12 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul // Perform the factorization - const k = 0 // TODO set close to an eigenvalue + const k = arr[n - 1][n - 1] // TODO this is apparently a somewhat + // old-fashioned choice; ideally set close to an eigenvalue, or + // perhaps better yet switch to the implicit QR version that is sometimes + // specifically called the "Francis algorithm" that is alluded to + // in the following TODO. (Or perhaps we switch to an independently + // optimized third-party package for the linear algebra operations...) for (let i = 0; i < n; i++) { arr[i][i] = subtract(arr[i][i], k) @@ -411,18 +416,18 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul const uniqueValues = [] const multiplicities = [] - for (const λ of values) { - const i = indexOf(uniqueValues, λ, equal) + for (const lambda of values) { + const i = indexOf(uniqueValues, lambda, equal) if (i === -1) { - uniqueValues.push(λ) + uniqueValues.push(lambda) multiplicities.push(1) } else { multiplicities[i] += 1 } } - // find eigenvectors by solving U − λE = 0 + // find eigenvectors by solving U − lambdaE = 0 // TODO replace with an iterative eigenvector algorithm // (this one might fail for imprecise eigenvalues) @@ -432,8 +437,8 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul const E = diag(Array(N).fill(one)) for (let i = 0; i < len; i++) { - const λ = uniqueValues[i] - const S = subtract(U, multiply(λ, E)) // the characteristic matrix + const lambda = uniqueValues[i] + const S = subtract(U, multiply(lambda, E)) // the characteristic matrix let solutions = usolveAll(S, b) solutions.shift() // ignore the null vector @@ -451,7 +456,8 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul const correction = multiply(inv(R), C) solutions = solutions.map(v => multiply(correction, v)) - vectors.push(...solutions.map(v => ({ value: λ, vector: flatten(v) }))) + vectors.push( + ...solutions.map(v => ({ value: lambda, vector: flatten(v) }))) } return vectors @@ -462,7 +468,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul * @return {[number,number]} */ function eigenvalues2x2 (a, b, c, d) { - // λ± = ½ trA ± ½ √( tr²A - 4 detA ) + // lambda_+- = 1/2 trA +- 1/2 sqrt( tr^2 A - 4 detA ) const trA = addScalar(a, d) const detA = subtract(multiplyScalar(a, d), multiplyScalar(b, c)) const x = multiplyScalar(trA, 0.5) @@ -473,7 +479,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul /** * For an 2x2 matrix compute the transformation matrix S, - * so that SAS⁻¹ is an upper triangular matrix + * so that SAS^-1 is an upper triangular matrix * @return {[[number,number],[number,number]]} * @see https://math.berkeley.edu/~ogus/old/Math_54-05/webfoils/jordan.pdf * @see http://people.math.harvard.edu/~knill/teaching/math21b2004/exhibits/2dmatrices/index.html @@ -498,12 +504,12 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul } // matrix is not diagonalizable - // compute diagonal elements of N = A - λI + // compute diagonal elements of N = A - lambdaI const na = subtract(a, l1) const nd = subtract(d, l1) - // N⃗₂ = 0 ⇒ S = ( N⃗₁, I⃗₁ ) - // N⃗₁ ≠ 0 ⇒ S = ( N⃗₂, I⃗₂ ) + // col(N,2) = 0 implies S = ( col(N,1), e_1 ) + // col(N,2) != 0 implies S = ( col(N,2), e_2 ) if (smaller(abs(b), prec) && smaller(abs(nd), prec)) { return [[na, one], [c, zero]] @@ -513,7 +519,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul } /** - * Enlarge the matrix from n×n to N×N, setting the new + * Enlarge the matrix from nxn to NxN, setting the new * elements to 1 on diagonal and 0 elsewhere */ function inflateMatrix (arr, N) { @@ -656,7 +662,7 @@ export function createComplexEigs ({ addScalar, subtract, flatten, multiply, mul const vectorShape = size(v) for (let w of orthog) { w = reshape(w, vectorShape) // make sure this is just a vector computation - // v := v − (w, v)/∥w∥² w + // v := v − (w, v)/|w|^2 w v = subtract(v, multiply(divideScalar(dot(w, v), dot(w, w)), w)) } diff --git a/test/typescript-tests/testTypes.ts b/test/typescript-tests/testTypes.ts index da9ac1970f..23a1a34fd9 100644 --- a/test/typescript-tests/testTypes.ts +++ b/test/typescript-tests/testTypes.ts @@ -1282,7 +1282,7 @@ Matrices examples assert.strictEqual(math.matrix([1, 2, 3]) instanceof math.Matrix, true) } - // Eigenvalues and eigenvectors (should this be extended?) + // Eigenvalues and eigenvectors { const D = [ [1, 1], diff --git a/test/unit-tests/function/matrix/eigs.test.js b/test/unit-tests/function/matrix/eigs.test.js index 8ce68a5ed6..f77862d651 100644 --- a/test/unit-tests/function/matrix/eigs.test.js +++ b/test/unit-tests/function/matrix/eigs.test.js @@ -4,6 +4,11 @@ import approx from '../../../../tools/approx.js' const { eigs, add, complex, divide, exp, fraction, matrix, matrixFromColumns, multiply, abs, size, transpose, bignumber: bignum, zeros, Matrix, Complex } = math describe('eigs', function () { + // helper to examine eigenvectors + function testEigenvectors (soln, predicate) { + soln.eigenvectors.forEach((ev, i) => predicate(ev.vector, i)) + } + it('only accepts a square matrix', function () { assert.throws(function () { eigs(matrix([[1, 2, 3], [4, 5, 6]])) }, /Matrix must be square/) assert.throws(function () { eigs([[1, 2, 3], [4, 5, 6]]) }, /Matrix must be square/) @@ -16,36 +21,30 @@ describe('eigs', function () { it('follows aiao-mimo', function () { const realSymArray = eigs([[1, 0], [0, 1]]) assert(Array.isArray(realSymArray.values) && typeof realSymArray.values[0] === 'number') - for (let ix = 0; ix < 2; ++ix) { - assert(Array.isArray(realSymArray.eigenvectors[ix].vector)) - } + testEigenvectors(realSymArray, vector => assert(Array.isArray(vector))) assert(typeof realSymArray.eigenvectors[0].vector[0] === 'number') const genericArray = eigs([[0, 1], [-1, 0]]) assert(Array.isArray(genericArray.values) && genericArray.values[0] instanceof Complex) - for (let ix = 0; ix < 2; ++ix) { - assert(Array.isArray(genericArray.eigenvectors[ix].vector)) - } - assert(genericArray.eigenvectors[0].vector[0] instanceof Complex) + testEigenvectors(genericArray, + vector => assert(Array.isArray(vector) && vector[0] instanceof Complex) + ) const realSymMatrix = eigs(matrix([[1, 0], [0, 1]])) assert(realSymMatrix.values instanceof Matrix) assert.deepStrictEqual(size(realSymMatrix.values), matrix([2])) - for (let ix = 0; ix < 2; ++ix) { - assert(realSymMatrix.eigenvectors[ix].vector instanceof Matrix) - assert.deepStrictEqual( - size(realSymMatrix.eigenvectors[ix].vector), - matrix([2])) - } + testEigenvectors(realSymMatrix, vector => { + assert(vector instanceof Matrix) + assert.deepStrictEqual(size(vector), matrix([2])) + }) const genericMatrix = eigs(matrix([[0, 1], [-1, 0]])) assert(genericMatrix.values instanceof Matrix) assert.deepStrictEqual(size(genericMatrix.values), matrix([2])) - for (let ix = 0; ix < 2; ++ix) { - assert(genericMatrix.eigenvectors[ix].vector instanceof Matrix) - assert.deepStrictEqual( - size(genericMatrix.eigenvectors[ix].vector), matrix([2])) - } + testEigenvectors(genericMatrix, vector => { + assert(vector instanceof Matrix) + assert.deepStrictEqual(size(vector), matrix([2])) + }) }) it('only accepts a matrix with valid element type', function () { @@ -148,10 +147,9 @@ describe('eigs', function () { [4.14, 4.27, 3.05, 2.24, 2.73, -4.47]] const ans = eigs(H) const E = ans.values - for (let j = 0; j < 6; j++) { - const v = ans.eigenvectors[j].vector - approx.deepEqual(multiply(E[j], v), multiply(H, v)) - } + testEigenvectors(ans, + (v, j) => approx.deepEqual(multiply(E[j], v), multiply(H, v)) + ) const Vcols = ans.eigenvectors.map(obj => obj.vector) const V = matrixFromColumns(...Vcols) const VtHV = multiply(transpose(V), H, V) @@ -168,11 +166,10 @@ describe('eigs', function () { const cnt = 0.1 const Ath = multiply(exp(multiply(complex(0, 1), -cnt)), A) const Hth = divide(add(Ath, transpose(Ath)), 2) - const { values, eigenvectors } = eigs(Hth) - for (const i of [0, 1, 2]) { - const v = eigenvectors[i].vector - approx.deepEqual(multiply(Hth, v), multiply(values[i], v)) - } + const example = eigs(Hth) + testEigenvectors(example, (v, i) => + approx.deepEqual(multiply(Hth, v), multiply(example.values[i], v)) + ) }) it('supports fractions', function () { @@ -225,17 +222,19 @@ describe('eigs', function () { // equal to 2 which has a unique eigenvector (up to scale, of course). // It is from https://web.uvic.ca/~tbazett/diffyqs/sec_multeigen.html // The iterative eigenvalue calculation currently being used has a - // great deal of difficulty converging. So we can get eigs to produce - // **something** by using a very low precision, but the results are - // basically garbage, so we don't actually check them. - const bad = [[2, 0, 0], [-1, -1, 9], [0, -1, 5]] - const junk = eigs(bad, 1e-4) - assert.strictEqual(junk.values.length, 3) - const junkm = eigs(matrix(bad), 1e-4) - assert.deepStrictEqual(junkm.values.size(), [3]) - // Hopefully in future iterations of mathjs we can ratchet down - // these precision values and actually test the return values for - // correctness. + // great deal of difficulty converging. We can use a fine precision, + // but it still doesn't produce good eigenvalues. Hopefully someday + // we'll be able to get closer. + 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) + // 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 + // only one. + const poorm = eigs(matrix(difficult), 1e-14) + assert.deepStrictEqual(poorm.values.size(), [3]) }) it('diagonalizes matrix with bigNumber', function () {