From 44513ec104e2e8df6696c3c73b10f39fd3327259 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 11 Jan 2026 12:35:24 +0000 Subject: [PATCH 1/6] Initial plan From cff49454fb5db9dd0ee659f74d2a529f036c3ab1 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 11 Jan 2026 12:44:23 +0000 Subject: [PATCH 2/6] Implement robust Schur decomposition with Francis QR algorithm - Add Hessenberg reduction using Householder reflections - Implement Francis double shift QR iteration for convergence - Handle 2x2 blocks for complex conjugate eigenvalue pairs - Add proper convergence detection and exceptional shifts - Update tests to verify decomposition properties (A=UTU', orthogonality, quasi-triangularity) - Add new tests for edge cases including rotation matrices with complex eigenvalues Co-authored-by: jacksonloper <2056977+jacksonloper@users.noreply.github.com> --- src/function/algebra/decomposition/schur.js | 476 +++++++++++++++++- .../algebra/decomposition/schur.test.js | 167 ++++-- 2 files changed, 583 insertions(+), 60 deletions(-) diff --git a/src/function/algebra/decomposition/schur.js b/src/function/algebra/decomposition/schur.js index 1beebd567c..927883b57d 100644 --- a/src/function/algebra/decomposition/schur.js +++ b/src/function/algebra/decomposition/schur.js @@ -8,7 +8,14 @@ const dependencies = [ 'multiply', 'qr', 'norm', - 'subtract' + 'subtract', + 'abs', + 'addScalar', + 'divideScalar', + 'multiplyScalar', + 'subtractScalar', + 'sqrt', + 'transpose' ] export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( @@ -19,13 +26,27 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( multiply, qr, norm, - subtract + subtract, + abs, + addScalar, + divideScalar, + multiplyScalar, + subtractScalar, + sqrt, + transpose } ) => { /** * * 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 and + * quasi-upper-triangular T such that A = U*T*U'. + * T is block upper triangular with 1x1 and 2x2 blocks on the diagonal. + * 1x1 blocks correspond to real eigenvalues, 2x2 blocks correspond to + * complex conjugate eigenvalue pairs. + * * https://en.wikipedia.org/wiki/Schur_decomposition * * Syntax: @@ -35,7 +56,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 {U: [[0, 1], [-1, 0]], T: [[3, 4], [0, 1]]} * * See also: * @@ -57,21 +78,440 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( return _schur(X) } }) + + /** + * Main Schur decomposition function using Francis QR algorithm + */ 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 } + const size = 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: matrix([]), T: matrix([]) } + } + if (n === 1) { + return { U: identity(1), T: X.clone() } + } + + // Convert to 2D array for internal processing + const arr = X.toArray() + + // Step 1: Reduce to upper Hessenberg form + // This is a similarity transformation: H = P' * A * P + const { H, P } = reduceToHessenberg(arr, 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' + const U = multiply(matrix(P), matrix(Q)) + + return { U, T: matrix(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) { + // Clone the array to avoid modifying the original + const H = arr.map(row => [...row]) + + // P will accumulate the orthogonal transformation + const P = [] + for (let i = 0; i < n; i++) { + P[i] = Array(n).fill(0) + P[i][i] = 1 + } + + 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 householder = computeHouseholderVector(x) + if (householder === null) { + continue // Column is already zero, skip + } + + const { v, beta } = householder + + // 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 = 0 + for (let i = 0; i < v.length; i++) { + sum = addScalar(sum, multiplyScalar(v[i], H[k + 1 + i][j])) + } + sum = multiplyScalar(sum, beta) + for (let i = 0; i < v.length; i++) { + H[k + 1 + i][j] = subtractScalar(H[k + 1 + i][j], multiplyScalar(v[i], sum)) + } + } + + // 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++) { + let sum = 0 + for (let j = 0; j < v.length; j++) { + sum = addScalar(sum, multiplyScalar(H[i][k + 1 + j], v[j])) + } + sum = multiplyScalar(sum, beta) + for (let j = 0; j < v.length; j++) { + H[i][k + 1 + j] = subtractScalar(H[i][k + 1 + j], multiplyScalar(sum, v[j])) + } + } + + // Accumulate P: P := P * (I - beta*v*v') + for (let i = 0; i < n; i++) { + let sum = 0 + for (let j = 0; j < v.length; j++) { + sum = addScalar(sum, multiplyScalar(P[i][k + 1 + j], v[j])) + } + sum = multiplyScalar(sum, beta) + for (let j = 0; j < v.length; j++) { + P[i][k + 1 + j] = subtractScalar(P[i][k + 1 + j], multiplyScalar(sum, v[j])) + } + } + } + + // Clean up small subdiagonal entries (should be zero due to Householder) + for (let i = 2; i < n; i++) { + for (let j = 0; j < i - 1; j++) { + H[i][j] = 0 + } + } + + 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 === 0) return null + + let sigma = 0 + for (let i = 1; 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, no transformation needed + if (abs(sigma) < 1e-14 && abs(x0) < 1e-14) { + return null + } + if (abs(sigma) < 1e-14) { + return null + } + + const normX = sqrt(addScalar(x0sq, sigma)) + + // Choose sign to avoid cancellation + let v0 + if (x0 <= 0) { + v0 = subtractScalar(x0, normX) + } else { + v0 = divideScalar(-sigma, addScalar(x0, normX)) + } + + const v0sq = multiplyScalar(v0, v0) + const beta = divideScalar(2, addScalar(1, divideScalar(sigma, v0sq))) + + // Construct v = [1, x[1]/v0, x[2]/v0, ...] + const v = [1] + 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 = Hin.map(row => [...row]) + + // Q accumulates the orthogonal transformations + const Q = [] + for (let i = 0; i < n; i++) { + Q[i] = Array(n).fill(0) + Q[i][i] = 1 + } + + const eps = Number.EPSILON // machine epsilon for convergence + 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 threshold = eps * (abs(H[q - 1][q - 1]) + abs(H[q][q])) + if (abs(H[q][q - 1]) <= threshold) { + H[q][q - 1] = 0 + 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: (a-d)^2 + 4bc (derived from characteristic polynomial) + const diff = subtractScalar(a, d) + const discriminant = addScalar(multiplyScalar(diff, diff), multiplyScalar(4, multiplyScalar(b, c))) + + if (discriminant < 0) { + // 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 + const shift = multiplyScalar(addScalar(abs(H[p][p - 1]), abs(H[p - 1][p - 2] || 0)), 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 threshold = eps * (abs(H[i - 1][i - 1]) + abs(H[i][i])) + if (abs(H[i][i - 1]) <= threshold) { + H[i][i - 1] = 0 + } + } + + 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] + + // 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 = 0 + 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 + const jStart = Math.max(0, k - 1) + for (let j = jStart; j < n; j++) { + let sum = 0 + for (let i = 0; i < r; i++) { + sum = addScalar(sum, multiplyScalar(v[i], H[k + i][j])) + } + sum = multiplyScalar(sum, beta) + for (let i = 0; i < r; i++) { + H[k + i][j] = subtractScalar(H[k + i][j], 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++) { + let sum = 0 + for (let j = 0; j < r; j++) { + sum = addScalar(sum, multiplyScalar(H[i][k + j], v[j])) + } + sum = multiplyScalar(sum, beta) + for (let j = 0; j < r; j++) { + H[i][k + j] = subtractScalar(H[i][k + j], multiplyScalar(sum, v[j])) + } + } + + // Accumulate Q + for (let i = 0; i < n; i++) { + let sum = 0 + for (let j = 0; j < r; j++) { + sum = addScalar(sum, multiplyScalar(Q[i][k + j], v[j])) + } + sum = multiplyScalar(sum, beta) + for (let j = 0; j < r; j++) { + Q[i][k + j] = subtractScalar(Q[i][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] + */ + function computeHouseholderVector3 (x, y, z) { + const norm = sqrt(addScalar(addScalar( + multiplyScalar(x, x), + multiplyScalar(y, y) + ), multiplyScalar(z, z))) + + if (abs(norm) < 1e-14) { + return null + } + + // Choose sign to avoid cancellation + const s = x >= 0 ? 1 : -1 + const u0 = addScalar(x, multiplyScalar(s, norm)) + + const v = [1, divideScalar(y, u0), divideScalar(z, u0)] + const vNormSq = addScalar(addScalar(1, multiplyScalar(v[1], v[1])), multiplyScalar(v[2], v[2])) + const beta = divideScalar(2, 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 (abs(norm) < 1e-14) { + return null + } + + // Choose sign to avoid cancellation + const s = x >= 0 ? 1 : -1 + const u0 = addScalar(x, multiplyScalar(s, norm)) + + const v = [1, divideScalar(y, u0)] + const vNormSq = addScalar(1, multiplyScalar(v[1], v[1])) + const beta = divideScalar(2, vNormSq) + + return { v, beta } } }) diff --git a/test/unit-tests/function/algebra/decomposition/schur.test.js b/test/unit-tests/function/algebra/decomposition/schur.test.js index 7af3c8b426..d94fa939d5 100644 --- a/test/unit-tests/function/algebra/decomposition/schur.test.js +++ b/test/unit-tests/function/algebra/decomposition/schur.test.js @@ -3,62 +3,145 @@ 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) + assert.ok(Math.abs(T[2][1]) > 0.1, '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/) }) }) From 904c21e061df15a893030fd0bbed4fe151f080a3 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 11 Jan 2026 12:47:13 +0000 Subject: [PATCH 3/6] Address code review feedback - Simplify Householder sigma check condition - Clarify discriminant comment for eigenvalue classification - Fix potential undefined access in exceptional shift calculation - Update docstring example to describe properties rather than specific values - Make 2x2 block test threshold relative to matrix norm Co-authored-by: jacksonloper <2056977+jacksonloper@users.noreply.github.com> --- src/function/algebra/decomposition/schur.js | 22 +++++++++++-------- .../algebra/decomposition/schur.test.js | 4 +++- 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/function/algebra/decomposition/schur.js b/src/function/algebra/decomposition/schur.js index 927883b57d..acabdb761a 100644 --- a/src/function/algebra/decomposition/schur.js +++ b/src/function/algebra/decomposition/schur.js @@ -55,8 +55,11 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( * * Examples: * - * const A = [[1, 0], [-4, 3]] - * math.schur(A) // returns {U: [[0, 1], [-1, 0]], T: [[3, 4], [0, 1]]} + * const A = [[1, 2], [0, 3]] + * const result = math.schur(A) + * // result.U is orthogonal: U * U' = I + * // result.T is upper triangular (quasi-upper-triangular) + * // A = U * T * U' * * See also: * @@ -210,10 +213,7 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( const x0 = x[0] const x0sq = multiplyScalar(x0, x0) - // If the vector is already a multiple of e_1, no transformation needed - if (abs(sigma) < 1e-14 && abs(x0) < 1e-14) { - return null - } + // If the vector is already a multiple of e_1 (sigma ≈ 0), no transformation needed if (abs(sigma) < 1e-14) { return null } @@ -292,7 +292,9 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( const c = H[p][p - 1] const d = H[p][p] - // Discriminant: (a-d)^2 + 4bc (derived from characteristic polynomial) + // 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))) @@ -317,8 +319,10 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( // Apply exceptional shift if convergence is slow if (iterCount === 10 || iterCount === 20) { - // Exceptional shift: use a random-ish perturbation - const shift = multiplyScalar(addScalar(abs(H[p][p - 1]), abs(H[p - 1][p - 2] || 0)), iterCount === 10 ? 1.5 : -1.5) + // 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]) : 0 + 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) } diff --git a/test/unit-tests/function/algebra/decomposition/schur.test.js b/test/unit-tests/function/algebra/decomposition/schur.test.js index d94fa939d5..054a879a6b 100644 --- a/test/unit-tests/function/algebra/decomposition/schur.test.js +++ b/test/unit-tests/function/algebra/decomposition/schur.test.js @@ -116,7 +116,9 @@ describe('schur', function () { 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) - assert.ok(Math.abs(T[2][1]) > 0.1, 'T[2][1] should be non-zero for complex eigenvalue 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 () { From 2d761e6023993dcec2f866f0f337348fd002ae4d Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Sun, 1 Mar 2026 01:56:48 +0100 Subject: [PATCH 4/6] chore: Avoid relying anywhere on matrix entries being JavaScript number Also ensures doc tests are actually run on the schur() function. --- AUTHORS | 1 + src/function/algebra/decomposition/schur.js | 249 ++++++++++---------- src/utils/array.js | 2 +- test/node-tests/doc.test.js | 16 +- 4 files changed, 142 insertions(+), 126 deletions(-) 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 acabdb761a..6bfaf32e66 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' @@ -10,6 +11,10 @@ const dependencies = [ 'norm', 'subtract', 'abs', + 'isZero', + 'isPositive', + 'isNegative', + 'equal', 'addScalar', 'divideScalar', 'multiplyScalar', @@ -28,6 +33,10 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( norm, subtract, abs, + isZero, + isPositive, + isNegative, + equal, addScalar, divideScalar, multiplyScalar, @@ -38,14 +47,16 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( ) => { /** * - * 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 and - * quasi-upper-triangular T such that A = U*T*U'. - * T is block upper triangular with 1x1 and 2x2 blocks on the diagonal. - * 1x1 blocks correspond to real eigenvalues, 2x2 blocks correspond to - * complex conjugate eigenvalue pairs. + * 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 * @@ -55,11 +66,8 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( * * Examples: * - * const A = [[1, 2], [0, 3]] - * const result = math.schur(A) - * // result.U is orthogonal: U * U' = I - * // result.T is upper triangular (quasi-upper-triangular) - * // A = U * T * U' + * const A = [[1, 0], [-4, 3]] + * math.schur(A) // returns {T: [[1, 0], [-4, 3]], U: [[1, 0], [0, 1]]} * * See also: * @@ -69,24 +77,17 @@ 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() - } - }, - - Matrix: function (X) { - return _schur(X) + 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) { - const size = X.size() + function _schur (X, size) { if (size.length !== 2 || size[0] !== size[1]) { throw new RangeError('Matrix must be square') } @@ -94,19 +95,12 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( const n = size[0] // Handle trivial cases - if (n === 0) { - return { U: matrix([]), T: matrix([]) } - } - if (n === 1) { - return { U: identity(1), T: X.clone() } - } - - // Convert to 2D array for internal processing - const arr = X.toArray() + 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(arr, n) + 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 @@ -114,9 +108,7 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( // Step 3: Combine transformations: U = P * Q // So that A = U * T * U' - const U = multiply(matrix(P), matrix(Q)) - - return { U, T: matrix(T) } + return { U: multiply(P, Q), T } } /** @@ -124,15 +116,9 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( * Returns H (upper Hessenberg) and P (orthogonal) such that H = P' * A * P */ function reduceToHessenberg (arr, n) { - // Clone the array to avoid modifying the original - const H = arr.map(row => [...row]) - + const H = clone(arr) // P will accumulate the orthogonal transformation - const P = [] - for (let i = 0; i < n; i++) { - P[i] = Array(n).fill(0) - P[i][i] = 1 - } + const P = identity(n, '') for (let k = 0; k < n - 2; k++) { // Compute Householder vector for column k, rows k+1 to n-1 @@ -141,56 +127,67 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( x.push(H[i][k]) } - const householder = computeHouseholderVector(x) - if (householder === null) { - continue // Column is already zero, skip - } - - const { v, beta } = householder + 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 = 0 + let sum = false + const Hcolj = [] for (let i = 0; i < v.length; i++) { - sum = addScalar(sum, multiplyScalar(v[i], H[k + 1 + i][j])) + 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(H[k + 1 + i][j], multiplyScalar(v[i], sum)) + H[k + 1 + i][j] = subtractScalar(Hcolj[i], multiplyScalar(v[i], sum)) } } // 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++) { - let sum = 0 + const Hrow = H[i] + let sum = false for (let j = 0; j < v.length; j++) { - sum = addScalar(sum, multiplyScalar(H[i][k + 1 + j], v[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(H[i][k + 1 + j], multiplyScalar(sum, v[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++) { - let sum = 0 + const Prow = P[i] + let sum = false for (let j = 0; j < v.length; j++) { - sum = addScalar(sum, multiplyScalar(P[i][k + 1 + j], v[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++) { - P[i][k + 1 + j] = subtractScalar(P[i][k + 1 + j], multiplyScalar(sum, v[j])) + Prow[k + 1 + j] = subtractScalar( + Prow[k + 1 + j], multiplyScalar(sum, v[j])) } } } // 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] = 0 + H[i][j] = zero } } @@ -203,36 +200,35 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( */ function computeHouseholderVector (x) { const m = x.length - if (m === 0) return null + if (m < 2) return { v: null, beta: null } - let sigma = 0 - for (let i = 1; i < m; i++) { + 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 (abs(sigma) < 1e-14) { - return null - } + // 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 - let v0 - if (x0 <= 0) { - v0 = subtractScalar(x0, normX) - } else { - v0 = divideScalar(-sigma, addScalar(x0, normX)) - } + 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 beta = divideScalar(2, addScalar(1, divideScalar(sigma, v0sq))) + 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 = [1] + const v = [one] for (let i = 1; i < m; i++) { v.push(divideScalar(x[i], v0)) } @@ -241,20 +237,16 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( } /** - * Francis QR algorithm with implicit double shift for upper Hessenberg matrices. + * 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 = Hin.map(row => [...row]) - + const H = clone(Hin) // Q accumulates the orthogonal transformations - const Q = [] - for (let i = 0; i < n; i++) { - Q[i] = Array(n).fill(0) - Q[i][i] = 1 - } + const Q = identity(n, '') - const eps = Number.EPSILON // machine epsilon for convergence + const zero = subtractScalar(H[0][0], H[0][0]) const maxIterationsPerEigenvalue = 30 // max iterations per eigenvalue const maxTotalIterations = 30 * n // safety limit @@ -268,9 +260,9 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( // Find the largest q such that H[q][q-1] is negligible while (q > 0) { - const threshold = eps * (abs(H[q - 1][q - 1]) + abs(H[q][q])) - if (abs(H[q][q - 1]) <= threshold) { - H[q][q - 1] = 0 + const scale = addScalar(abs(H[q - 1][q - 1]), abs(H[q][q])) + if (equal(scale, addScalar(scale, H[q][q - 1]))) { + H[q][q - 1] = zero break } q-- @@ -296,9 +288,10 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( // 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))) + const discriminant = addScalar( + multiplyScalar(diff, diff), multiplyScalar(4, multiplyScalar(b, c))) - if (discriminant < 0) { + if (isNegative(discriminant)) { // Complex eigenvalues - keep the 2x2 block p -= 2 iterCount = 0 @@ -321,8 +314,9 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( 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]) : 0 - const shift = multiplyScalar(addScalar(subdiagVal, prevSubdiag), iterCount === 10 ? 1.5 : -1.5) + 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) } @@ -338,9 +332,9 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( // Clean up tiny subdiagonal elements for (let i = 1; i < n; i++) { // Use a relative threshold based on nearby diagonal elements - const threshold = eps * (abs(H[i - 1][i - 1]) + abs(H[i][i])) - if (abs(H[i][i - 1]) <= threshold) { - H[i][i - 1] = 0 + const scale = addScalar(abs(H[i - 1][i - 1]), abs(H[i][i])) + if (equal(scale, addScalar(scale, H[i][i - 1]))) { + H[i][i - 1] = zero } } @@ -359,6 +353,8 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( 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) @@ -369,7 +365,8 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( // = 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) + // 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 @@ -397,7 +394,7 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( ) // M[q+2][q] = H[q+2][q+1]*H[q+1][q] (only nonzero element from Hessenberg structure) - let z = 0 + let z = zero if (q + 2 <= p) { z = multiplyScalar(H[q + 2][q + 1], Hq1q) } @@ -428,40 +425,48 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( 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 = 0 + let sum = false + const Hcolj = [] for (let i = 0; i < r; i++) { - sum = addScalar(sum, multiplyScalar(v[i], H[k + i][j])) + 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(H[k + i][j], multiplyScalar(v[i], sum)) + 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(H[i][k + j], v[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(H[i][k + j], multiplyScalar(sum, v[j])) + H[i][k + j] = subtractScalar(Hrow[k + j], multiplyScalar(sum, v[j])) } } // Accumulate Q for (let i = 0; i < n; i++) { - let sum = 0 - for (let j = 0; j < r; j++) { + 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++) { - Q[i][k + j] = subtractScalar(Q[i][k + j], multiplyScalar(sum, v[j])) + Qrow[k + j] = subtractScalar(Qrow[k + j], multiplyScalar(sum, v[j])) } } @@ -476,6 +481,7 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( /** * 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( @@ -483,17 +489,18 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( multiplyScalar(y, y) ), multiplyScalar(z, z))) - if (abs(norm) < 1e-14) { - return null - } + if (isZero(norm)) return null + const zero = subtractScalar(x, x) // Choose sign to avoid cancellation - const s = x >= 0 ? 1 : -1 - const u0 = addScalar(x, multiplyScalar(s, norm)) + const signedNorm = isNegative(x) ? subtractScalar(zero, norm) : norm + const u0 = addScalar(x, signedNorm) - const v = [1, divideScalar(y, u0), divideScalar(z, u0)] - const vNormSq = addScalar(addScalar(1, multiplyScalar(v[1], v[1])), multiplyScalar(v[2], v[2])) - const beta = divideScalar(2, vNormSq) + 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 } } @@ -504,17 +511,17 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( function computeHouseholderVector2 (x, y) { const norm = sqrt(addScalar(multiplyScalar(x, x), multiplyScalar(y, y))) - if (abs(norm) < 1e-14) { - return null - } + if (isZero(norm)) return null + const zero = subtractScalar(x, x) // Choose sign to avoid cancellation - const s = x >= 0 ? 1 : -1 - const u0 = addScalar(x, multiplyScalar(s, norm)) + const signedNorm = isNegative(x) ? subtractScalar(zero, norm) : norm + const u0 = addScalar(x, signedNorm) - const v = [1, divideScalar(y, u0)] - const vNormSq = addScalar(1, multiplyScalar(v[1], v[1])) - const beta = divideScalar(2, vNormSq) + 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' From 374778674d4522c7bf7236bbff1cf0ca0e5c86a7 Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Sun, 1 Mar 2026 02:03:34 +0100 Subject: [PATCH 5/6] chore: remove unused dependencies --- src/function/algebra/decomposition/schur.js | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/src/function/algebra/decomposition/schur.js b/src/function/algebra/decomposition/schur.js index 6bfaf32e66..eb4d80d5ef 100644 --- a/src/function/algebra/decomposition/schur.js +++ b/src/function/algebra/decomposition/schur.js @@ -7,9 +7,6 @@ const dependencies = [ 'matrix', 'identity', 'multiply', - 'qr', - 'norm', - 'subtract', 'abs', 'isZero', 'isPositive', @@ -19,8 +16,7 @@ const dependencies = [ 'divideScalar', 'multiplyScalar', 'subtractScalar', - 'sqrt', - 'transpose' + 'sqrt' ] export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( @@ -29,9 +25,6 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( matrix, identity, multiply, - qr, - norm, - subtract, abs, isZero, isPositive, @@ -41,8 +34,7 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( divideScalar, multiplyScalar, subtractScalar, - sqrt, - transpose + sqrt } ) => { /** From 9164f99f921efd065fcafa2c3d91f5f17665e038 Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Sun, 1 Mar 2026 08:32:08 +0100 Subject: [PATCH 6/6] refactor: use the simpler equalScalar instead of equal --- src/function/algebra/decomposition/schur.js | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/function/algebra/decomposition/schur.js b/src/function/algebra/decomposition/schur.js index eb4d80d5ef..c69d7c30d2 100644 --- a/src/function/algebra/decomposition/schur.js +++ b/src/function/algebra/decomposition/schur.js @@ -11,7 +11,7 @@ const dependencies = [ 'isZero', 'isPositive', 'isNegative', - 'equal', + 'equalScalar', 'addScalar', 'divideScalar', 'multiplyScalar', @@ -29,7 +29,7 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( isZero, isPositive, isNegative, - equal, + equalScalar, addScalar, divideScalar, multiplyScalar, @@ -253,7 +253,7 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( // 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 (equal(scale, addScalar(scale, H[q][q - 1]))) { + if (equalScalar(scale, addScalar(scale, H[q][q - 1]))) { H[q][q - 1] = zero break } @@ -325,7 +325,7 @@ export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( 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 (equal(scale, addScalar(scale, H[i][i - 1]))) { + if (equalScalar(scale, addScalar(scale, H[i][i - 1]))) { H[i][i - 1] = zero } }