spb/metrika Public

1//2//  LinearAlgebra.swift3//  Metrika4//5//  Author:  Simon-Pierre Boucher6//  Contact: contact@spboucher.ai7//  Copyright © 2026 Simon-Pierre Boucher. All rights reserved.8//910import Accelerate1112/// Thin LAPACK/BLAS wrappers used by the CPU estimators. All matrices are13/// stored column-major (LAPACK convention) as flat [Double] buffers.14enum LinearAlgebra {15    struct Failure: Error, CustomStringConvertible {16        let routine: String17        let info: Int18        var description: String { "\(routine) failed with info=\(info)" }19    }2021    /// QR decomposition of an n×k matrix (n ≥ k) via `dgeqrf`. CLAUDE.md §5:22    /// never form X'X — coefficient solves and (X'X)⁻¹ all go through R.23    struct QR {24        /// Factored matrix as returned by dgeqrf: R in the upper triangle,25        /// Householder vectors below the diagonal.26        private(set) var factored: [Double]27        private(set) var tau: [Double]28        let rows: Int29        let cols: Int3031        init(matrix: [Double], rows: Int, cols: Int) throws {32            precondition(matrix.count == rows * cols, "matrix size mismatch")33            precondition(rows >= cols, "QR requires rows ≥ cols")34            var a = matrix35            var m = __CLPK_integer(rows)36            var n = __CLPK_integer(cols)37            var lda = m38            var tau = [Double](repeating: 0, count: cols)39            var info: __CLPK_integer = 04041            // Workspace query, then factorization.42            var lwork: __CLPK_integer = -143            var workQuery = [Double](repeating: 0, count: 1)44            dgeqrf_(&m, &n, &a, &lda, &tau, &workQuery, &lwork, &info)45            lwork = __CLPK_integer(workQuery[0])46            var work = [Double](repeating: 0, count: max(1, Int(lwork)))47            dgeqrf_(&m, &n, &a, &lda, &tau, &work, &lwork, &info)48            guard info == 0 else { throw Failure(routine: "dgeqrf", info: Int(info)) }4950            self.factored = a51            self.tau = tau52            self.rows = rows53            self.cols = cols54        }5556        /// Solves min ‖Ax − b‖ using the stored factorization:57        /// x = R⁻¹ (Q'b)[0..<k].58        func solve(rhs: [Double]) throws -> [Double] {59            precondition(rhs.count == rows, "rhs size mismatch")60            var qtb = rhs61            var side: Int8 = Int8(UInt8(ascii: "L"))62            var trans: Int8 = Int8(UInt8(ascii: "T"))63            var m = __CLPK_integer(rows)64            var nrhs: __CLPK_integer = 165            var k = __CLPK_integer(cols)66            var lda = m67            var ldc = m68            var info: __CLPK_integer = 069            var a = factored70            var tauCopy = tau7172            var lwork: __CLPK_integer = -173            var workQuery = [Double](repeating: 0, count: 1)74            dormqr_(&side, &trans, &m, &nrhs, &k, &a, &lda, &tauCopy, &qtb, &ldc,75                    &workQuery, &lwork, &info)76            lwork = __CLPK_integer(workQuery[0])77            var work = [Double](repeating: 0, count: max(1, Int(lwork)))78            dormqr_(&side, &trans, &m, &nrhs, &k, &a, &lda, &tauCopy, &qtb, &ldc,79                    &work, &lwork, &info)80            guard info == 0 else { throw Failure(routine: "dormqr", info: Int(info)) }8182            // Back-substitute R x = (Q'b)[0..<k].83            var uplo: Int8 = Int8(UInt8(ascii: "U"))84            var transN: Int8 = Int8(UInt8(ascii: "N"))85            var diag: Int8 = Int8(UInt8(ascii: "N"))86            var n = __CLPK_integer(cols)87            var solution = Array(qtb[0..<cols])88            var ldb = n89            dtrtrs_(&uplo, &transN, &diag, &n, &nrhs, &a, &lda, &solution, &ldb, &info)90            guard info == 0 else { throw Failure(routine: "dtrtrs", info: Int(info)) }91            return solution92        }9394        /// (X'X)⁻¹ = R⁻¹ R⁻ᵀ, computed from the k×k triangular factor.95        func crossProductInverse() throws -> [Double] {96            // Extract the k×k upper triangle of R (column-major).97            var r = [Double](repeating: 0, count: cols * cols)98            for j in 0..<cols {99                for i in 0...j {100                    r[j * cols + i] = factored[j * rows + i]101                }102            }103            var uplo: Int8 = Int8(UInt8(ascii: "U"))104            var diag: Int8 = Int8(UInt8(ascii: "N"))105            var n = __CLPK_integer(cols)106            var lda = n107            var info: __CLPK_integer = 0108            dtrtri_(&uplo, &diag, &n, &r, &lda, &info)109            guard info == 0 else { throw Failure(routine: "dtrtri", info: Int(info)) }110111            // XtXinv = Rinv · Rinvᵀ (symmetric k×k).112            var result = [Double](repeating: 0, count: cols * cols)113            for j in 0..<cols {114                for i in 0..<cols {115                    var sum = 0.0116                    // Rinv is upper triangular: Rinv[i, l] nonzero for l ≥ i.117                    for l in max(i, j)..<cols {118                        sum += r[l * cols + i] * r[l * cols + j]119                    }120                    result[j * cols + i] = sum121                }122            }123            return result124        }125126        /// Explicit thin Q (n×k, column-major) via `dorgqr`, used for127        /// leverage values h_i = ‖q_i‖².128        func thinQ() throws -> [Double] {129            var a = factored130            var m = __CLPK_integer(rows)131            var n = __CLPK_integer(cols)132            var k = n133            var lda = m134            var info: __CLPK_integer = 0135            var tauCopy = tau136137            var lwork: __CLPK_integer = -1138            var workQuery = [Double](repeating: 0, count: 1)139            dorgqr_(&m, &n, &k, &a, &lda, &tauCopy, &workQuery, &lwork, &info)140            lwork = __CLPK_integer(workQuery[0])141            var work = [Double](repeating: 0, count: max(1, Int(lwork)))142            dorgqr_(&m, &n, &k, &a, &lda, &tauCopy, &work, &lwork, &info)143            guard info == 0 else { throw Failure(routine: "dorgqr", info: Int(info)) }144            return a145        }146    }147148    /// y ← A·x for column-major A (rows×cols).149    static func multiply(150        matrix a: [Double], rows: Int, cols: Int, vector x: [Double]151    ) -> [Double] {152        precondition(x.count == cols)153        var result = [Double](repeating: 0, count: rows)154        cblas_dgemv(155            CblasColMajor, CblasNoTrans,156            Int32(rows), Int32(cols),157            1.0, a, Int32(rows), x, 1, 0.0, &result, 1158        )159        return result160    }161162    /// Symmetric sandwich product B · M · B for k×k column-major matrices.163    static func sandwich(bread: [Double], meat: [Double], k: Int) -> [Double] {164        precondition(bread.count == k * k && meat.count == k * k)165        var bm = [Double](repeating: 0, count: k * k)166        cblas_dgemm(167            CblasColMajor, CblasNoTrans, CblasNoTrans,168            Int32(k), Int32(k), Int32(k),169            1.0, bread, Int32(k), meat, Int32(k),170            0.0, &bm, Int32(k)171        )172        var result = [Double](repeating: 0, count: k * k)173        cblas_dgemm(174            CblasColMajor, CblasNoTrans, CblasNoTrans,175            Int32(k), Int32(k), Int32(k),176            1.0, bm, Int32(k), bread, Int32(k),177            0.0, &result, Int32(k)178        )179        return result180    }181182    /// Solves the symmetric positive-definite system A·x = b (k×k) via183    /// Cholesky — used for Wald tests on coefficient subsets.184    static func solveSymmetric(185        _ a: [Double], k: Int, rhs: [Double]186    ) throws -> [Double] {187        var factor = a188        var solution = rhs189        var uplo: Int8 = Int8(UInt8(ascii: "U"))190        var n = __CLPK_integer(k)191        var nrhs: __CLPK_integer = 1192        var lda = n193        var ldb = n194        var info: __CLPK_integer = 0195        dposv_(&uplo, &n, &nrhs, &factor, &lda, &solution, &ldb, &info)196        guard info == 0 else { throw Failure(routine: "dposv", info: Int(info)) }197        return solution198    }199}200