spb/metrika Public

1//2//  CholeskySolver.swift3//  Metrika4//5//  Author:  Simon-Pierre Boucher6//  Contact: contact@spboucher.ai7//  Copyright © 2026 Simon-Pierre Boucher. All rights reserved.8//910/// Plain-Swift Cholesky solve for the k×k normal equations of the GPU11/// bootstrap. k is the regressor count (single digits in practice), so a12/// dependency-free O(k³) routine beats reaching for LAPACK — and ZQStats13/// is the only module allowed to import Accelerate (CLAUDE.md §3).14enum CholeskySolver {15    struct SingularMatrix: Error {}1617    /// Solves A·x = b for symmetric positive-definite A (column-major,18    /// k×k). Throws on non-positive-definite input.19    static func solve(_ a: [Double], k: Int, rhs: [Double]) throws -> [Double] {20        precondition(a.count == k * k && rhs.count == k)2122        // Lower-triangular factor L with A = L·Lᵀ.23        var l = [Double](repeating: 0, count: k * k)24        for j in 0..<k {25            var diagonal = a[j * k + j]26            for p in 0..<j {27                diagonal -= l[p * k + j] * l[p * k + j]28            }29            guard diagonal > 0 else { throw SingularMatrix() }30            let root = diagonal.squareRoot()31            l[j * k + j] = root32            for i in (j + 1)..<k {33                var value = a[j * k + i]34                for p in 0..<j {35                    value -= l[p * k + i] * l[p * k + j]36                }37                l[j * k + i] = value / root38            }39        }4041        // Forward substitution L·z = b.42        var z = [Double](repeating: 0, count: k)43        for i in 0..<k {44            var value = rhs[i]45            for p in 0..<i {46                value -= l[p * k + i] * z[p]47            }48            z[i] = value / l[i * k + i]49        }5051        // Back substitution Lᵀ·x = z.52        var x = [Double](repeating: 0, count: k)53        for i in stride(from: k - 1, through: 0, by: -1) {54            var value = z[i]55            for p in (i + 1)..<k {56                value -= l[i * k + p] * x[p]57            }58            x[i] = value / l[i * k + i]59        }60        return x61    }62}63