mirror of
https://github.com/josdejong/mathjs.git
synced 2026-02-01 16:07:46 +00:00
General eigenproblem algorithm (#1743)
* split the eigs function into multiple algorithms * moved checks and coersions to eigs.js, made them more robust * fix little bugs, make im and re more robust * Implemented matrix balancing algorithm * fix typos * a draft of reduction to Hessenberg matrix * finished implementation of reduction to Hessenberg * fix Hessenberg elimination for complex numbers * implemented non-shifted explicit QR algorithm for real matrices * implemented vector computation, won't work untill usolve is fixed * refactored to match yarn lint * some minor changes * solve merge conflicts * refactored and re-fixed #1789 * some old uncommited changes * fix small problems introduced by merging * done some polishing * improved jsdoc description of eigs * little changes in jsdoc Co-authored-by: Jos de Jong <wjosdejong@gmail.com>
This commit is contained in:
Michal Grňo
GitHub
parent
c23d312db6
commit
9e8deb5c86
@ -240,7 +240,7 @@ export const createNorm = /* #__PURE__ */ factory(
|
||||
const tx = ctranspose(x)
|
||||
const squaredX = multiply(tx, x)
|
||||
const eigenVals = eigs(squaredX).values
|
||||
const rho = eigenVals.get([eigenVals.size()[0] - 1])
|
||||
const rho = eigenVals[eigenVals.length - 1]
|
||||
return abs(sqrt(rho))
|
||||
}
|
||||
|
||||
|
||||
@ -41,6 +41,10 @@ export const createIm = /* #__PURE__ */ factory(name, dependencies, ({ typed })
|
||||
return x.mul(0)
|
||||
},
|
||||
|
||||
Fraction: function (x) {
|
||||
return x.mul(0)
|
||||
},
|
||||
|
||||
Complex: function (x) {
|
||||
return x.im
|
||||
},
|
||||
|
||||
@ -41,6 +41,10 @@ export const createRe = /* #__PURE__ */ factory(name, dependencies, ({ typed })
|
||||
return x
|
||||
},
|
||||
|
||||
Fraction: function (x) {
|
||||
return x
|
||||
},
|
||||
|
||||
Complex: function (x) {
|
||||
return x.re
|
||||
},
|
||||
|
||||
@ -1,373 +1,196 @@
|
||||
import { clone } from '../../utils/object.js'
|
||||
import { factory } from '../../utils/factory.js'
|
||||
import { format } from '../../utils/string.js'
|
||||
import { createComplex } from './eigs/complex.js'
|
||||
import { createRealSymmetric } from './eigs/realSymetric.js'
|
||||
import { typeOf, isNumber, isBigNumber, isComplex, isFraction } from '../../utils/is.js'
|
||||
|
||||
const name = 'eigs'
|
||||
const dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'inv', 'bignumber', 'multiply', 'add']
|
||||
|
||||
export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) => {
|
||||
// The absolute state of math.js's dependency system:
|
||||
const dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'divideScalar', 'inv', 'bignumber', 'multiply', 'add', 'larger', 'column', 'flatten', 'number', 'complex', 'sqrt', 'diag', 'qr', 'usolveAll', 'im', 're', 'smaller', 'round', 'log10', 'transpose']
|
||||
export const createEigs = /* #__PURE__ */ factory(name, dependencies, ({ config, typed, matrix, addScalar, subtract, equal, abs, atan, cos, sin, multiplyScalar, divideScalar, inv, bignumber, multiply, add, larger, column, flatten, number, complex, sqrt, diag, qr, usolveAll, im, re, smaller, round, log10, transpose }) => {
|
||||
const doRealSymetric = createRealSymmetric({ config, addScalar, subtract, column, flatten, equal, abs, atan, cos, sin, multiplyScalar, inv, bignumber, complex, multiply, add })
|
||||
const doComplex = createComplex({ config, addScalar, subtract, multiply, multiplyScalar, flatten, divideScalar, sqrt, abs, bignumber, diag, qr, inv, usolveAll, equal, complex, larger, smaller, round, log10, transpose })
|
||||
|
||||
/**
|
||||
* Compute eigenvalue and eigenvector of a real symmetric matrix.
|
||||
* Only applicable to two dimensional symmetric matrices. Uses Jacobi
|
||||
* Algorithm. Matrix containing mixed type ('number', 'bignumber', 'fraction')
|
||||
* of elements are not supported. Input matrix or 2D array should contain all elements
|
||||
* of either 'number', 'bignumber' or 'fraction' type. For 'number' and 'fraction', the
|
||||
* eigenvalues are of 'number' type. For 'bignumber' the eigenvalues are of ''bignumber' type.
|
||||
* Eigenvectors are always of 'number' type.
|
||||
* Compute eigenvalues and eigenvectors of a matrix. The eigenvalues are sorted by their absolute value, ascending.
|
||||
* An eigenvalue with multiplicity k will be listed k times. The eigenvectors are returned as columns of a matrix –
|
||||
* the eigenvector that belongs to the j-th eigenvalue in the list (eg. `values[j]`) is the j-th column (eg. `column(vectors, j)`).
|
||||
* If the algorithm fails to converge, it will throw an error – in that case, however, you may still find useful information
|
||||
* in `err.values` and `err.vectors`.
|
||||
*
|
||||
* Syntax:
|
||||
*
|
||||
* math.eigs(x)
|
||||
* math.eigs(x, [prec])
|
||||
*
|
||||
* Examples:
|
||||
*
|
||||
* const { eigs, multiply, column, transpose } = math
|
||||
* const H = [[5, 2.3], [2.3, 1]]
|
||||
* const ans = math.eigs(H) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}
|
||||
* const ans = eigs(H) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}
|
||||
* const E = ans.values
|
||||
* const U = ans.vectors
|
||||
* math.multiply(H, math.column(U, 0)) // returns math.multiply(E[0], math.column(U, 0))
|
||||
* const UTxHxU = math.multiply(math.transpose(U), H, U) // rotates H to the eigen-representation
|
||||
* multiply(H, column(U, 0)) // returns multiply(E[0], column(U, 0))
|
||||
* const UTxHxU = multiply(transpose(U), H, U) // diagonalizes H
|
||||
* E[0] == UTxHxU[0][0] // returns true
|
||||
*
|
||||
* See also:
|
||||
*
|
||||
* inv
|
||||
*
|
||||
* @param {Array | Matrix} x Matrix to be diagonalized
|
||||
* @return {{values: Array, vectors: Array} | {values: Matrix, vectors: Matrix}} Object containing eigenvalues (Array or Matrix) and eigenvectors (2D Array/Matrix with eigenvectors as columns).
|
||||
*
|
||||
* @param {number | BigNumber} [prec] Precision, default value: 1e-15
|
||||
* @return {{values: Array, vectors: Matrix}} Object containing an array of eigenvalues and a matrix with eigenvectors as columns.
|
||||
*
|
||||
*/
|
||||
return typed('eigs', {
|
||||
|
||||
Array: function (x) {
|
||||
// check array size
|
||||
const mat = matrix(x)
|
||||
const size = mat.size()
|
||||
if (size.length !== 2 || size[0] !== size[1]) {
|
||||
throw new RangeError('Matrix must be square ' +
|
||||
'(size: ' + format(size) + ')')
|
||||
}
|
||||
|
||||
// use dense 2D matrix implementation
|
||||
const ans = checkAndSubmit(mat, size[0])
|
||||
return { values: ans[0], vectors: ans[1] }
|
||||
return computeValuesAndVectors(mat)
|
||||
},
|
||||
|
||||
Matrix: function (x) {
|
||||
// use dense 2D array implementation
|
||||
// dense matrix
|
||||
const size = x.size()
|
||||
if (size.length !== 2 || size[0] !== size[1]) {
|
||||
throw new RangeError('Matrix must be square ' +
|
||||
'(size: ' + format(size) + ')')
|
||||
}
|
||||
const ans = checkAndSubmit(x, size[0])
|
||||
return { values: matrix(ans[0]), vectors: matrix(ans[1]) }
|
||||
'Array, number|BigNumber': function (x, prec) {
|
||||
const mat = matrix(x)
|
||||
return computeValuesAndVectors(mat, prec)
|
||||
},
|
||||
|
||||
Matrix: function (mat) {
|
||||
return computeValuesAndVectors(mat)
|
||||
},
|
||||
|
||||
'Matrix, number|BigNumber': function (mat, prec) {
|
||||
return computeValuesAndVectors(mat, prec)
|
||||
}
|
||||
})
|
||||
|
||||
// Is the matrix
|
||||
// symmetric ?
|
||||
function isSymmetric (x, n) {
|
||||
for (let i = 0; i < n; i++) {
|
||||
for (let j = i; j < n; j++) {
|
||||
// not symmtric
|
||||
if (!equal(x[i][j], x[j][i])) {
|
||||
throw new TypeError('Input matrix is not symmetric')
|
||||
function computeValuesAndVectors (mat, prec) {
|
||||
if (prec === undefined) {
|
||||
prec = config.epsilon
|
||||
}
|
||||
|
||||
const size = mat.size()
|
||||
|
||||
if (size.length !== 2 || size[0] !== size[1]) {
|
||||
throw new RangeError('Matrix must be square (size: ' + format(size) + ')')
|
||||
}
|
||||
|
||||
const arr = mat.toArray()
|
||||
const N = size[0]
|
||||
|
||||
if (isReal(arr, N, prec)) {
|
||||
coerceReal(arr, N)
|
||||
|
||||
if (isSymmetric(arr, N, prec)) {
|
||||
const type = coerceTypes(mat, arr, N)
|
||||
return doRealSymetric(arr, N, prec, type)
|
||||
}
|
||||
}
|
||||
|
||||
const type = coerceTypes(mat, arr, N)
|
||||
return doComplex(arr, N, prec, type)
|
||||
}
|
||||
|
||||
/** @return {boolean} */
|
||||
function isSymmetric (arr, N, prec) {
|
||||
for (let i = 0; i < N; i++) {
|
||||
for (let j = i; j < N; j++) {
|
||||
// TODO proper comparison of bignum and frac
|
||||
if (larger(bignumber(abs(subtract(arr[i][j], arr[j][i]))), prec)) {
|
||||
return false
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return true
|
||||
}
|
||||
|
||||
/** @return {boolean} */
|
||||
function isReal (arr, N, prec) {
|
||||
for (let i = 0; i < N; i++) {
|
||||
for (let j = 0; j < N; j++) {
|
||||
// TODO proper comparison of bignum and frac
|
||||
if (larger(bignumber(abs(im(arr[i][j]))), prec)) {
|
||||
return false
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return true
|
||||
}
|
||||
|
||||
function coerceReal (arr, N) {
|
||||
for (let i = 0; i < N; i++) {
|
||||
for (let j = 0; j < N; j++) {
|
||||
arr[i][j] = re(arr[i][j])
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// check input for possible problems
|
||||
// and perform diagonalization efficiently for
|
||||
// specific type of number
|
||||
function checkAndSubmit (x, n) {
|
||||
let type = x.datatype()
|
||||
// type check
|
||||
if (type === undefined) {
|
||||
type = x.getDataType()
|
||||
/** @return {'number' | 'BigNumber' | 'Complex'} */
|
||||
function coerceTypes (mat, arr, N) {
|
||||
/** @type {string} */
|
||||
const type = mat.datatype()
|
||||
|
||||
if (type === 'number' || type === 'BigNumber' || type === 'Complex') {
|
||||
return type
|
||||
}
|
||||
if (type !== 'number' && type !== 'BigNumber' && type !== 'Fraction') {
|
||||
if (type === 'mixed') {
|
||||
throw new TypeError('Mixed matrix element type is not supported')
|
||||
} else {
|
||||
throw new TypeError('Matrix element type not supported (' + type + ')')
|
||||
|
||||
let hasNumber = false
|
||||
let hasBig = false
|
||||
let hasComplex = false
|
||||
|
||||
for (let i = 0; i < N; i++) {
|
||||
for (let j = 0; j < N; j++) {
|
||||
const el = arr[i][j]
|
||||
|
||||
if (isNumber(el) || isFraction(el)) {
|
||||
hasNumber = true
|
||||
} else if (isBigNumber(el)) {
|
||||
hasBig = true
|
||||
} else if (isComplex(el)) {
|
||||
hasComplex = true
|
||||
} else {
|
||||
throw TypeError('Unsupported type in Matrix: ' + typeOf(el))
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (hasBig && hasComplex) {
|
||||
console.warn('Complex BigNumbers not supported, this operation will lose precission.')
|
||||
}
|
||||
|
||||
if (hasComplex) {
|
||||
for (let i = 0; i < N; i++) {
|
||||
for (let j = 0; j < N; j++) {
|
||||
arr[i][j] = complex(arr[i][j])
|
||||
}
|
||||
}
|
||||
|
||||
return 'Complex'
|
||||
}
|
||||
|
||||
if (hasBig) {
|
||||
for (let i = 0; i < N; i++) {
|
||||
for (let j = 0; j < N; j++) {
|
||||
arr[i][j] = bignumber(arr[i][j])
|
||||
}
|
||||
}
|
||||
|
||||
return 'BigNumber'
|
||||
}
|
||||
|
||||
if (hasNumber) {
|
||||
for (let i = 0; i < N; i++) {
|
||||
for (let j = 0; j < N; j++) {
|
||||
arr[i][j] = number(arr[i][j])
|
||||
}
|
||||
}
|
||||
|
||||
return 'number'
|
||||
} else {
|
||||
isSymmetric(x.toArray(), n)
|
||||
throw TypeError('Matrix contains unsupported types only.')
|
||||
}
|
||||
|
||||
// perform efficient calculation for 'numbers'
|
||||
if (type === 'number') {
|
||||
return diag(x.toArray())
|
||||
} else if (type === 'Fraction') {
|
||||
const xArr = x.toArray()
|
||||
// convert fraction to numbers
|
||||
for (let i = 0; i < n; i++) {
|
||||
for (let j = i; j < n; j++) {
|
||||
xArr[i][j] = xArr[i][j].valueOf()
|
||||
xArr[j][i] = xArr[i][j]
|
||||
}
|
||||
}
|
||||
return diag(x.toArray())
|
||||
} else if (type === 'BigNumber') {
|
||||
return diagBig(x.toArray())
|
||||
}
|
||||
}
|
||||
|
||||
// diagonalization implementation for number (efficient)
|
||||
function diag (x) {
|
||||
const N = x.length
|
||||
const e0 = Math.abs(config.epsilon / N)
|
||||
let psi
|
||||
let Sij = new Array(N)
|
||||
// Sij is Identity Matrix
|
||||
for (let i = 0; i < N; i++) {
|
||||
Sij[i] = createArray(N, 0)
|
||||
Sij[i][i] = 1.0
|
||||
}
|
||||
// initial error
|
||||
let Vab = getAij(x)
|
||||
while (Math.abs(Vab[1]) >= Math.abs(e0)) {
|
||||
const i = Vab[0][0]
|
||||
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)
|
||||
Vab = getAij(x)
|
||||
}
|
||||
const Ei = createArray(N, 0) // eigenvalues
|
||||
for (let i = 0; i < N; i++) {
|
||||
Ei[i] = x[i][i]
|
||||
}
|
||||
return sorting(clone(Ei), clone(Sij))
|
||||
}
|
||||
|
||||
// diagonalization implementation for bigNumber
|
||||
function diagBig (x) {
|
||||
const N = x.length
|
||||
const e0 = abs(config.epsilon / N)
|
||||
let psi
|
||||
let Sij = new Array(N)
|
||||
// Sij is Identity Matrix
|
||||
for (let i = 0; i < N; i++) {
|
||||
Sij[i] = createArray(N, 0)
|
||||
Sij[i][i] = 1.0
|
||||
}
|
||||
// initial error
|
||||
let Vab = getAijBig(x)
|
||||
while (abs(Vab[1]) >= abs(e0)) {
|
||||
const i = Vab[0][0]
|
||||
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)
|
||||
Vab = getAijBig(x)
|
||||
}
|
||||
const Ei = createArray(N, 0) // eigenvalues
|
||||
for (let i = 0; i < N; i++) {
|
||||
Ei[i] = x[i][i]
|
||||
}
|
||||
// return [clone(Ei), clone(Sij)]
|
||||
return sorting(clone(Ei), clone(Sij))
|
||||
}
|
||||
|
||||
// get angle
|
||||
function getTheta (aii, ajj, aij) {
|
||||
const denom = (ajj - aii)
|
||||
if (Math.abs(denom) <= config.epsilon) {
|
||||
return Math.PI / 4
|
||||
} else {
|
||||
return 0.5 * Math.atan(2 * aij / (ajj - aii))
|
||||
}
|
||||
}
|
||||
|
||||
// get angle
|
||||
function getThetaBig (aii, ajj, aij) {
|
||||
const denom = subtract(ajj, aii)
|
||||
if (abs(denom) <= config.epsilon) {
|
||||
return bignumber(-1).acos().div(4)
|
||||
} else {
|
||||
return multiplyScalar(0.5, atan(multiply(2, aij, inv(denom))))
|
||||
}
|
||||
}
|
||||
|
||||
// update eigvec
|
||||
function Sij1 (Sij, theta, i, j) {
|
||||
const N = Sij.length
|
||||
const c = Math.cos(theta)
|
||||
const s = Math.sin(theta)
|
||||
const Ski = createArray(N, 0)
|
||||
const Skj = createArray(N, 0)
|
||||
for (let k = 0; k < N; k++) {
|
||||
Ski[k] = c * Sij[k][i] - s * Sij[k][j]
|
||||
Skj[k] = s * Sij[k][i] + c * Sij[k][j]
|
||||
}
|
||||
for (let k = 0; k < N; k++) {
|
||||
Sij[k][i] = Ski[k]
|
||||
Sij[k][j] = Skj[k]
|
||||
}
|
||||
return Sij
|
||||
}
|
||||
// update eigvec for overlap
|
||||
function Sij1Big (Sij, theta, i, j) {
|
||||
const N = Sij.length
|
||||
const c = cos(theta)
|
||||
const s = sin(theta)
|
||||
const Ski = createArray(N, bignumber(0))
|
||||
const Skj = createArray(N, bignumber(0))
|
||||
for (let k = 0; k < N; k++) {
|
||||
Ski[k] = subtract(multiplyScalar(c, Sij[k][i]), multiplyScalar(s, Sij[k][j]))
|
||||
Skj[k] = addScalar(multiplyScalar(s, Sij[k][i]), multiplyScalar(c, Sij[k][j]))
|
||||
}
|
||||
for (let k = 0; k < N; k++) {
|
||||
Sij[k][i] = Ski[k]
|
||||
Sij[k][j] = Skj[k]
|
||||
}
|
||||
return Sij
|
||||
}
|
||||
|
||||
// update matrix
|
||||
function x1Big (Hij, theta, i, j) {
|
||||
const N = Hij.length
|
||||
const c = bignumber(cos(theta))
|
||||
const s = bignumber(sin(theta))
|
||||
const c2 = multiplyScalar(c, c)
|
||||
const s2 = multiplyScalar(s, s)
|
||||
const Aki = createArray(N, bignumber(0))
|
||||
const Akj = createArray(N, bignumber(0))
|
||||
// 2cs Hij
|
||||
const csHij = multiply(bignumber(2), c, s, Hij[i][j])
|
||||
// Aii
|
||||
const Aii = addScalar(subtract(multiplyScalar(c2, Hij[i][i]), csHij), multiplyScalar(s2, Hij[j][j]))
|
||||
const Ajj = add(multiplyScalar(s2, Hij[i][i]), csHij, multiplyScalar(c2, Hij[j][j]))
|
||||
// 0 to i
|
||||
for (let k = 0; k < N; k++) {
|
||||
Aki[k] = subtract(multiplyScalar(c, Hij[i][k]), multiplyScalar(s, Hij[j][k]))
|
||||
Akj[k] = addScalar(multiplyScalar(s, Hij[i][k]), multiplyScalar(c, Hij[j][k]))
|
||||
}
|
||||
// Modify Hij
|
||||
Hij[i][i] = Aii
|
||||
Hij[j][j] = Ajj
|
||||
Hij[i][j] = bignumber(0)
|
||||
Hij[j][i] = bignumber(0)
|
||||
// 0 to i
|
||||
for (let k = 0; k < N; k++) {
|
||||
if (k !== i && k !== j) {
|
||||
Hij[i][k] = Aki[k]
|
||||
Hij[k][i] = Aki[k]
|
||||
Hij[j][k] = Akj[k]
|
||||
Hij[k][j] = Akj[k]
|
||||
}
|
||||
}
|
||||
return Hij
|
||||
}
|
||||
|
||||
// update matrix
|
||||
function x1 (Hij, theta, i, j) {
|
||||
const N = Hij.length
|
||||
const c = Math.cos(theta)
|
||||
const s = Math.sin(theta)
|
||||
const c2 = c * c
|
||||
const s2 = s * s
|
||||
const Aki = createArray(N, 0)
|
||||
const Akj = createArray(N, 0)
|
||||
// Aii
|
||||
const Aii = c2 * Hij[i][i] - 2 * c * s * Hij[i][j] + s2 * Hij[j][j]
|
||||
const Ajj = s2 * Hij[i][i] + 2 * c * s * Hij[i][j] + c2 * Hij[j][j]
|
||||
// 0 to i
|
||||
for (let k = 0; k < N; k++) {
|
||||
Aki[k] = c * Hij[i][k] - s * Hij[j][k]
|
||||
Akj[k] = s * Hij[i][k] + c * Hij[j][k]
|
||||
}
|
||||
// Modify Hij
|
||||
Hij[i][i] = Aii
|
||||
Hij[j][j] = Ajj
|
||||
Hij[i][j] = 0
|
||||
Hij[j][i] = 0
|
||||
// 0 to i
|
||||
for (let k = 0; k < N; k++) {
|
||||
if (k !== i && k !== j) {
|
||||
Hij[i][k] = Aki[k]
|
||||
Hij[k][i] = Aki[k]
|
||||
Hij[j][k] = Akj[k]
|
||||
Hij[k][j] = Akj[k]
|
||||
}
|
||||
}
|
||||
return Hij
|
||||
}
|
||||
|
||||
// get max off-diagonal value from Upper Diagonal
|
||||
function getAij (Mij) {
|
||||
const N = Mij.length
|
||||
let maxMij = 0
|
||||
let maxIJ = [0, 1]
|
||||
for (let i = 0; i < N; i++) {
|
||||
for (let j = i + 1; j < N; j++) {
|
||||
if (Math.abs(maxMij) < Math.abs(Mij[i][j])) {
|
||||
maxMij = Math.abs(Mij[i][j])
|
||||
maxIJ = [i, j]
|
||||
}
|
||||
}
|
||||
}
|
||||
return [maxIJ, maxMij]
|
||||
}
|
||||
|
||||
// get max off-diagonal value from Upper Diagonal
|
||||
function getAijBig (Mij) {
|
||||
const N = Mij.length
|
||||
let maxMij = 0
|
||||
let maxIJ = [0, 1]
|
||||
for (let i = 0; i < N; i++) {
|
||||
for (let j = i + 1; j < N; j++) {
|
||||
if (abs(maxMij) < abs(Mij[i][j])) {
|
||||
maxMij = abs(Mij[i][j])
|
||||
maxIJ = [i, j]
|
||||
}
|
||||
}
|
||||
}
|
||||
return [maxIJ, maxMij]
|
||||
}
|
||||
|
||||
// sort results
|
||||
function sorting (E, S) {
|
||||
const N = E.length
|
||||
const Ef = Array(N)
|
||||
const Sf = Array(N)
|
||||
for (let k = 0; k < N; k++) {
|
||||
Sf[k] = Array(N)
|
||||
}
|
||||
for (let i = 0; i < N; i++) {
|
||||
let minID = 0
|
||||
let minE = E[0]
|
||||
for (let j = 0; j < E.length; j++) {
|
||||
if (E[j] < minE) {
|
||||
minID = j
|
||||
minE = E[minID]
|
||||
}
|
||||
}
|
||||
Ef[i] = E.splice(minID, 1)[0]
|
||||
for (let k = 0; k < N; k++) {
|
||||
Sf[k][i] = S[k][minID]
|
||||
S[k].splice(minID, 1)
|
||||
}
|
||||
}
|
||||
return [clone(Ef), clone(Sf)]
|
||||
}
|
||||
|
||||
/**
|
||||
* Create an array of a certain size and fill all items with an initial value
|
||||
* @param {number} size
|
||||
* @param {number} value
|
||||
* @return {number[]}
|
||||
*/
|
||||
function createArray (size, value) {
|
||||
// TODO: as soon as all browsers support Array.fill, use that instead (IE doesn't support it)
|
||||
const array = new Array(size)
|
||||
|
||||
for (let i = 0; i < size; i++) {
|
||||
array[i] = value
|
||||
}
|
||||
|
||||
return array
|
||||
}
|
||||
})
|
||||
|
||||
576
src/function/matrix/eigs/complex.js
Normal file
576
src/function/matrix/eigs/complex.js
Normal file
@ -0,0 +1,576 @@
|
||||
import { clone } from '../../../utils/object.js'
|
||||
|
||||
export function createComplex ({ addScalar, subtract, flatten, multiply, multiplyScalar, divideScalar, sqrt, abs, bignumber, diag, inv, qr, usolveAll, equal, complex, larger, smaller, round, log10, transpose }) {
|
||||
/**
|
||||
* @param {number[][]} arr the matrix to find eigenvalues of
|
||||
* @param {number} N size of the matrix
|
||||
* @param {number|BigNumber} prec precision, anything lower will be considered zero
|
||||
* @param {'number'|'BigNumber'|'Complex'} type
|
||||
* @param {boolean} findVectors should we find eigenvectors?
|
||||
*/
|
||||
function main (arr, N, prec, type, findVectors) {
|
||||
if (findVectors === undefined) {
|
||||
findVectors = true
|
||||
}
|
||||
|
||||
// TODO check if any row/col are zero except the diagonal
|
||||
|
||||
// make sure corresponding rows and columns have similar magnitude
|
||||
// important because of numerical stability
|
||||
const R = balance(arr, N, prec, type, findVectors)
|
||||
|
||||
// R is the row transformation matrix
|
||||
// A' = R A R⁻¹, A is the original matrix
|
||||
// (if findVectors is false, R is undefined)
|
||||
|
||||
// TODO if magnitudes of elements vary over many orders,
|
||||
// move greatest elements to the top left corner
|
||||
|
||||
// using similarity transformations, reduce the matrix
|
||||
// to Hessenberg form (upper triangular plus one subdiagonal row)
|
||||
// updates the transformation matrix R with new row operationsq
|
||||
reduceToHessenberg(arr, N, prec, type, findVectors, R)
|
||||
|
||||
// find eigenvalues
|
||||
let { values, C } = iterateUntilTriangular(arr, N, prec, type, findVectors)
|
||||
|
||||
// values is the list of eigenvalues, C is the column
|
||||
// transformation matrix that transforms the hessenberg
|
||||
// matrix to upper triangular
|
||||
|
||||
// compose transformations A → hess. and hess. → triang.
|
||||
C = multiply(inv(R), C)
|
||||
|
||||
let vectors
|
||||
|
||||
if (findVectors) {
|
||||
vectors = findEigenvectors(arr, N, C, values, prec)
|
||||
vectors = transpose((vectors)) // vectors are columns of a matrix
|
||||
}
|
||||
|
||||
return { values, vectors }
|
||||
}
|
||||
|
||||
/**
|
||||
* @param {number[][]} arr
|
||||
* @param {number} N
|
||||
* @param {number} prec
|
||||
* @param {'number'|'BigNumber'|'Complex'} type
|
||||
* @returns {number[][]}
|
||||
*/
|
||||
function balance (arr, N, prec, type, findVectors) {
|
||||
const big = type === 'BigNumber'
|
||||
const cplx = type === 'Complex'
|
||||
|
||||
const zero = big ? bignumber(0) : cplx ? complex(0) : 0
|
||||
const one = big ? bignumber(1) : cplx ? complex(1) : 1
|
||||
|
||||
// base of the floating-point arithmetic
|
||||
const radix = big ? bignumber(10) : 2
|
||||
const radixSq = multiplyScalar(radix, radix)
|
||||
|
||||
// the diagonal transformation matrix R
|
||||
let Rdiag
|
||||
if (findVectors) {
|
||||
Rdiag = Array(N).fill(one)
|
||||
}
|
||||
|
||||
// this isn't the only time we loop thru the matrix...
|
||||
let last = false
|
||||
|
||||
while (!last) {
|
||||
// ...haha I'm joking! unless...
|
||||
last = true
|
||||
|
||||
for (let i = 0; i < N; i++) {
|
||||
// compute the taxicab norm of i-th column and row
|
||||
// TODO optimize for complex numbers
|
||||
let colNorm = zero
|
||||
let rowNorm = zero
|
||||
|
||||
for (let j = 0; j < N; j++) {
|
||||
if (i === j) continue
|
||||
const c = abs(arr[i][j])
|
||||
colNorm = addScalar(colNorm, c)
|
||||
rowNorm = addScalar(rowNorm, c)
|
||||
}
|
||||
|
||||
if (!equal(colNorm, 0) && !equal(rowNorm, 0)) {
|
||||
// find integer power closest to balancing the matrix
|
||||
// (we want to scale only by integer powers of radix,
|
||||
// so that we don't lose any precision due to round-off)
|
||||
|
||||
let f = one
|
||||
let c = colNorm
|
||||
|
||||
const rowDivRadix = divideScalar(rowNorm, radix)
|
||||
const rowMulRadix = multiplyScalar(rowNorm, radix)
|
||||
|
||||
while (smaller(c, rowDivRadix)) {
|
||||
c = multiplyScalar(c, radixSq)
|
||||
f = multiplyScalar(f, radix)
|
||||
}
|
||||
while (larger(c, rowMulRadix)) {
|
||||
c = divideScalar(c, radixSq)
|
||||
f = divideScalar(f, radix)
|
||||
}
|
||||
|
||||
// check whether balancing is needed
|
||||
// condition = (c + rowNorm) / f < 0.95 * (colNorm + rowNorm)
|
||||
const condition = smaller(divideScalar(addScalar(c, rowNorm), f), multiplyScalar(addScalar(colNorm, rowNorm), 0.95))
|
||||
|
||||
// apply balancing similarity transformation
|
||||
if (condition) {
|
||||
// we should loop once again to check whether
|
||||
// another rebalancing is needed
|
||||
last = false
|
||||
|
||||
const g = divideScalar(1, f)
|
||||
|
||||
for (let j = 0; j < N; j++) {
|
||||
if (i === j) {
|
||||
continue
|
||||
}
|
||||
arr[i][j] = multiplyScalar(arr[i][j], f)
|
||||
arr[j][i] = multiplyScalar(arr[j][i], g)
|
||||
}
|
||||
|
||||
// keep track of transformations
|
||||
if (findVectors) {
|
||||
Rdiag[i] = multiplyScalar(Rdiag[i], f)
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// return the diagonal row transformation matrix
|
||||
return diag(Rdiag)
|
||||
}
|
||||
|
||||
/**
|
||||
* @param {number[][]} arr
|
||||
* @param {number} N
|
||||
* @param {number} prec
|
||||
* @param {'number'|'BigNumber'|'Complex'} type
|
||||
* @param {boolean} findVectors
|
||||
* @param {number[][]} R the row transformation matrix that will be modified
|
||||
*/
|
||||
function reduceToHessenberg (arr, N, prec, type, findVectors, R) {
|
||||
const big = type === 'BigNumber'
|
||||
const cplx = type === 'Complex'
|
||||
|
||||
const zero = big ? bignumber(0) : cplx ? complex(0) : 0
|
||||
|
||||
if (big) { prec = bignumber(prec) }
|
||||
|
||||
for (let i = 0; i < N - 2; i++) {
|
||||
// Find the largest subdiag element in the i-th col
|
||||
|
||||
let maxIndex = 0
|
||||
let max = zero
|
||||
|
||||
for (let j = i + 1; j < N; j++) {
|
||||
const el = arr[j][i]
|
||||
if (smaller(abs(max), abs(el))) {
|
||||
max = el
|
||||
maxIndex = j
|
||||
}
|
||||
}
|
||||
|
||||
// This col is pivoted, no need to do anything
|
||||
if (smaller(abs(max), prec)) {
|
||||
continue
|
||||
}
|
||||
|
||||
if (maxIndex !== i + 1) {
|
||||
// Interchange maxIndex-th and (i+1)-th row
|
||||
const tmp1 = arr[maxIndex]
|
||||
arr[maxIndex] = arr[i + 1]
|
||||
arr[i + 1] = tmp1
|
||||
|
||||
// Interchange maxIndex-th and (i+1)-th column
|
||||
for (let j = 0; j < N; j++) {
|
||||
const tmp2 = arr[j][maxIndex]
|
||||
arr[j][maxIndex] = arr[j][i + 1]
|
||||
arr[j][i + 1] = tmp2
|
||||
}
|
||||
|
||||
// keep track of transformations
|
||||
if (findVectors) {
|
||||
const tmp3 = R[maxIndex]
|
||||
R[maxIndex] = R[i + 1]
|
||||
R[i + 1] = tmp3
|
||||
}
|
||||
}
|
||||
|
||||
// Reduce following rows and columns
|
||||
for (let j = i + 2; j < N; j++) {
|
||||
const n = divideScalar(arr[j][i], max)
|
||||
|
||||
if (n === 0) {
|
||||
continue
|
||||
}
|
||||
|
||||
// from j-th row subtract n-times (i+1)th row
|
||||
for (let k = 0; k < N; k++) {
|
||||
arr[j][k] = subtract(arr[j][k], multiplyScalar(n, arr[i + 1][k]))
|
||||
}
|
||||
|
||||
// to (i+1)th column add n-times j-th column
|
||||
for (let k = 0; k < N; k++) {
|
||||
arr[k][i + 1] = addScalar(arr[k][i + 1], multiplyScalar(n, arr[k][j]))
|
||||
}
|
||||
|
||||
// keep track of transformations
|
||||
if (findVectors) {
|
||||
for (let k = 0; k < N; k++) {
|
||||
R[j][k] = subtract(R[j][k], multiplyScalar(n, R[i + 1][k]))
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return R
|
||||
}
|
||||
|
||||
/**
|
||||
* @returns {{values: values, C: Matrix}}
|
||||
* @see Press, Wiliams: Numerical recipes in Fortran 77
|
||||
* @see https://en.wikipedia.org/wiki/QR_algorithm
|
||||
*/
|
||||
function iterateUntilTriangular (A, N, prec, type, findVectors) {
|
||||
const big = type === 'BigNumber'
|
||||
const cplx = type === 'Complex'
|
||||
|
||||
const one = big ? bignumber(1) : cplx ? complex(1) : 1
|
||||
|
||||
if (big) { prec = bignumber(prec) }
|
||||
|
||||
// The Francis Algorithm
|
||||
// The core idea of this algorithm is that doing successive
|
||||
// A' = Q⁺AQ 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
|
||||
// the eigenvalues of its diagonal blocks and we know how to find
|
||||
// eigenvalues of a 2x2 matrix, we know the eigenvalues of A.
|
||||
|
||||
let arr = clone(A)
|
||||
|
||||
// the list of converged eigenvalues
|
||||
const lambdas = []
|
||||
|
||||
// size of arr, which will get smaller as eigenvalues converge
|
||||
let n = N
|
||||
|
||||
// the diagonal of the block-diagonal matrix that turns
|
||||
// converged 2x2 matrices into upper triangular matrices
|
||||
const Sdiag = []
|
||||
|
||||
// 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
|
||||
let Qpartial = findVectors ? diag(Array(n).fill(one)) : undefined
|
||||
|
||||
// last eigenvalue converged before this many steps
|
||||
let lastConvergenceBefore = 0
|
||||
|
||||
while (lastConvergenceBefore <= 100) {
|
||||
lastConvergenceBefore += 1
|
||||
|
||||
// TODO if the convergence is slow, do something clever
|
||||
|
||||
// Perform the factorization
|
||||
|
||||
const k = 0 // TODO set close to an eigenvalue
|
||||
|
||||
for (let i = 0; i < n; i++) {
|
||||
arr[i][i] = subtract(arr[i][i], k)
|
||||
}
|
||||
|
||||
// TODO do an implicit QR transformation
|
||||
const { Q, R } = qr(arr)
|
||||
arr = multiply(R, Q)
|
||||
|
||||
for (let i = 0; i < n; i++) {
|
||||
arr[i][i] = addScalar(arr[i][i], k)
|
||||
}
|
||||
|
||||
// keep track of transformations
|
||||
if (findVectors) {
|
||||
Qpartial = multiply(Qpartial, Q)
|
||||
}
|
||||
|
||||
// The rightmost diagonal element converged to an eigenvalue
|
||||
if (n === 1 || smaller(abs(arr[n - 1][n - 2]), prec)) {
|
||||
lastConvergenceBefore = 0
|
||||
lambdas.push(arr[n - 1][n - 1])
|
||||
|
||||
// keep track of transformations
|
||||
if (findVectors) {
|
||||
Sdiag.unshift([[1]])
|
||||
inflateMatrix(Qpartial, N)
|
||||
Qtotal = multiply(Qtotal, Qpartial)
|
||||
|
||||
if (n > 1) {
|
||||
Qpartial = diag(Array(n - 1).fill(one))
|
||||
}
|
||||
}
|
||||
|
||||
// reduce the matrix size
|
||||
n -= 1
|
||||
arr.pop()
|
||||
for (let i = 0; i < n; i++) {
|
||||
arr[i].pop()
|
||||
}
|
||||
|
||||
// The rightmost diagonal 2x2 block converged
|
||||
} else if (n === 2 || smaller(abs(arr[n - 2][n - 3]), prec)) {
|
||||
lastConvergenceBefore = 0
|
||||
const ll = eigenvalues2x2(
|
||||
arr[n - 2][n - 2], arr[n - 2][n - 1],
|
||||
arr[n - 1][n - 2], arr[n - 1][n - 1]
|
||||
)
|
||||
lambdas.push(...ll)
|
||||
|
||||
// keep track of transformations
|
||||
if (findVectors) {
|
||||
Sdiag.unshift(jordanBase2x2(
|
||||
arr[n - 2][n - 2], arr[n - 2][n - 1],
|
||||
arr[n - 1][n - 2], arr[n - 1][n - 1],
|
||||
ll[0], ll[1], prec, type
|
||||
))
|
||||
inflateMatrix(Qpartial, N)
|
||||
Qtotal = multiply(Qtotal, Qpartial)
|
||||
|
||||
if (n > 2) {
|
||||
Qpartial = diag(Array(n - 2).fill(one))
|
||||
}
|
||||
}
|
||||
|
||||
// reduce the matrix size
|
||||
n -= 2
|
||||
arr.pop()
|
||||
arr.pop()
|
||||
for (let i = 0; i < n; i++) {
|
||||
arr[i].pop()
|
||||
arr[i].pop()
|
||||
}
|
||||
}
|
||||
|
||||
if (n === 0) {
|
||||
break
|
||||
}
|
||||
}
|
||||
|
||||
// standard sorting
|
||||
lambdas.sort((a, b) => +subtract(abs(a), abs(b)))
|
||||
|
||||
// the algorithm didn't converge
|
||||
if (lastConvergenceBefore > 100) {
|
||||
const err = Error('The eigenvalues failed to converge. Only found these eigenvalues: ' + lambdas.join(', '))
|
||||
err.values = lambdas
|
||||
err.vectors = []
|
||||
throw err
|
||||
}
|
||||
|
||||
// combine the overall QR transformation Qtotal with the subsequent
|
||||
// transformation S that turns the diagonal 2x2 blocks to upper triangular
|
||||
const C = findVectors ? multiply(Qtotal, blockDiag(Sdiag, N)) : undefined
|
||||
|
||||
return { values: lambdas, C }
|
||||
}
|
||||
|
||||
/**
|
||||
* @param {Matrix} A original matrix
|
||||
* @param {number} N size of A
|
||||
* @param {Matrix} C column transformation matrix that turns A into upper triangular
|
||||
* @param {number[]} values array of eigenvalues of A
|
||||
* @returns {Matrix[]} eigenvalues
|
||||
*/
|
||||
function findEigenvectors (A, N, C, values, prec) {
|
||||
const Cinv = inv(C)
|
||||
const U = multiply(Cinv, A, C)
|
||||
|
||||
// turn values into a kind of "multiset"
|
||||
// this way it is easier to find eigenvectors
|
||||
const uniqueValues = []
|
||||
const multiplicities = []
|
||||
|
||||
for (const λ of values) {
|
||||
const i = indexOf(uniqueValues, λ, equal)
|
||||
|
||||
// a dirty trick that helps us find more vectors
|
||||
// TODO with iterative algorithm this can be removed
|
||||
const roundedλ = round(λ, subtract(-1, log10(prec)))
|
||||
|
||||
if (i === -1) {
|
||||
uniqueValues.push(roundedλ)
|
||||
multiplicities.push(1)
|
||||
} else {
|
||||
multiplicities[i] += 1
|
||||
}
|
||||
}
|
||||
|
||||
// find eigenvectors by solving U − λE = 0
|
||||
// TODO replace with an iterative eigenvector algorithm
|
||||
// (this one might fail for imprecise eigenvalues)
|
||||
|
||||
const vectors = []
|
||||
const len = uniqueValues.length
|
||||
const b = Array(N).fill(0)
|
||||
const E = diag(Array(N).fill(1))
|
||||
|
||||
// eigenvalues for which usolve failed (due to numerical error)
|
||||
const failedLambdas = []
|
||||
|
||||
for (let i = 0; i < len; i++) {
|
||||
const λ = uniqueValues[i]
|
||||
|
||||
let solutions = usolveAll(subtract(U, multiply(λ, E)), b)
|
||||
solutions = solutions.map(v => multiply(C, v))
|
||||
|
||||
solutions.shift() // ignore the null vector
|
||||
|
||||
// looks like we missed something
|
||||
if (solutions.length < multiplicities[i]) {
|
||||
failedLambdas.push(λ)
|
||||
}
|
||||
|
||||
vectors.push(...solutions.map(v => flatten(v)))
|
||||
}
|
||||
|
||||
if (failedLambdas.length !== 0) {
|
||||
const err = new Error('Failed to find eigenvectors for the following eigenvalues: ' + failedLambdas.join(', '))
|
||||
err.values = values
|
||||
err.vectors = vectors
|
||||
throw err
|
||||
}
|
||||
|
||||
return vectors
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute the eigenvalues of an 2x2 matrix
|
||||
* @return {[number,number]}
|
||||
*/
|
||||
function eigenvalues2x2 (a, b, c, d) {
|
||||
// λ± = ½ trA ± ½ √( tr²A - 4 detA )
|
||||
const trA = addScalar(a, d)
|
||||
const detA = subtract(multiplyScalar(a, d), multiplyScalar(b, c))
|
||||
const x = multiplyScalar(trA, 0.5)
|
||||
const y = multiplyScalar(sqrt(subtract(multiplyScalar(trA, trA), multiplyScalar(4, detA))), 0.5)
|
||||
|
||||
return [addScalar(x, y), subtract(x, y)]
|
||||
}
|
||||
|
||||
/**
|
||||
* For an 2x2 matrix compute the transformation matrix S,
|
||||
* so that SAS⁻¹ 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
|
||||
*/
|
||||
function jordanBase2x2 (a, b, c, d, l1, l2, prec, type) {
|
||||
const big = type === 'BigNumber'
|
||||
const cplx = type === 'Complex'
|
||||
|
||||
const zero = big ? bignumber(0) : cplx ? complex(0) : 0
|
||||
const one = big ? bignumber(1) : cplx ? complex(1) : 1
|
||||
|
||||
// matrix is already upper triangular
|
||||
// return an identity matrix
|
||||
if (smaller(abs(c), prec)) {
|
||||
return [[one, zero], [zero, one]]
|
||||
}
|
||||
|
||||
// matrix is diagonalizable
|
||||
// return its eigenvectors as columns
|
||||
if (larger(abs(subtract(l1, l2)), prec)) {
|
||||
return [[subtract(l1, d), subtract(l2, d)], [c, c]]
|
||||
}
|
||||
|
||||
// matrix is not diagonalizable
|
||||
// compute off-diagonal elements of N = A - λI
|
||||
// N₁₂ = 0 ⇒ S = ( N⃗₁, I⃗₁ )
|
||||
// N₁₂ ≠ 0 ⇒ S = ( N⃗₂, I⃗₂ )
|
||||
|
||||
const na = subtract(a, l1)
|
||||
const nb = subtract(b, l1)
|
||||
const nc = subtract(c, l1)
|
||||
const nd = subtract(d, l1)
|
||||
|
||||
if (smaller(abs(nb), prec)) {
|
||||
return [[na, one], [nc, zero]]
|
||||
} else {
|
||||
return [[nb, zero], [nd, one]]
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Enlarge the matrix from n×n to N×N, setting the new
|
||||
* elements to 1 on diagonal and 0 elsewhere
|
||||
*/
|
||||
function inflateMatrix (arr, N) {
|
||||
// add columns
|
||||
for (let i = 0; i < arr.length; i++) {
|
||||
arr[i].push(...Array(N - arr[i].length).fill(0))
|
||||
}
|
||||
|
||||
// add rows
|
||||
for (let i = arr.length; i < N; i++) {
|
||||
arr.push(Array(N).fill(0))
|
||||
arr[i][i] = 1
|
||||
}
|
||||
|
||||
return arr
|
||||
}
|
||||
|
||||
/**
|
||||
* Create a block-diagonal matrix with the given square matrices on the diagonal
|
||||
* @param {Matrix[] | number[][][]} arr array of matrices to be placed on the diagonal
|
||||
* @param {number} N the size of the resulting matrix
|
||||
*/
|
||||
function blockDiag (arr, N) {
|
||||
const M = []
|
||||
for (let i = 0; i < N; i++) {
|
||||
M[i] = Array(N).fill(0)
|
||||
}
|
||||
|
||||
let I = 0
|
||||
for (const sub of arr) {
|
||||
const n = sub.length
|
||||
|
||||
for (let i = 0; i < n; i++) {
|
||||
for (let j = 0; j < n; j++) {
|
||||
M[I + i][I + j] = sub[i][j]
|
||||
}
|
||||
}
|
||||
|
||||
I += n
|
||||
}
|
||||
|
||||
return M
|
||||
}
|
||||
|
||||
/**
|
||||
* Finds the index of an element in an array using a custom equality function
|
||||
* @template T
|
||||
* @param {Array<T>} arr array in which to search
|
||||
* @param {T} el the element to find
|
||||
* @param {function(T, T): boolean} fn the equality function, first argument is an element of `arr`, the second is always `el`
|
||||
* @returns {number} the index of `el`, or -1 when it's not in `arr`
|
||||
*/
|
||||
function indexOf (arr, el, fn) {
|
||||
for (let i = 0; i < arr.length; i++) {
|
||||
if (fn(arr[i], el)) {
|
||||
return i
|
||||
}
|
||||
}
|
||||
return -1
|
||||
}
|
||||
|
||||
return main
|
||||
}
|
||||
282
src/function/matrix/eigs/realSymetric.js
Normal file
282
src/function/matrix/eigs/realSymetric.js
Normal file
@ -0,0 +1,282 @@
|
||||
import { clone } from '../../../utils/object.js'
|
||||
|
||||
export function createRealSymmetric ({ config, addScalar, subtract, abs, atan, cos, sin, multiplyScalar, inv, bignumber, multiply, add }) {
|
||||
/**
|
||||
* @param {number[] | BigNumber[]} arr
|
||||
* @param {number} N
|
||||
* @param {number} prec
|
||||
* @param {'number' | 'BigNumber'} type
|
||||
*/
|
||||
function main (arr, N, prec = config.epsilon, type) {
|
||||
if (type === 'number') {
|
||||
return diag(arr, prec)
|
||||
}
|
||||
|
||||
if (type === 'BigNumber') {
|
||||
return diagBig(arr, prec)
|
||||
}
|
||||
|
||||
throw TypeError('Unsupported data type: ' + type)
|
||||
}
|
||||
|
||||
// diagonalization implementation for number (efficient)
|
||||
function diag (x, precision) {
|
||||
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] = createArray(N, 0)
|
||||
Sij[i][i] = 1.0
|
||||
}
|
||||
// initial error
|
||||
let Vab = getAij(x)
|
||||
while (Math.abs(Vab[1]) >= Math.abs(e0)) {
|
||||
const i = Vab[0][0]
|
||||
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)
|
||||
Vab = getAij(x)
|
||||
}
|
||||
const Ei = createArray(N, 0) // eigenvalues
|
||||
for (let i = 0; i < N; i++) {
|
||||
Ei[i] = x[i][i]
|
||||
}
|
||||
return sorting(clone(Ei), clone(Sij))
|
||||
}
|
||||
|
||||
// diagonalization implementation for bigNumber
|
||||
function diagBig (x, precision) {
|
||||
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] = createArray(N, 0)
|
||||
Sij[i][i] = 1.0
|
||||
}
|
||||
// initial error
|
||||
let Vab = getAijBig(x)
|
||||
while (abs(Vab[1]) >= abs(e0)) {
|
||||
const i = Vab[0][0]
|
||||
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)
|
||||
Vab = getAijBig(x)
|
||||
}
|
||||
const Ei = createArray(N, 0) // eigenvalues
|
||||
for (let i = 0; i < N; i++) {
|
||||
Ei[i] = x[i][i]
|
||||
}
|
||||
// return [clone(Ei), clone(Sij)]
|
||||
return sorting(clone(Ei), clone(Sij))
|
||||
}
|
||||
|
||||
// get angle
|
||||
function getTheta (aii, ajj, aij) {
|
||||
const denom = (ajj - aii)
|
||||
if (Math.abs(denom) <= config.epsilon) {
|
||||
return Math.PI / 4.0
|
||||
} else {
|
||||
return 0.5 * Math.atan(2.0 * aij / (ajj - aii))
|
||||
}
|
||||
}
|
||||
|
||||
// get angle
|
||||
function getThetaBig (aii, ajj, aij) {
|
||||
const denom = subtract(ajj, aii)
|
||||
if (abs(denom) <= config.epsilon) {
|
||||
return bignumber(-1).acos().div(4)
|
||||
} else {
|
||||
return multiplyScalar(0.5, atan(multiply(2.0, aij, inv(denom))))
|
||||
}
|
||||
}
|
||||
|
||||
// update eigvec
|
||||
function Sij1 (Sij, theta, i, j) {
|
||||
const N = Sij.length
|
||||
const c = Math.cos(theta)
|
||||
const s = Math.sin(theta)
|
||||
const Ski = createArray(N, 0)
|
||||
const Skj = createArray(N, 0)
|
||||
for (let k = 0; k < N; k++) {
|
||||
Ski[k] = c * Sij[k][i] - s * Sij[k][j]
|
||||
Skj[k] = s * Sij[k][i] + c * Sij[k][j]
|
||||
}
|
||||
for (let k = 0; k < N; k++) {
|
||||
Sij[k][i] = Ski[k]
|
||||
Sij[k][j] = Skj[k]
|
||||
}
|
||||
return Sij
|
||||
}
|
||||
// update eigvec for overlap
|
||||
function Sij1Big (Sij, theta, i, j) {
|
||||
const N = Sij.length
|
||||
const c = cos(theta)
|
||||
const s = sin(theta)
|
||||
const Ski = createArray(N, bignumber(0))
|
||||
const Skj = createArray(N, bignumber(0))
|
||||
for (let k = 0; k < N; k++) {
|
||||
Ski[k] = subtract(multiplyScalar(c, Sij[k][i]), multiplyScalar(s, Sij[k][j]))
|
||||
Skj[k] = addScalar(multiplyScalar(s, Sij[k][i]), multiplyScalar(c, Sij[k][j]))
|
||||
}
|
||||
for (let k = 0; k < N; k++) {
|
||||
Sij[k][i] = Ski[k]
|
||||
Sij[k][j] = Skj[k]
|
||||
}
|
||||
return Sij
|
||||
}
|
||||
|
||||
// update matrix
|
||||
function x1Big (Hij, theta, i, j) {
|
||||
const N = Hij.length
|
||||
const c = bignumber(cos(theta))
|
||||
const s = bignumber(sin(theta))
|
||||
const c2 = multiplyScalar(c, c)
|
||||
const s2 = multiplyScalar(s, s)
|
||||
const Aki = createArray(N, bignumber(0))
|
||||
const Akj = createArray(N, bignumber(0))
|
||||
// 2cs Hij
|
||||
const csHij = multiply(bignumber(2), c, s, Hij[i][j])
|
||||
// Aii
|
||||
const Aii = addScalar(subtract(multiplyScalar(c2, Hij[i][i]), csHij), multiplyScalar(s2, Hij[j][j]))
|
||||
const Ajj = add(multiplyScalar(s2, Hij[i][i]), csHij, multiplyScalar(c2, Hij[j][j]))
|
||||
// 0 to i
|
||||
for (let k = 0; k < N; k++) {
|
||||
Aki[k] = subtract(multiplyScalar(c, Hij[i][k]), multiplyScalar(s, Hij[j][k]))
|
||||
Akj[k] = addScalar(multiplyScalar(s, Hij[i][k]), multiplyScalar(c, Hij[j][k]))
|
||||
}
|
||||
// Modify Hij
|
||||
Hij[i][i] = Aii
|
||||
Hij[j][j] = Ajj
|
||||
Hij[i][j] = bignumber(0)
|
||||
Hij[j][i] = bignumber(0)
|
||||
// 0 to i
|
||||
for (let k = 0; k < N; k++) {
|
||||
if (k !== i && k !== j) {
|
||||
Hij[i][k] = Aki[k]
|
||||
Hij[k][i] = Aki[k]
|
||||
Hij[j][k] = Akj[k]
|
||||
Hij[k][j] = Akj[k]
|
||||
}
|
||||
}
|
||||
return Hij
|
||||
}
|
||||
|
||||
// update matrix
|
||||
function x1 (Hij, theta, i, j) {
|
||||
const N = Hij.length
|
||||
const c = Math.cos(theta)
|
||||
const s = Math.sin(theta)
|
||||
const c2 = c * c
|
||||
const s2 = s * s
|
||||
const Aki = createArray(N, 0)
|
||||
const Akj = createArray(N, 0)
|
||||
// Aii
|
||||
const Aii = c2 * Hij[i][i] - 2 * c * s * Hij[i][j] + s2 * Hij[j][j]
|
||||
const Ajj = s2 * Hij[i][i] + 2 * c * s * Hij[i][j] + c2 * Hij[j][j]
|
||||
// 0 to i
|
||||
for (let k = 0; k < N; k++) {
|
||||
Aki[k] = c * Hij[i][k] - s * Hij[j][k]
|
||||
Akj[k] = s * Hij[i][k] + c * Hij[j][k]
|
||||
}
|
||||
// Modify Hij
|
||||
Hij[i][i] = Aii
|
||||
Hij[j][j] = Ajj
|
||||
Hij[i][j] = 0
|
||||
Hij[j][i] = 0
|
||||
// 0 to i
|
||||
for (let k = 0; k < N; k++) {
|
||||
if (k !== i && k !== j) {
|
||||
Hij[i][k] = Aki[k]
|
||||
Hij[k][i] = Aki[k]
|
||||
Hij[j][k] = Akj[k]
|
||||
Hij[k][j] = Akj[k]
|
||||
}
|
||||
}
|
||||
return Hij
|
||||
}
|
||||
|
||||
// get max off-diagonal value from Upper Diagonal
|
||||
function getAij (Mij) {
|
||||
const N = Mij.length
|
||||
let maxMij = 0
|
||||
let maxIJ = [0, 1]
|
||||
for (let i = 0; i < N; i++) {
|
||||
for (let j = i + 1; j < N; j++) {
|
||||
if (Math.abs(maxMij) < Math.abs(Mij[i][j])) {
|
||||
maxMij = Math.abs(Mij[i][j])
|
||||
maxIJ = [i, j]
|
||||
}
|
||||
}
|
||||
}
|
||||
return [maxIJ, maxMij]
|
||||
}
|
||||
|
||||
// get max off-diagonal value from Upper Diagonal
|
||||
function getAijBig (Mij) {
|
||||
const N = Mij.length
|
||||
let maxMij = 0
|
||||
let maxIJ = [0, 1]
|
||||
for (let i = 0; i < N; i++) {
|
||||
for (let j = i + 1; j < N; j++) {
|
||||
if (abs(maxMij) < abs(Mij[i][j])) {
|
||||
maxMij = abs(Mij[i][j])
|
||||
maxIJ = [i, j]
|
||||
}
|
||||
}
|
||||
}
|
||||
return [maxIJ, maxMij]
|
||||
}
|
||||
|
||||
// sort results
|
||||
function sorting (E, S) {
|
||||
const N = E.length
|
||||
const values = Array(N)
|
||||
const vectors = Array(N)
|
||||
|
||||
for (let k = 0; k < N; k++) {
|
||||
vectors[k] = Array(N)
|
||||
}
|
||||
for (let i = 0; i < N; i++) {
|
||||
let minID = 0
|
||||
let minE = E[0]
|
||||
for (let j = 0; j < E.length; j++) {
|
||||
if (abs(E[j]) < abs(minE)) {
|
||||
minID = j
|
||||
minE = E[minID]
|
||||
}
|
||||
}
|
||||
values[i] = E.splice(minID, 1)[0]
|
||||
for (let k = 0; k < N; k++) {
|
||||
vectors[k][i] = S[k][minID]
|
||||
S[k].splice(minID, 1)
|
||||
}
|
||||
}
|
||||
|
||||
return { values, vectors }
|
||||
}
|
||||
|
||||
/**
|
||||
* Create an array of a certain size and fill all items with an initial value
|
||||
* @param {number} size
|
||||
* @param {number} value
|
||||
* @return {number[]}
|
||||
*/
|
||||
function createArray (size, value) {
|
||||
// TODO: as soon as all browsers support Array.fill, use that instead (IE doesn't support it)
|
||||
const array = new Array(size)
|
||||
|
||||
for (let i = 0; i < size; i++) {
|
||||
array[i] = value
|
||||
}
|
||||
|
||||
return array
|
||||
}
|
||||
|
||||
return main
|
||||
}
|
||||
@ -551,7 +551,7 @@ export function generalize (a) {
|
||||
* This method does not validate Array Matrix shape
|
||||
* @param {Array} array
|
||||
* @param {function} typeOf Callback function to use to determine the type of a value
|
||||
* @return string
|
||||
* @return {string}
|
||||
*/
|
||||
export function getArrayDataType (array, typeOf) {
|
||||
let type // to hold type info
|
||||
|
||||
@ -2,22 +2,20 @@ import assert from 'assert'
|
||||
import math from '../../../../src/defaultInstance.js'
|
||||
import approx from '../../../../tools/approx.js'
|
||||
const eigs = math.eigs
|
||||
const complex = math.complex
|
||||
|
||||
describe('eigs', function () {
|
||||
it('should only accept a square matrix', function () {
|
||||
it('only accepts a square matrix', function () {
|
||||
assert.throws(function () { eigs(math.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/)
|
||||
assert.throws(function () { eigs([[1, 2], [4, 5, 6]]) }, /DimensionError: Dimension mismatch/)
|
||||
assert.throws(function () { eigs([4, 5, 6]) }, /Matrix must be square/)
|
||||
assert.throws(function () { eigs(1.0) }, /TypeError: Unexpected type of argument/)
|
||||
assert.throws(function () { eigs('random') }, /TypeError: Unexpected type of argument/)
|
||||
assert.throws(function () { eigs(math.matrix([[1, 2], [2.1, 3]])) }, /Input matrix is not symmetric/)
|
||||
})
|
||||
|
||||
it('should only accept a matrix with valid element type', function () {
|
||||
assert.throws(function () { eigs([['x', 2], [4, 5]]) }, /Mixed matrix element type is not supported/)
|
||||
assert.throws(function () { eigs([[1, 2], [2, '5']]) }, /Mixed matrix element type is not supported/)
|
||||
assert.throws(function () { eigs([['1', '2'], ['2', '5']]) }, /Matrix element type not supported/)
|
||||
it('only accepts a matrix with valid element type', function () {
|
||||
assert.throws(function () { eigs([['x', 2], [4, 5]]) }, /Cannot convert "x" to a number/)
|
||||
})
|
||||
|
||||
it('eigenvalue check for diagonal matrix', function () {
|
||||
@ -28,22 +26,25 @@ describe('eigs', function () {
|
||||
approx.deepEqual(eigs(
|
||||
[[2, 0, 0], [0, 1, 0], [0, 0, 5]]).values, [1, 2, 5]
|
||||
)
|
||||
approx.deepEqual(eigs(
|
||||
[[complex(2, 1), 0, 0], [0, 1, 0], [0, 0, complex(0, 5)]]).values, [complex(1, 0), complex(2, 1), complex(0, 5)]
|
||||
)
|
||||
})
|
||||
|
||||
it('eigenvalue check for 2x2 simple matrix', function () {
|
||||
it('calculates eigenvalues for 2x2 simple matrix', function () {
|
||||
// 2x2 test
|
||||
approx.deepEqual(eigs(
|
||||
[[1, 0.1], [0.1, 1]]).values, [0.9, 1.1]
|
||||
)
|
||||
approx.deepEqual(eigs(
|
||||
math.matrix([[1, 0.1], [0.1, 1]])).values, math.matrix([0.9, 1.1])
|
||||
math.matrix([[1, 0.1], [0.1, 1]])).values, [0.9, 1.1]
|
||||
)
|
||||
approx.deepEqual(eigs(
|
||||
[[5, 2.3], [2.3, 1]]).values, [-0.04795013082563382, 6.047950130825635]
|
||||
)
|
||||
})
|
||||
|
||||
it('eigenvalue check for 3x3 and 4x4 matrix', function () {
|
||||
it('calculates eigenvalues for 3x3 and 4x4 matrix', function () {
|
||||
// 3x3 test and 4x4
|
||||
approx.deepEqual(eigs(
|
||||
[[1.0, 1.0, 1.0],
|
||||
@ -56,7 +57,7 @@ describe('eigs', function () {
|
||||
[-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,
|
||||
[-8.687249803623432, -0.9135495807127523, 2.26552473288741, 5.6502090685149735]
|
||||
[-0.9135495807127523, 2.26552473288741, 5.6502090685149735, -8.687249803623432]
|
||||
)
|
||||
})
|
||||
|
||||
@ -78,7 +79,7 @@ describe('eigs', function () {
|
||||
approx.deepEqual(Ei, E)
|
||||
})
|
||||
|
||||
it('fractions are supported', function () {
|
||||
it('supports fractions', function () {
|
||||
const aij = math.fraction('1/2')
|
||||
approx.deepEqual(eigs(
|
||||
[[aij, aij, aij],
|
||||
@ -88,7 +89,7 @@ describe('eigs', function () {
|
||||
)
|
||||
})
|
||||
|
||||
it('bigNumber diagonalization is supported', function () {
|
||||
it('diagonalizes matrix with bigNumber', function () {
|
||||
const x = [[math.bignumber(1), math.bignumber(0)], [math.bignumber(0), math.bignumber(1)]]
|
||||
approx.deepEqual(eigs(x).values, [math.bignumber(1), math.bignumber(1)])
|
||||
const y = [[math.bignumber(1), math.bignumber(1.0)], [math.bignumber(1.0), math.bignumber(1)]]
|
||||
@ -112,7 +113,7 @@ describe('eigs', function () {
|
||||
approx.deepEqual(Ei, E)
|
||||
})
|
||||
|
||||
it('make sure BigNumbers input is actually calculated with BigNumber precision', function () {
|
||||
it('actually calculates BigNumbers input with BigNumber precision', function () {
|
||||
const B = math.bignumber([
|
||||
[0, 1],
|
||||
[1, 0]
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user