spb/uqo-eval

1// UQO Éval — algèbre linéaire dense (Float64Array, ligne-major) pour l'estimation hédonique.2// Aucune dépendance : Cholesky, résolution, inverse, produits — assez pour k ≤ ~350.34export type Mat = { n: number; k: number; a: Float64Array }; // n lignes × k colonnes, a[i*k+j]56export const zeros = (n: number, k: number): Mat => ({ n, k, a: new Float64Array(n * k) });78/** G = X'X (symétrique k×k) et c = X'y. Coût n·k²/2. */9export function gram(X: Mat, y?: Float64Array, w?: Float64Array): { G: Float64Array; c: Float64Array } {10  const { n, k, a } = X;11  const G = new Float64Array(k * k);12  const c = new Float64Array(k);13  for (let i = 0; i < n; i++) {14    const off = i * k;15    const wi = w ? w[i] : 1;16    for (let j = 0; j < k; j++) {17      const xij = a[off + j] * wi;18      if (xij === 0) continue;19      const rowj = j * k;20      for (let l = j; l < k; l++) G[rowj + l] += xij * a[off + l];21      if (y) c[j] += xij * y[i];22    }23  }24  for (let j = 0; j < k; j++) for (let l = 0; l < j; l++) G[j * k + l] = G[l * k + j];25  return { G, c };26}2728/** Cholesky in-place : A = L·L' (A symétrique définie positive, k×k). Retourne L (triangulaire inf.) ou null. */29export function cholesky(A: Float64Array, k: number, jitter = 0): Float64Array | null {30  const L = new Float64Array(k * k);31  for (let j = 0; j < k; j++) {32    let s = A[j * k + j] + jitter;33    for (let p = 0; p < j; p++) s -= L[j * k + p] * L[j * k + p];34    if (!(s > 1e-12)) return null;35    const ljj = Math.sqrt(s);36    L[j * k + j] = ljj;37    for (let i = j + 1; i < k; i++) {38      let t = A[i * k + j];39      for (let p = 0; p < j; p++) t -= L[i * k + p] * L[j * k + p];40      L[i * k + j] = t / ljj;41    }42  }43  return L;44}4546/** Résout L L' x = b. */47export function cholSolve(L: Float64Array, k: number, b: Float64Array): Float64Array {48  const y = new Float64Array(k);49  for (let i = 0; i < k; i++) {50    let s = b[i];51    for (let p = 0; p < i; p++) s -= L[i * k + p] * y[p];52    y[i] = s / L[i * k + i];53  }54  const x = new Float64Array(k);55  for (let i = k - 1; i >= 0; i--) {56    let s = y[i];57    for (let p = i + 1; p < k; p++) s -= L[p * k + i] * x[p];58    x[i] = s / L[i * k + i];59  }60  return x;61}6263/** Inverse d'une matrice SDP via Cholesky. */64export function cholInverse(L: Float64Array, k: number): Float64Array {65  const inv = new Float64Array(k * k);66  const e = new Float64Array(k);67  for (let j = 0; j < k; j++) {68    e.fill(0);69    e[j] = 1;70    const col = cholSolve(L, k, e);71    for (let i = 0; i < k; i++) inv[i * k + j] = col[i];72  }73  return inv;74}7576/**77 * Résout (A + ridge·I) x = b avec repli : si Cholesky échoue (colinéarité parfaite),78 * on ajoute une régularisation infinitésimale croissante. Retourne aussi l'inverse.79 */80export function solveSPD(A: Float64Array, k: number, b: Float64Array, ridge = 0): { x: Float64Array; inv: Float64Array; jitter: number } {81  let jitter = ridge;82  for (let attempt = 0; attempt < 8; attempt++) {83    const L = cholesky(A, k, jitter);84    if (L) return { x: cholSolve(L, k, b), inv: cholInverse(L, k), jitter };85    // échelle du problème : trace moyenne86    let tr = 0;87    for (let j = 0; j < k; j++) tr += A[j * k + j];88    jitter = Math.max(jitter * 10, (tr / k) * 1e-10, 1e-10);89  }90  throw new Error("matrice X'X singulière — retirez une variable ou un effet fixe colinéaire");91}9293export function matVec(X: Mat, beta: Float64Array): Float64Array {94  const { n, k, a } = X;95  const out = new Float64Array(n);96  for (let i = 0; i < n; i++) {97    let s = 0;98    const off = i * k;99    for (let j = 0; j < k; j++) s += a[off + j] * beta[j];100    out[i] = s;101  }102  return out;103}104105/** Sous-ensemble de lignes. */106export function takeRows(X: Mat, idx: Int32Array | number[]): Mat {107  const k = X.k;108  const out = zeros(idx.length, k);109  for (let r = 0; r < idx.length; r++) out.a.set(X.a.subarray(idx[r] * k, idx[r] * k + k), r * k);110  return out;111}112113export function takeVec(v: Float64Array, idx: Int32Array | number[]): Float64Array {114  const out = new Float64Array(idx.length);115  for (let r = 0; r < idx.length; r++) out[r] = v[idx[r]];116  return out;117}118119/** Colonnes sélectionnées. */120export function takeCols(X: Mat, cols: number[]): Mat {121  const out = zeros(X.n, cols.length);122  for (let i = 0; i < X.n; i++) for (let j = 0; j < cols.length; j++) out.a[i * cols.length + j] = X.a[i * X.k + cols[j]];123  return out;124}125126/** Concatène deux matrices côte à côte. */127export function hcat(A: Mat, B: Mat): Mat {128  const k = A.k + B.k;129  const out = zeros(A.n, k);130  for (let i = 0; i < A.n; i++) {131    out.a.set(A.a.subarray(i * A.k, i * A.k + A.k), i * k);132    out.a.set(B.a.subarray(i * B.k, i * B.k + B.k), i * k + A.k);133  }134  return out;135}136137export const mean = (v: ArrayLike<number>): number => {138  let s = 0;139  for (let i = 0; i < v.length; i++) s += v[i];140  return v.length ? s / v.length : NaN;141};142143export function quantile(sorted: ArrayLike<number>, q: number): number {144  if (!sorted.length) return NaN;145  const pos = (sorted.length - 1) * q;146  const lo = Math.floor(pos);147  const hi = Math.ceil(pos);148  return sorted[lo] + (sorted[hi] - sorted[lo]) * (pos - lo);149}150151export const sortedCopy = (v: ArrayLike<number>): Float64Array => Float64Array.from(v as ArrayLike<number>).sort();152153export const median = (v: ArrayLike<number>): number => quantile(sortedCopy(v), 0.5);154155export function sd(v: ArrayLike<number>): number {156  const m = mean(v);157  let s = 0;158  for (let i = 0; i < v.length; i++) s += (v[i] - m) ** 2;159  return Math.sqrt(s / Math.max(1, v.length - 1));160}161