diff --git a/AUTHORS b/AUTHORS index dd5006d29b..6a38bc1974 100644 --- a/AUTHORS +++ b/AUTHORS @@ -284,5 +284,6 @@ anslem chibuike <144047596+AnslemHack@users.noreply.github.com> Ayomide Bamigbade Anadian Dheemanth D <165369664+Dheemanth07@users.noreply.github.com> +Jackson Loper # Generated by tools/update-authors.js diff --git a/src/function/algebra/decomposition/schur.js b/src/function/algebra/decomposition/schur.js index 1beebd567c..c69d7c30d2 100644 --- a/src/function/algebra/decomposition/schur.js +++ b/src/function/algebra/decomposition/schur.js @@ -1,3 +1,4 @@ +import { arraySize, clone } from '../../../utils/array.js' import { factory } from '../../../utils/factory.js' const name = 'schur' @@ -6,9 +7,16 @@ const dependencies = [ 'matrix', 'identity', 'multiply', - 'qr', - 'norm', - 'subtract' + 'abs', + 'isZero', + 'isPositive', + 'isNegative', + 'equalScalar', + 'addScalar', + 'divideScalar', + 'multiplyScalar', + 'subtractScalar', + 'sqrt' ] export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( @@ -17,15 +25,31 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( matrix, identity, multiply, - qr, - norm, - subtract + abs, + isZero, + isPositive, + isNegative, + equalScalar, + addScalar, + divideScalar, + multiplyScalar, + subtractScalar, + sqrt } ) => { /** * - * Performs a real Schur decomposition of the real matrix A = UTU' where U is orthogonal - * and T is upper quasi-triangular. + * Performs a real Schur decomposition of the real matrix A = UTU' where + * U is orthogonal and T is upper quasi-triangular. + * + * Real Schur decomposition: For a real square matrix A, returns orthogonal + * U (which is to say, U * U' = I), and quasi-upper-triangular T such that + * A = U*T*U'. In more detail, T is block upper triangular with 1x1 and 2x2 + * blocks on the diagonal. 1x1 blocks correspond to real eigenvalues and + * 2x2 blocks correspond to complex conjugate eigenvalue pairs. + * The two matrices are returned as a plain JavaScript object with properties + * 'T' and 'U' whose values are the respective matrices. + * * https://en.wikipedia.org/wiki/Schur_decomposition * * Syntax: @@ -35,7 +59,7 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( * Examples: * * const A = [[1, 0], [-4, 3]] - * math.schur(A) // returns {T: [[3, 4], [0, 1]], R: [[0, 1], [-1, 0]]} + * math.schur(A) // returns {T: [[1, 0], [-4, 3]], U: [[1, 0], [0, 1]]} * * See also: * @@ -45,33 +69,452 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( * @return {{U: Array | Matrix, T: Array | Matrix}} Object containing both matrix U and T of the Schur Decomposition A=UTU' */ return typed(name, { - Array: function (X) { - const r = _schur(matrix(X)) - return { - U: r.U.valueOf(), - T: r.T.valueOf() + Array: X => _schur(X, arraySize(X)), + Matrix: X => { + const { U, T } = _schur(X.toArray(), X.size()) + return { U: matrix(U), T: matrix(T) } + } + }) + + /** + * Main Schur decomposition function using Francis QR algorithm + */ + function _schur (X, size) { + if (size.length !== 2 || size[0] !== size[1]) { + throw new RangeError('Matrix must be square') + } + + const n = size[0] + + // Handle trivial cases + if (n === 0) return { U: [], T: [] } + if (n === 1) return { U: [[1]], T: clone(X) } + + // Step 1: Reduce to upper Hessenberg form + // This is a similarity transformation: H = P' * A * P + const { H, P } = reduceToHessenberg(X, n) + + // Step 2: Apply Francis QR algorithm to get quasi-triangular form + // This computes the Schur form: T = Q' * H * Q + const { T, Q } = francisQR(H, n) + + // Step 3: Combine transformations: U = P * Q + // So that A = U * T * U' + return { U: multiply(P, Q), T } + } + + /** + * Reduce matrix to upper Hessenberg form using Householder reflections. + * Returns H (upper Hessenberg) and P (orthogonal) such that H = P' * A * P + */ + function reduceToHessenberg (arr, n) { + const H = clone(arr) + // P will accumulate the orthogonal transformation + const P = identity(n, '') + + for (let k = 0; k < n - 2; k++) { + // Compute Householder vector for column k, rows k+1 to n-1 + const x = [] + for (let i = k + 1; i < n; i++) { + x.push(H[i][k]) + } + + const { v, beta } = computeHouseholderVector(x) + if (!v) continue // Column is already zero, skip + + // Apply Householder reflection from the left: H := (I - beta*v*v') * H + // Only affects rows k+1 to n-1 + for (let j = k; j < n; j++) { + let sum = false + const Hcolj = [] + for (let i = 0; i < v.length; i++) { + const Hentry = H[k + 1 + i][j] + Hcolj.push(Hentry) + const term = multiplyScalar(v[i], Hentry) + if (sum !== false) sum = addScalar(sum, term) + else sum = term + } + sum = multiplyScalar(sum, beta) + for (let i = 0; i < v.length; i++) { + H[k + 1 + i][j] = subtractScalar(Hcolj[i], multiplyScalar(v[i], sum)) + } } - }, - Matrix: function (X) { - return _schur(X) + // Apply Householder reflection from the right: H := H * (I - beta*v*v') + // Affects all rows, columns k+1 to n-1 + for (let i = 0; i < n; i++) { + const Hrow = H[i] + let sum = false + for (let j = 0; j < v.length; j++) { + const term = multiplyScalar(Hrow[k + 1 + j], v[j]) + if (sum !== false) sum = addScalar(sum, term) + else sum = term + } + sum = multiplyScalar(sum, beta) + for (let j = 0; j < v.length; j++) { + H[i][k + 1 + j] = subtractScalar( + Hrow[k + 1 + j], multiplyScalar(sum, v[j])) + } + } + + // Accumulate P: P := P * (I - beta*v*v') + for (let i = 0; i < n; i++) { + const Prow = P[i] + let sum = false + for (let j = 0; j < v.length; j++) { + const term = multiplyScalar(Prow[k + 1 + j], v[j]) + if (sum !== false) sum = addScalar(sum, term) + else sum = term + } + sum = multiplyScalar(sum, beta) + for (let j = 0; j < v.length; j++) { + Prow[k + 1 + j] = subtractScalar( + Prow[k + 1 + j], multiplyScalar(sum, v[j])) + } + } } - }) - function _schur (X) { - const n = X.size()[0] - let A = X - let U = identity(n) - let k = 0 - let A0 - do { - A0 = A - const QR = qr(A) - const Q = QR.Q - const R = QR.R - A = multiply(R, Q) - U = multiply(U, Q) - if ((k++) > 100) { break } - } while (norm(subtract(A, A0)) > 1e-4) - return { U, T: A } + + // Clean up small subdiagonal entries (should be zero due to Householder) + const entry = H[0][0] + const zero = subtractScalar(entry, entry) + for (let i = 2; i < n; i++) { + for (let j = 0; j < i - 1; j++) { + H[i][j] = zero + } + } + + return { H, P } + } + + /** + * Compute Householder vector v and scalar beta such that + * (I - beta * v * v') * x = ||x|| * e_1 + */ + function computeHouseholderVector (x) { + const m = x.length + if (m < 2) return { v: null, beta: null } + + let sigma = multiplyScalar(x[1], x[1]) + for (let i = 2; i < m; i++) { + sigma = addScalar(sigma, multiplyScalar(x[i], x[i])) + } + + const x0 = x[0] + const x0sq = multiplyScalar(x0, x0) + + // If the vector is already a multiple of e_1 (sigma ≈ 0), + // no transformation needed + if (isZero(sigma)) return { v: null, beta: null } + + const normX = sqrt(addScalar(x0sq, sigma)) + + // Choose sign to avoid cancellation + const zero = subtractScalar(x0, x0) + const v0 = isPositive(x0) + ? divideScalar(subtractScalar(zero, sigma), addScalar(x0, normX)) + : subtractScalar(x0, normX) + + const v0sq = multiplyScalar(v0, v0) + const one = divideScalar(x0, x0) + const two = addScalar(one, one) + const beta = divideScalar(two, addScalar(one, divideScalar(sigma, v0sq))) + + // Construct v = [1, x[1]/v0, x[2]/v0, ...] + const v = [one] + for (let i = 1; i < m; i++) { + v.push(divideScalar(x[i], v0)) + } + + return { v, beta } + } + + /** + * Francis QR algorithm with implicit double shift for upper + * Hessenberg matrices. + * Computes the real Schur form T and orthogonal Q such that T = Q' * H * Q + */ + function francisQR (Hin, n) { + const H = clone(Hin) + // Q accumulates the orthogonal transformations + const Q = identity(n, '') + + const zero = subtractScalar(H[0][0], H[0][0]) + const maxIterationsPerEigenvalue = 30 // max iterations per eigenvalue + const maxTotalIterations = 30 * n // safety limit + + let p = n - 1 // Index of last unconverged eigenvalue + let iterCount = 0 + let totalIter = 0 + + while (p > 0 && totalIter < maxTotalIterations) { + totalIter++ + let q = p - 1 + + // Find the largest q such that H[q][q-1] is negligible + while (q > 0) { + const scale = addScalar(abs(H[q - 1][q - 1]), abs(H[q][q])) + if (equalScalar(scale, addScalar(scale, H[q][q - 1]))) { + H[q][q - 1] = zero + break + } + q-- + } + + // q is the start of the unreduced block [q, p] + // If q == p, we have a 1x1 block (real eigenvalue) + if (q === p) { + p-- + iterCount = 0 + continue + } + + // If q == p-1, we have a 2x2 block + if (q === p - 1) { + // Check if eigenvalues are complex + const a = H[p - 1][p - 1] + const b = H[p - 1][p] + const c = H[p][p - 1] + const d = H[p][p] + + // Discriminant of characteristic polynomial λ² - (a+d)λ + (ad-bc) + // discriminant = (a+d)² - 4(ad-bc) = (a-d)² + 4bc + // If discriminant < 0, eigenvalues are complex conjugates + const diff = subtractScalar(a, d) + const discriminant = addScalar( + multiplyScalar(diff, diff), multiplyScalar(4, multiplyScalar(b, c))) + + if (isNegative(discriminant)) { + // Complex eigenvalues - keep the 2x2 block + p -= 2 + iterCount = 0 + continue + } + + // Real eigenvalues in a 2x2 block + // If we've tried many times and it won't split, accept it + if (iterCount > maxIterationsPerEigenvalue) { + p -= 2 + iterCount = 0 + continue + } + } + + // Perform Francis double shift QR step + iterCount++ + + // Apply exceptional shift if convergence is slow + if (iterCount === 10 || iterCount === 20) { + // Exceptional shift: use a random-ish perturbation based on subdiagonal elements + const subdiagVal = abs(H[p][p - 1]) + const prevSubdiag = (p >= 2) ? abs(H[p - 1][p - 2]) : zero + const shift = multiplyScalar( + addScalar(subdiagVal, prevSubdiag), iterCount === 10 ? 1.5 : -1.5) + for (let i = q; i <= p; i++) { + H[i][i] = addScalar(H[i][i], shift) + } + francisStep(H, Q, n, q, p) + for (let i = q; i <= p; i++) { + H[i][i] = subtractScalar(H[i][i], shift) + } + } else { + francisStep(H, Q, n, q, p) + } + } + + // Clean up tiny subdiagonal elements + for (let i = 1; i < n; i++) { + // Use a relative threshold based on nearby diagonal elements + const scale = addScalar(abs(H[i - 1][i - 1]), abs(H[i][i])) + if (equalScalar(scale, addScalar(scale, H[i][i - 1]))) { + H[i][i - 1] = zero + } + } + + return { T: H, Q } + } + + /** + * Perform one Francis double shift QR step on the active block [q, p] of H. + * This implements the implicit double shift QR iteration. + */ + function francisStep (H, Q, n, q, p) { + // Compute the Wilkinson shift from the bottom 2x2 submatrix + // The shift is chosen as the eigenvalue of the 2x2 block closest to H[p][p] + const a = H[p - 1][p - 1] + const b = H[p - 1][p] + const c = H[p][p - 1] + const d = H[p][p] + + const zero = subtractScalar(a, a) + + // Compute the eigenvalues of the 2x2 matrix [[a,b],[c,d]] + // trace and determinant + const trace = addScalar(a, d) + const det = subtractScalar(multiplyScalar(a, d), multiplyScalar(b, c)) + + // For the implicit double shift, we use both eigenvalues + // First column of (H - s1*I)(H - s2*I) where s1, s2 are the eigenvalues + // = H^2 - trace*H + det*I + // First column is [H^2]_0 - trace*H_0 + det*e_0 for active block + + // Compute first column of H^2 - trace*H + det*I + // (restricted to active block) + // For row q: sum_k H[q][k]*H[k][q] - trace*H[q][q] + det + // Since H is Hessenberg, H[k][q] = 0 for k > q+1 + + const Hqq = H[q][q] + const Hqq1 = H[q][q + 1] + const Hq1q = H[q + 1][q] + const Hq1q1 = H[q + 1][q + 1] + + // First column of M = H^2 - trace*H + det*I + // M[q][q] = H[q][q]*H[q][q] + H[q][q+1]*H[q+1][q] - trace*H[q][q] + det + let x = addScalar( + addScalar( + multiplyScalar(Hqq, Hqq), + multiplyScalar(Hqq1, Hq1q) + ), + subtractScalar(det, multiplyScalar(trace, Hqq)) + ) + + // M[q+1][q] = H[q+1][q]*H[q][q] + H[q+1][q+1]*H[q+1][q] - trace*H[q+1][q] + // = H[q+1][q] * (H[q][q] + H[q+1][q+1] - trace) + // = H[q+1][q] * (H[q][q] + H[q+1][q+1] - a - d) + let y = multiplyScalar( + Hq1q, + subtractScalar(addScalar(Hqq, Hq1q1), trace) + ) + + // M[q+2][q] = H[q+2][q+1]*H[q+1][q] (only nonzero element from Hessenberg structure) + let z = zero + if (q + 2 <= p) { + z = multiplyScalar(H[q + 2][q + 1], Hq1q) + } + + // Perform bulge chasing: apply Householder reflections to eliminate the bulge + for (let k = q; k <= p - 1; k++) { + // Determine the size of the Householder reflection (3 or 2) + const r = Math.min(3, p - k + 1) + + // Compute Householder vector for [x, y, z] or [x, y] + let householder + if (r === 3) { + householder = computeHouseholderVector3(x, y, z) + } else { + householder = computeHouseholderVector2(x, y) + } + + if (householder === null) { + // Small values, skip this iteration + if (k < p - 1) { + x = H[k + 1][k] + y = H[k + 2][k] + z = (k + 3 <= p) ? H[k + 3][k] : 0 + } + continue + } + + const { v, beta } = householder + + // Apply Householder reflection from the left + // TODO: Unify this code with the similar code in reduceToHessenberg + const jStart = Math.max(0, k - 1) + for (let j = jStart; j < n; j++) { + let sum = false + const Hcolj = [] + for (let i = 0; i < r; i++) { + const Hentry = H[k + i][j] + Hcolj.push(Hentry) + const term = multiplyScalar(v[i], Hentry) + if (sum !== false) sum = addScalar(sum, term) + else sum = term + } + sum = multiplyScalar(sum, beta) + for (let i = 0; i < r; i++) { + H[k + i][j] = subtractScalar(Hcolj[i], multiplyScalar(v[i], sum)) + } + } + + // Apply Householder reflection from the right + const iEnd = Math.min(n, k + r + 1) + for (let i = 0; i < iEnd; i++) { + const Hrow = H[i] + let sum = 0 + for (let j = 0; j < r; j++) { + sum = addScalar(sum, multiplyScalar(Hrow[k + j], v[j])) + } + sum = multiplyScalar(sum, beta) + for (let j = 0; j < r; j++) { + H[i][k + j] = subtractScalar(Hrow[k + j], multiplyScalar(sum, v[j])) + } + } + + // Accumulate Q + for (let i = 0; i < n; i++) { + const Qrow = Q[i] + let sum = multiplyScalar(Qrow[k], v[0]) + for (let j = 1; j < r; j++) { + sum = addScalar(sum, multiplyScalar(Q[i][k + j], v[j])) + } + sum = multiplyScalar(sum, beta) + for (let j = 0; j < r; j++) { + Qrow[k + j] = subtractScalar(Qrow[k + j], multiplyScalar(sum, v[j])) + } + } + + // Prepare for next iteration + if (k < p - 1) { + x = H[k + 1][k] + y = H[k + 2][k] + z = (k + 3 <= p) ? H[k + 3][k] : 0 + } + } + } + + /** + * Compute Householder vector for 3-element vector [x, y, z] + * Can this be unified with computeHouseholderVector? + */ + function computeHouseholderVector3 (x, y, z) { + const norm = sqrt(addScalar(addScalar( + multiplyScalar(x, x), + multiplyScalar(y, y) + ), multiplyScalar(z, z))) + + if (isZero(norm)) return null + + const zero = subtractScalar(x, x) + // Choose sign to avoid cancellation + const signedNorm = isNegative(x) ? subtractScalar(zero, norm) : norm + const u0 = addScalar(x, signedNorm) + + const one = divideScalar(norm, norm) + const v = [one, divideScalar(y, u0), divideScalar(z, u0)] + const vNormSq = addScalar( + addScalar(one, multiplyScalar(v[1], v[1])), multiplyScalar(v[2], v[2])) + const beta = divideScalar(addScalar(one, one), vNormSq) + + return { v, beta } + } + + /** + * Compute Householder vector for 2-element vector [x, y] + */ + function computeHouseholderVector2 (x, y) { + const norm = sqrt(addScalar(multiplyScalar(x, x), multiplyScalar(y, y))) + + if (isZero(norm)) return null + + const zero = subtractScalar(x, x) + // Choose sign to avoid cancellation + const signedNorm = isNegative(x) ? subtractScalar(zero, norm) : norm + const u0 = addScalar(x, signedNorm) + + const one = divideScalar(norm, norm) + const v = [one, divideScalar(y, u0)] + const vNormSq = addScalar(one, multiplyScalar(v[1], v[1])) + const beta = divideScalar(addScalar(one, one), vNormSq) + + return { v, beta } } }) diff --git a/src/utils/array.js b/src/utils/array.js index 3c8b8813fa..e15054b6ee 100644 --- a/src/utils/array.js +++ b/src/utils/array.js @@ -951,5 +951,5 @@ export function deepForEach (array, callback, skipIndex = false) { * @returns {Array} cloned array */ export function clone (array) { - return Object.assign([], array) + return array.map(elt => Array.isArray(elt) ? clone(elt) : elt) } diff --git a/test/node-tests/doc.test.js b/test/node-tests/doc.test.js index aa2ca84027..49306c4c34 100644 --- a/test/node-tests/doc.test.js +++ b/test/node-tests/doc.test.js @@ -82,10 +82,18 @@ function extractValue (spec) { try { value = eval(spec) // eslint-disable-line no-eval } catch (err) { - if (spec[0] === '[') { - // maybe it was an array with mathjs expressions in it + if ('[{'.includes(spec[0])) { + // maybe it was an array or object with mathjs expressions in it try { - value = math.evaluate(spec).toArray() + value = math.evaluate(spec) + if (spec[0] === '[') value = value.toArray() + else { + for (const key in value) { + if (math.isMatrix(value[key])) { + value[key] = value[key].toArray() + } + } + } } catch (newError) { value = spec } @@ -113,7 +121,7 @@ const knownProblems = new Set([ 'rotate', 'reshape', 'partitionSelect', 'matrixFromFunction', 'matrixFromColumns', 'getMatrixDataType', 'eigs', 'diff', 'slu', 'rationalize', 'qr', 'lusolve', 'lup', 'derivative', - 'symbolicEqual', 'schur', 'sylvester', 'freqz', 'round', + 'symbolicEqual', 'sylvester', 'freqz', 'round', 'import', 'typed', 'unit', 'sparse', 'matrix', 'index', 'bignumber', 'fraction', 'complex', 'parse' diff --git a/test/unit-tests/function/algebra/decomposition/schur.test.js b/test/unit-tests/function/algebra/decomposition/schur.test.js index 7af3c8b426..054a879a6b 100644 --- a/test/unit-tests/function/algebra/decomposition/schur.test.js +++ b/test/unit-tests/function/algebra/decomposition/schur.test.js @@ -3,62 +3,147 @@ import assert from 'assert' import math from '../../../../../src/defaultInstance.js' +/** + * Helper function to verify Schur decomposition properties: + * 1. A = U*T*U' (decomposition is accurate) + * 2. U is orthogonal (U*U' = I) + * 3. T is quasi-upper-triangular (lower triangular elements are zero, except + * for 2x2 blocks on diagonal which represent complex conjugate eigenvalues) + */ +function verifySchurDecomposition (A, result, tolerance = 1e-10) { + const { U, T } = result + const n = Array.isArray(A) ? A.length : A.size()[0] + + // Verify A = U*T*U' + const reconstructed = math.multiply(math.multiply(U, T), math.transpose(U)) + const errorNorm = math.norm(math.subtract(A, reconstructed)) + assert.ok(errorNorm < tolerance, `Decomposition error too large: ${errorNorm}`) + + // Verify U is orthogonal: U*U' = I + const UUT = math.multiply(U, math.transpose(U)) + const orthogonalError = math.norm(math.subtract(UUT, math.identity(n))) + assert.ok(orthogonalError < tolerance, `U is not orthogonal: ${orthogonalError}`) + + // Verify T is quasi-upper-triangular + // Elements below the first subdiagonal must be zero + // Elements on the first subdiagonal may be non-zero only for 2x2 blocks + const Tarr = Array.isArray(T) ? T : T.valueOf() + for (let i = 2; i < n; i++) { + for (let j = 0; j < i - 1; j++) { + assert.ok( + Math.abs(Tarr[i][j]) < tolerance, + `T[${i}][${j}] = ${Tarr[i][j]} should be zero (quasi-upper-triangular)` + ) + } + } +} + describe('schur', function () { it('should calculate schur decomposition of order 5 Array with numbers', function () { - assert.ok(math.norm(math.subtract(math.schur([ - [-5.3, -1.4, -0.2, 0.7, 1.0], - [-0.4, -1.0, -0.1, -1.2, 0.7], - [0.3, 0.7, -2.5, 0.7, -0.3], - [3.6, -0.1, 1.4, -2.4, 0.3], - [2.8, 0.7, 1.4, 0.5, -4.8] - ]).T, [ - [-6.97747746169558, 0.5046853036888738, -0.5269551218982134, 2.9902479419087253, -2.2914719859941908], - [7.686296479877504e-28, -3.667202530573731, -1.3776362163231233, 0.4680921120934126, -0.374760141366345], - [-4.92421882171775e-28, 0.0859443307361262, -3.798852922420336, 0.6595326982269121, 0.38704916773017245], - [-1.8971150504442668e-71, -1.3481560449785515e-44, -8.960727665382592e-44, -1.3867791434503591, 0.5989746088924175], - [5.3038070596554604e-164, -7.015716167009369e-138, 1.0415178925287816e-136, 3.3306222609902884e-93, -0.16968794185997693] - ])) < 1e-3) - assert.ok(math.norm(math.subtract(math.schur([ + const A = [ [-5.3, -1.4, -0.2, 0.7, 1.0], [-0.4, -1.0, -0.1, -1.2, 0.7], [0.3, 0.7, -2.5, 0.7, -0.3], [3.6, -0.1, 1.4, -2.4, 0.3], [2.8, 0.7, 1.4, 0.5, -4.8] - ]).U, [ - [0.6039524392527362, -0.11955248228665324, 0.5309978859071411, -0.3239619623824945, 0.48377530651243317], - [0.034196874004165004, -0.3407725193032822, -0.05878741847605931, -0.7490924342670748, -0.5640117270489927], - [-0.024973926752867345, 0.6355774530247909, -0.45835474572889245, -0.5151114453511857, 0.3463938944685342], - [-0.42238909820980175, -0.6350073886392834, -0.26854304887535935, -0.13970317163522103, 0.5715948922294664], - [-0.6745633977804205, 0.24977632323107185, 0.6575567219332444, -0.22147775183161147, 0.03380493473031572] - ])) < 1e-3) + ] + const result = math.schur(A) + + // Verify decomposition properties + verifySchurDecomposition(A, result) + + // Verify result types + assert.ok(Array.isArray(result.T)) + assert.ok(Array.isArray(result.U)) }) it('should calculate schur decomposition of order 5 Matrix with numbers', function () { - assert.ok(math.norm(math.subtract(math.schur(math.matrix([ - [-5.3, -1.4, -0.2, 0.7, 1.0], - [-0.4, -1.0, -0.1, -1.2, 0.7], - [0.3, 0.7, -2.5, 0.7, -0.3], - [3.6, -0.1, 1.4, -2.4, 0.3], - [2.8, 0.7, 1.4, 0.5, -4.8] - ])).T, math.matrix([ - [-6.97747746169558, 0.5046853036888738, -0.5269551218982134, 2.9902479419087253, -2.2914719859941908], - [7.686296479877504e-28, -3.667202530573731, -1.3776362163231233, 0.4680921120934126, -0.374760141366345], - [-4.92421882171775e-28, 0.0859443307361262, -3.798852922420336, 0.6595326982269121, 0.38704916773017245], - [-1.8971150504442668e-71, -1.3481560449785515e-44, -8.960727665382592e-44, -1.3867791434503591, 0.5989746088924175], - [5.3038070596554604e-164, -7.015716167009369e-138, 1.0415178925287816e-136, 3.3306222609902884e-93, -0.16968794185997693] - ]))) < 1e-3) - assert.ok(math.norm(math.subtract(math.schur(math.matrix([ + const A = math.matrix([ [-5.3, -1.4, -0.2, 0.7, 1.0], [-0.4, -1.0, -0.1, -1.2, 0.7], [0.3, 0.7, -2.5, 0.7, -0.3], [3.6, -0.1, 1.4, -2.4, 0.3], [2.8, 0.7, 1.4, 0.5, -4.8] - ])).U, math.matrix([ - [0.6039524392527362, -0.11955248228665324, 0.5309978859071411, -0.3239619623824945, 0.48377530651243317], - [0.034196874004165004, -0.3407725193032822, -0.05878741847605931, -0.7490924342670748, -0.5640117270489927], - [-0.024973926752867345, 0.6355774530247909, -0.45835474572889245, -0.5151114453511857, 0.3463938944685342], - [-0.42238909820980175, -0.6350073886392834, -0.26854304887535935, -0.13970317163522103, 0.5715948922294664], - [-0.6745633977804205, 0.24977632323107185, 0.6575567219332444, -0.22147775183161147, 0.03380493473031572] - ]))) < 1e-3) + ]) + const result = math.schur(A) + + // Verify decomposition properties + verifySchurDecomposition(A, result) + + // Verify result types are matrices + assert.ok(math.isMatrix(result.T)) + assert.ok(math.isMatrix(result.U)) + }) + + it('should handle 2x2 matrix', function () { + const A = [[1, 2], [3, 4]] + const result = math.schur(A) + verifySchurDecomposition(A, result) + }) + + it('should handle 1x1 matrix', function () { + const A = [[5]] + const result = math.schur(A) + assert.deepStrictEqual(result.T, [[5]]) + assert.deepStrictEqual(result.U, [[1]]) + }) + + it('should handle identity matrix', function () { + const A = math.identity(3) + const result = math.schur(A) + verifySchurDecomposition(A, result) + }) + + it('should handle orthogonal/rotation matrix with complex eigenvalues', function () { + // This matrix has complex eigenvalues and was previously problematic + const A = math.matrix([ + [-0.03591206220229135, -0.09100469507870354, 0.9952027277203429], + [-0.3802068171617618, -0.9197139315803332, -0.09782157349362577], + [0.9242040358990549, -0.38189583596928406, -0.0015717815434243287] + ]) + const result = math.schur(A) + + // Verify decomposition properties + verifySchurDecomposition(A, result) + + // Verify T has a 2x2 block for complex eigenvalue pair + const T = result.T.valueOf() + // T should have form: + // [real eigenvalue, *, *] + // [0, 2x2 block] + // [0, 2x2 block] + // The 2x2 block will have a non-zero T[2][1] element + assert.ok(Math.abs(T[1][0]) < 1e-10, 'T[1][0] should be zero') + assert.ok(Math.abs(T[2][0]) < 1e-10, 'T[2][0] should be zero') + // T[2][1] should be non-zero (complex eigenvalue 2x2 block) + // The magnitude should be significant relative to the matrix norm + const matrixNorm = math.norm(A) + assert.ok(Math.abs(T[2][1]) > 1e-10 * matrixNorm, 'T[2][1] should be non-zero for complex eigenvalue block') + }) + + it('should handle symmetric matrix', function () { + const A = [ + [4, 2, 2], + [2, 5, 1], + [2, 1, 6] + ] + const result = math.schur(A) + verifySchurDecomposition(A, result) + }) + + it('should handle diagonal matrix', function () { + const A = [ + [1, 0, 0], + [0, 2, 0], + [0, 0, 3] + ] + const result = math.schur(A) + verifySchurDecomposition(A, result) + }) + + it('should throw error for non-square matrix', function () { + assert.throws(function () { + math.schur([[1, 2, 3], [4, 5, 6]]) + }, /Matrix must be square/) }) })