Skip to content

Commit

Permalink
feat: Use simple shifts in QR eigenvalue iterations that improve conv…
Browse files Browse the repository at this point in the history 
…ergence

  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 josdejong#2178.

  Also responds to the review feedback provided in PR josdejong#3037.
  • Loading branch information
@gwhitney
gwhitney committed Oct 4, 2023
1 parent bc82bfc commit 7c77bc2
Show file tree
Hide file tree
Showing 4 changed files with 95 additions and 87 deletions.
57 changes: 30 additions & 27 deletions src/function/matrix/eigs.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]) {
Expand All @@ -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} */
Expand Down
48 changes: 27 additions & 21 deletions src/function/matrix/eigs/complexEigs.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -35,15 +35,15 @@ 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)

// 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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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]]
Expand All @@ -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) {
Expand Down Expand Up @@ -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))
}

Expand Down
2 changes: 1 addition & 1 deletion test/typescript-tests/testTypes.ts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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],
Expand Down
75 changes: 37 additions & 38 deletions test/unit-tests/function/matrix/eigs.test.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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/)
Expand All @@ -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 () {
Expand Down Expand Up @@ -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)
Expand All @@ -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 () {
Expand Down Expand Up @@ -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 () {
Expand Down

0 comments on commit 7c77bc2

Please sign in to comment.