spb/earth-now Public License

1/**2 * earth-now.co3 * Author:  Simon-Pierre Boucher4 * Contact: contact@spboucher.ai5 * File:    packages/models/src/linalg.ts6 * Purpose: Minimal dense linear algebra — least-squares via normal equations with partial-pivot Gaussian elimination7 */89/** Solve A x = b (A: n×n, row-major) with partial pivoting. Throws on singular systems. */10export function solveLinearSystem(A: number[][], b: number[]): number[] {11  const n = b.length;12  // Augmented copy so callers' arrays are never mutated.13  const M = A.map((row, i) => [...row, b[i]!]);1415  for (let col = 0; col < n; col++) {16    let pivot = col;17    for (let r = col + 1; r < n; r++) {18      if (Math.abs(M[r]![col]!) > Math.abs(M[pivot]![col]!)) pivot = r;19    }20    if (Math.abs(M[pivot]![col]!) < 1e-12) throw new Error("solveLinearSystem: singular matrix");21    if (pivot !== col) {22      const tmp = M[col]!;23      M[col] = M[pivot]!;24      M[pivot] = tmp;25    }26    for (let r = col + 1; r < n; r++) {27      const f = M[r]![col]! / M[col]![col]!;28      for (let c = col; c <= n; c++) M[r]![c]! -= f * M[col]![c]!;29    }30  }3132  const x = new Array<number>(n).fill(0);33  for (let r = n - 1; r >= 0; r--) {34    let acc = M[r]![n]!;35    for (let c = r + 1; c < n; c++) acc -= M[r]![c]! * x[c]!;36    x[r] = acc / M[r]![r]!;37  }38  return x;39}4041/**42 * Ordinary least squares: minimize ||X β − y||² via normal equations XᵀX β = Xᵀy.43 * X is m×k row-major (m observations, k basis functions).44 */45export function leastSquares(X: number[][], y: number[]): number[] {46  const m = X.length;47  if (m === 0 || m !== y.length) throw new Error("leastSquares: dimension mismatch");48  const k = X[0]!.length;49  const XtX: number[][] = Array.from({ length: k }, () => new Array<number>(k).fill(0));50  const Xty = new Array<number>(k).fill(0);51  for (let i = 0; i < m; i++) {52    const row = X[i]!;53    for (let a = 0; a < k; a++) {54      Xty[a]! += row[a]! * y[i]!;55      for (let b = a; b < k; b++) XtX[a]![b]! += row[a]! * row[b]!;56    }57  }58  for (let a = 0; a < k; a++) for (let b = 0; b < a; b++) XtX[a]![b] = XtX[b]![a]!;59  return solveLinearSystem(XtX, Xty);60}61