Advanced Linear Algebra Toolkit

1/**2 * Example 49: Advanced Linear Algebra Toolkit3 *4 * Documents the v1.0.0 linear algebra expansion beyond the earlier SVD/QR/LU5 * walkthrough: Hessenberg reduction, Schur/polar decompositions, matrix6 * functions, structured solvers, sparse CSR solving, and matrix equations.7 */89import {10  block_diag,11  denseToCSR,12  expm,13  hadamard,14  hessenberg,15  logm,16  lyapunov,17  matrix_power,18  polar,19  schur,20  solve_banded,21  sparseSolve,22  sqrtm,23  sylvester,24  toeplitz,25} from "deepbox/linalg";26import { type Tensor, tensor } from "deepbox/ndarray";2728function toMatrix(t: Tensor): number[][] {29  const rows = t.shape[0] ?? 0;30  const cols = t.shape[1] ?? 0;31  return Array.from({ length: rows }, (_, i) =>32    Array.from({ length: cols }, (_, j) => Number(t.at(i, j)))33  );34}3536function matmul(a: number[][], b: number[][]): number[][] {37  const rows = a.length;38  const cols = b[0]?.length ?? 0;39  const inner = b.length;40  return Array.from({ length: rows }, (_, i) =>41    Array.from({ length: cols }, (_, j) => {42      let sum = 0;43      for (let k = 0; k < inner; k++) {44        sum += (a[i]?.[k] ?? 0) * (b[k]?.[j] ?? 0);45      }46      return sum;47    })48  );49}5051function transpose(a: number[][]): number[][] {52  const rows = a.length;53  const cols = a[0]?.length ?? 0;54  return Array.from({ length: cols }, (_, j) =>55    Array.from({ length: rows }, (_, i) => a[i]?.[j] ?? 0)56  );57}5859function frobeniusDiff(a: number[][], b: number[][]): number {60  let sum = 0;61  for (let i = 0; i < a.length; i++) {62    for (let j = 0; j < (a[i]?.length ?? 0); j++) {63      const diff = (a[i]?.[j] ?? 0) - (b[i]?.[j] ?? 0);64      sum += diff * diff;65    }66  }67  return Math.sqrt(sum);68}6970function matvec(a: number[][], x: number[]): number[] {71  return a.map((row) => row.reduce((sum, value, index) => sum + value * (x[index] ?? 0), 0));72}7374console.log("=".repeat(72));75console.log("Example 49: Advanced Linear Algebra Toolkit");76console.log("=".repeat(72));7778// ============================================================================79// Part 1: Hessenberg and Schur decompositions80// ============================================================================81console.log("\n🏗️  Part 1: Hessenberg + Schur");82console.log("-".repeat(72));8384const systemMatrix = tensor([85  [4, 1, 0],86  [1, 3, 1],87  [0, 1, 2],88]);8990const [hessenbergForm, hessenbergQ] = hessenberg(systemMatrix);91const [schurT, schurQ] = schur(systemMatrix);9293const systemDense = toMatrix(systemMatrix);94const hessenbergReconstruction = matmul(95  matmul(toMatrix(hessenbergQ), toMatrix(hessenbergForm)),96  transpose(toMatrix(hessenbergQ))97);98const schurReconstruction = matmul(99  matmul(toMatrix(schurQ), toMatrix(schurT)),100  transpose(toMatrix(schurQ))101);102103console.log(104  `Hessenberg reconstruction error: ${frobeniusDiff(hessenbergReconstruction, systemDense).toExponential(3)}`105);106console.log(107  `Schur reconstruction error:      ${frobeniusDiff(schurReconstruction, systemDense).toExponential(3)}`108);109console.log("Upper Hessenberg form:");110console.log(hessenbergForm.toString());111112// ============================================================================113// Part 2: Polar decomposition114// ============================================================================115console.log("\n🧭 Part 2: Polar Decomposition");116console.log("-".repeat(72));117118const featureTransform = tensor([119  [1.2, 0.3],120  [-0.4, 0.9],121]);122const [orthogonalPart, positivePart] = polar(featureTransform);123const polarReconstruction = matmul(toMatrix(orthogonalPart), toMatrix(positivePart));124125console.log(126  `Polar reconstruction error: ${frobeniusDiff(polarReconstruction, toMatrix(featureTransform)).toExponential(3)}`127);128console.log("Orthogonal factor U:");129console.log(orthogonalPart.toString());130console.log("Positive-semidefinite factor P:");131console.log(positivePart.toString());132133// ============================================================================134// Part 3: Matrix functions and powers135// ============================================================================136console.log("\n🧮 Part 3: Matrix Functions");137console.log("-".repeat(72));138139const diagonalDynamics = tensor([140  [1.1, 0],141  [0, 0.85],142]);143const diagonalExp = expm(diagonalDynamics);144const diagonalLog = logm(diagonalExp);145const covariance = tensor([146  [4, 1],147  [1, 3],148]);149const covarianceSqrt = sqrtm(covariance);150const covarianceRecovered = matmul(toMatrix(covarianceSqrt), toMatrix(covarianceSqrt));151const transition = tensor([152  [0.92, 0.08],153  [0.05, 0.95],154]);155const fiveStepTransition = matrix_power(transition, 5);156157console.log("expm(A) for diagonal dynamics:");158console.log(diagonalExp.toString());159console.log("logm(expm(A)) recovers:");160console.log(diagonalLog.toString());161console.log(162  `sqrtm(C) * sqrtm(C) reconstruction error: ${frobeniusDiff(covarianceRecovered, toMatrix(covariance)).toExponential(3)}`163);164console.log("Five-step Markov transition:");165console.log(fiveStepTransition.toString());166167// ============================================================================168// Part 4: Structured dense and sparse solvers169// ============================================================================170console.log("\n🪜 Part 4: Structured Solvers");171console.log("-".repeat(72));172173const tridiagonalBands = tensor([174  [0, -1, -1, -1],175  [4, 4, 4, 4],176  [-1, -1, -1, 0],177]);178const tridiagonalRhs = tensor([15, 10, 10, 15]);179const tridiagonalSolution = solve_banded([1, 1], tridiagonalBands, tridiagonalRhs);180181const sparseSystem = tensor([182  [4, -1, 0, 0],183  [-1, 4, -1, 0],184  [0, -1, 4, -1],185  [0, 0, -1, 3],186]);187const sparseRhs = tensor([15, 10, 10, 10]);188const sparseSolution = sparseSolve(denseToCSR(sparseSystem), sparseRhs);189190console.log("Banded tridiagonal solution:");191console.log(tridiagonalSolution.toString());192console.log("Sparse CSR solution:");193console.log(sparseSolution.toString());194195const denseSparseSystem = toMatrix(sparseSystem);196const sparseResidual = matvec(197  denseSparseSystem,198  Array.from({ length: sparseSolution.shape[0] ?? 0 }, (_, i) => Number(sparseSolution.at(i)))199);200console.log(201  `Sparse residual preview: ${sparseResidual.map((value, index) => (value - Number(sparseRhs.at(index))).toFixed(6)).join(", ")}`202);203204// ============================================================================205// Part 5: Sylvester and Lyapunov equations206// ============================================================================207console.log("\n📐 Part 5: Matrix Equations");208console.log("-".repeat(72));209210const aSylvester = tensor([211  [1, 2],212  [0, 3],213]);214const bSylvester = tensor([215  [4, 1],216  [0, 5],217]);218const cSylvester = tensor([219  [10, 15],220  [12, 24],221]);222const sylvesterSolution = sylvester(aSylvester, bSylvester, cSylvester);223224const stableA = tensor([225  [-1, 0.2],226  [0, -0.6],227]);228const lyapunovQ = tensor([229  [2, 0.5],230  [0.5, 1],231]);232const lyapunovSolution = lyapunov(stableA, lyapunovQ);233234console.log("Sylvester solution X:");235console.log(sylvesterSolution.toString());236console.log("Lyapunov solution X:");237console.log(lyapunovSolution.toString());238239// ============================================================================240// Part 6: Special matrix builders241// ============================================================================242console.log("\n🧱 Part 6: Special Matrices");243console.log("-".repeat(72));244245const smoothingToeplitz = toeplitz([4, -1, 0]);246const hadamardBasis = hadamard(4);247const blockSystem = block_diag(smoothingToeplitz, tensor([[2]]));248249console.log("Toeplitz smoothing kernel:");250console.log(smoothingToeplitz.toString());251console.log("Hadamard basis (n=4):");252console.log(hadamardBasis.toString());253console.log(`Block-diagonal system shape: [${blockSystem.shape.join(", ")}]`);254255// ============================================================================256// Summary257// ============================================================================258console.log("\n💡 Key Takeaways");259console.log("-".repeat(72));260console.log(261  "• Hessenberg and Schur factorizations are the backbone for many advanced eigensolver workflows."262);263console.log("• Polar decomposition splits a transform into rotation-like and scale-like factors.");264console.log(265  "• expm/logm/sqrtm/matrix_power let you move between discrete and continuous matrix dynamics."266);267console.log(268  "• solve_banded and sparseSolve are the practical path once dense solves become structured or sparse."269);270console.log(271  "• Sylvester and Lyapunov solvers are useful for control, filtering, and state-space tooling."272);273console.log(274  "• Special matrix constructors help build repeatable test systems and numerical examples quickly."275);276277console.log("\n✅ Advanced Linear Algebra Toolkit Complete!");278console.log("=".repeat(72));279