linalg

package
v1.0.0-alpha.13 Latest Latest
Warning

This package is not in the latest version of its module.

Go to latest
Published: Oct 15, 2026 License: MIT Imports: 5 Imported by: 0

Documentation

Overview

Package linalg is Pulse's one linear-algebra core: Pulse-owned dense (Matrix), symmetric (Sym) and vector (Vec) types plus the kernels the engine and the synthetic-data generator factor, solve and invert with.

Two backends

The REFERENCE kernels — Cholesky, CholeskyRidge, SolveSPD, InverseSPD, the general LU solve (FactorLU, LUFactor, Solve) and the CoMoment accumulator's Add and Merge — are pure Go and FMA-free: every product that feeds an addition or subtraction is written as an explicit float64(a*b) conversion, which the Go specification forbids the compiler from fusing into a fused multiply-add. Their output is therefore bit-identical on every architecture, and a caller with a bit-level contract (same seed, same bytes) can rely on it.

Routines with no bit contract — SymEigen, SVD, QR, Rank and ConditionNumber, plus the gonum-backed SPD path FactorSPD (an SPDFactor solved against and inverted) and the dense product Mul — are backed by gonum and may differ in the last ulps across architectures. On one machine the SPD path and Mul equal the raw gonum calls bit for bit (they are thin wrappers, with no policy layered on top), which is what lets the regression engine route through them with byte-identical output. No gonum type ever appears on this package's exported surface: gonum is pre-1.0, and Pulse's public API is frozen.

Order, sign and rank policy

gonum's raw output is canonicalised so every decision above rounding is the same on every architecture:

  • Order: eigenvalues and singular values are descending. Values within the rank tolerance of their neighbour form a tie cluster, ordered by original variable order — the index of each vector's dominant component, lowest first.
  • Sign: every eigenvector, and every right singular vector, is flipped so its largest-magnitude (dominant) component is positive; a magnitude tie (within DominanceTolerance, relative) goes to the lowest index. A singular pair flips U and V together. QR flips so R's diagonal is non-negative.
  • Rank: one tolerance, RankTolerance = max(rows, cols)·Epsilon·σ_max; a singular value at or below it is zero. Rank uses it unless the caller passes an explicit tol > 0, and ConditionNumber reports +Inf for a matrix that is rank-deficient under it.

Reference Cholesky

Cholesky is the lower-triangular Banachiewicz (row-by-row) factor S = L·Lᵀ. For row i and column j ≤ i it starts from S[i][j], subtracts L[i][k]·L[j][k] for k ascending one term at a time, and then either takes the square root (diagonal) or divides by L[j][j] with a true division (off-diagonal). It fails iff a diagonal pivot is ≤ 0. That operation order is a contract: changing it changes bits.

CholeskyRidge retries a failed factorisation with a growing diagonal ridge (DefaultRidgeSchedule: eight tries, increments 1e-6 … 1e1) and returns the ACCUMULATED ridge, so the figure reads as "how far from positive definite was this", not as an iteration count.

General solve

FactorLU factors a general square matrix as P·A = L·U by Doolittle elimination with partial pivoting (largest pivot per column, ties to the lowest row); Solve, LUFactor.Solve, SolveMatrix and Inverse solve against it. They are reference kernels — pure Go, FMA-free — so the factor and every solution are bit-identical on every architecture. A matrix that is rank-deficient under RankTolerance (ConditionNumber +Inf) is refused with PULSE_MATRIX_SINGULAR, never solved to NaN; the refusal is the one decision that rides the gonum singular values.

Co-moments

CoMoment accumulates weighted means and covariances of p variables, Listwise (complete rows) or Pairwise (per-pair complete rows), and merges exactly (Chan–Golub–LeVeque), so partial accumulators built by workers combine to the serial figures. Weights follow Pulse's weighting contract: NaN, ±Inf and negative weights are skipped and counted (NWeightInvalid); a zero weight counts toward N but carries no mass.

For a result that is bit-identical however many workers produced it, fold rows into one CoMoment per MergeBlockSize-record block keyed by ABSOLUTE record index and combine the block sequence with MergeTree, a fixed binary tree by block index: the result is a pure function of the block sequence, not of the segmentation.

Errors

Every failure is a *errors.CodedError: PULSE_MATRIX_SHAPE_MISMATCH for operands whose dimensions disagree (or a nil operand), and PULSE_MATRIX_SINGULAR for a matrix the kernel cannot factor or invert. A gonum-backed routine handed a NaN or infinite element, or whose iteration does not converge, also reports PULSE_MATRIX_SINGULAR, with details "reason" = "non_finite" / "no_convergence".

Imports

linalg imports only the standard library, gonum and Pulse's errors package (gated by TestLinalgImportBoundary).

Index

Constants

View Source
const ConditionTolerance = 1e16

ConditionTolerance is the condition-number ceiling of the gonum-backed SPD path: a factor whose estimated condition number exceeds it still factors, but its Solve and Inverse report PULSE_MATRIX_SINGULAR because their results carry no reliable digit. It equals gonum's own threshold.

View Source
const DominanceTolerance = 0x1p-26

DominanceTolerance is the relative slack used to decide which component of a vector is its LARGEST-MAGNITUDE (dominant) one: every component whose magnitude is at least (1 - DominanceTolerance) times the largest magnitude is tied for dominance, and the tie goes to the lowest index. The slack (2⁻²⁶ = √Epsilon) keeps a decision between two components that are equal in exact arithmetic — 1/√2 and -1/√2, say — from flipping on a last-ulp difference between architectures.

View Source
const Epsilon = 0x1p-52

Epsilon is the float64 machine epsilon, 2⁻⁵², the one constant behind every rank decision in this package.

View Source
const MergeBlockSize = 4096

MergeBlockSize is the blocked-merge block length in records. A block-keyed reducer folds the row at ABSOLUTE record index r into block r / MergeBlockSize — never into "this worker's partial" — so the block sequence, and therefore MergeTree's result, does not depend on how the records were split across workers. The engine cuts its parallel-decode segments at multiples of MergeBlockSize so no block straddles two workers. The value is part of the bit contract: changing it changes merged bits.

Variables

This section is empty.

Functions

func ConditionNumber

func ConditionNumber(m *Matrix) (float64, error)

ConditionNumber returns the 2-norm condition number σ_max / σ_min of m over its min(rows, cols) singular values. A matrix that is rank-deficient under the default tolerance (Rank(m, 0) < min(rows, cols)) — the zero matrix included — has condition number +Inf; a near-singular matrix above the tolerance returns its large finite ratio. A nil or empty m is PULSE_MATRIX_SHAPE_MISMATCH (the condition number of an empty matrix is undefined); a non-finite element is PULSE_MATRIX_SINGULAR (details "reason").

func Rank

func Rank(m *Matrix, tol float64) (int, error)

Rank returns the numerical rank of m: the number of singular values strictly greater than tol. A tol that is not > 0 (zero, negative or NaN) selects the ONE documented default, RankTolerance(rows, cols, σ_max). A nil m is PULSE_MATRIX_SHAPE_MISMATCH; a non-finite element is PULSE_MATRIX_SINGULAR (details "reason"). An empty m has rank 0.

func RankTolerance

func RankTolerance(rows, cols int, sigmaMax float64) float64

RankTolerance is the ONE rank tolerance Pulse uses:

max(rows, cols) · Epsilon · sigmaMax

where sigmaMax is the largest singular value (for a symmetric positive-semidefinite matrix, its largest eigenvalue). A singular value at or below the tolerance counts as zero. This is the LAPACK / NumPy matrix_rank default, so a rank Pulse reports agrees with what an analyst reproduces elsewhere.

Types

type CoMoment

type CoMoment struct {
	// contains filtered or unexported fields
}

CoMoment is a weighted, mergeable accumulator of the means and the covariance (co-moment) matrix of p variables — the one streaming first/second-moment primitive the engine, synth and embedders share.

Update

Add folds one row by the weighted Welford–West recurrence. Merge folds another accumulator by the exact Chan–Golub–LeVeque combine,

W = Wa + Wb,  δ = mean_b − mean_a,
mean = mean_a + δ·Wb/W,
M2 = M2a + M2b + δδᵀ·Wa·Wb/W,

per pair in Pairwise mode. Both are REFERENCE kernels: pure Go and FMA-free (every product feeding an addition is an explicit float64(a*b) conversion), so the same sequence of Add and Merge calls gives the same bits on every architecture. A fixed merge tree over fixed blocks is therefore reproducible bit for bit, whatever runs it.

The first row to carry mass sets the mean to x exactly, and a later row equal to the mean adds exactly zero, so a constant column's variance is exactly 0. With unit weights the recurrence is the classic unweighted Welford one.

Missing values and weights

A NaN in x is a missing value: Listwise skips the row, Pairwise skips it per pair. A row the mode admits nowhere (Listwise: any NaN; Pairwise: every value NaN) is not counted and its weight is never judged. Otherwise the weight is validated as Pulse's weighting contract does: a NaN, ±Inf or negative w skips the row and increments NWeightInvalid; w = 0 is valid — the row counts toward N (and the pair counts) but adds no mass to W, the means or the co-moments; a positive finite w is used as given.

Tolerance

Floating-point addition is not associative, so merging the pieces of a split is not bit-identical to one serial pass. It agrees within 1e-10 relative to the figure's scale: |Δmean_i| ≤ 1e-10·(|mean_i| + sd_i + 1), |ΔCov_ij| ≤ 1e-10·sd_i·sd_j, |ΔCorr_ij| ≤ 1e-10, the weight sums (W, PairW, PairW2) and the Kish sizes (NEff, PairNEff) to 1e-10 relative. Counts (N, NWeightInvalid, PairN, PairNMass) are exact. Merge is associative to the same tolerance, and merging with an accumulator that carries no mass (empty, or zero-weight rows only) leaves every moment bit-identical on either side.

A CoMoment is not safe for concurrent use; give each worker its own and Merge them.

func MergeTree

func MergeTree(blocks []*CoMoment) (*CoMoment, error)

MergeTree folds a sequence of per-block accumulators, ordered by block index, in a FIXED binary tree and returns the root; no input is mutated, and a single block comes back as an independent copy.

The tree is level-wise pairwise by position: blocks (0,1), (2,3), … merge left.Merge(right); an odd tail is carried up to the next level unmerged; repeat until one remains. Five blocks are therefore ((b0·b1)·(b2·b3))·b4. The shape is a pure function of the sequence length, and Merge is an FMA-free reference kernel, so the result is a pure function of the block sequence — bit-identical for any grouping of the same blocks into worker segments, on every architecture.

Every slot must hold an accumulator (an empty one is fine — it is a valid block that saw no rows). An empty sequence, a nil block, or blocks differing in P or Mode is PULSE_MATRIX_SHAPE_MISMATCH.

func NewCoMoment

func NewCoMoment(p int, mode CoMomentMode) (*CoMoment, error)

NewCoMoment returns an empty accumulator over p variables. A negative p or an unknown mode is PULSE_MATRIX_SHAPE_MISMATCH.

func NewCoMomentFromListwise

func NewCoMomentFromListwise(s ListwiseState) (*CoMoment, error)

NewCoMomentFromListwise rebuilds the Listwise accumulator s was captured from: every output, and every later Add or Merge, is bit-identical to the original's. A nil M2, a Mean whose length is not M2's order, a negative count, NMass above N, or a W / W2 that is not finite and non-negative is PULSE_MATRIX_SHAPE_MISMATCH.

func (*CoMoment) Add

func (c *CoMoment) Add(x []float64, w float64)

Add folds one row. len(x) must equal P(); anything else is a programming error and panics, like an out-of-range index. See the type comment for missing values and weights.

func (*CoMoment) Clone

func (c *CoMoment) Clone() *CoMoment

Clone returns an independent deep copy.

func (*CoMoment) Corr

func (c *CoMoment) Corr() *Sym

Corr returns the correlation matrix C_ij/√(M2_ii·M2_jj) — each pair's own co-moments in Pairwise mode — clamped to [−1, 1] against rounding. The one-root form is TEST_PEARSON_R's arithmetic, so a single block of unit-weight rows correlates to the same bits as that test; when the product M2_ii·M2_jj leaves the normal range (overflows to +Inf or falls below the smallest normal float64) the entry takes the two-root form C_ij/(√M2_ii·√M2_jj) instead, so no spread overflows or underflows. The diagonal is exactly 1 for a variable with spread. A variable with zero spread (a constant column, a single massed row, no rows) has NO defined correlation: every entry touching it, its own diagonal included, is NaN — never 0.

func (*CoMoment) Cov

func (c *CoMoment) Cov(ddof int) *Sym

Cov returns the covariance matrix M2/(W − ddof) — per pair in Pairwise mode. ddof = 0 is the weighted population covariance; ddof = 1 the frequency-weight sample covariance (the unweighted sample covariance under unit weights). An entry whose W − ddof ≤ 0 is NaN. Probability-weight corrections are the caller's: NEff is reported beside it.

func (*CoMoment) ListwiseState

func (c *CoMoment) ListwiseState() (ListwiseState, error)

ListwiseState captures c's stored state (independent copies). A Pairwise accumulator has no single state and is PULSE_MATRIX_SHAPE_MISMATCH.

func (*CoMoment) Mean

func (c *CoMoment) Mean() *Vec

Mean returns the weighted means (Pairwise: variable i over the rows where x[i] is present). A variable with no mass has a NaN mean.

func (*CoMoment) Merge

func (c *CoMoment) Merge(other *CoMoment) error

Merge folds other into the receiver (the receiver is MUTATED; other is not) by the exact Chan–Golub–LeVeque combine — so a merge tree is a sequence of left.Merge(right) calls. Counts and NWeightInvalid add. A nil other, or one whose P or Mode differs, is PULSE_MATRIX_SHAPE_MISMATCH and leaves the receiver untouched. c.Merge(c) is allowed and doubles the mass: every element of other is read before the receiver's copy of it is written.

func (*CoMoment) Mode

func (c *CoMoment) Mode() CoMomentMode

Mode returns the missing-value mode.

func (*CoMoment) N

func (c *CoMoment) N() int64

N returns the rows counted: Listwise, the complete rows with a valid weight (zero included); Pairwise, the rows with at least one present value and a valid weight. Rows skipped for an invalid weight are NWeightInvalid, never N.

func (*CoMoment) NEff

func (c *CoMoment) NEff() float64

NEff returns Kish's effective sample size (Σw)²/Σw² over the rows N counts — N itself under unit weights — and 0 when no row carries mass.

func (*CoMoment) NWeightInvalid

func (c *CoMoment) NWeightInvalid() int64

NWeightInvalid returns the rows skipped for a NaN, ±Inf or negative weight (among rows the mode would otherwise have admitted).

func (*CoMoment) P

func (c *CoMoment) P() int

P returns the number of variables.

func (*CoMoment) PairM2

func (c *CoMoment) PairM2(i, j int) float64

PairM2 returns the raw weighted co-moment Σw(xi−mi)(xj−mj) of the pair (i, j) over the rows the pair admits — the M2 numerator Cov divides: Cov(ddof) at (i, j) is exactly PairM2(i, j)/(PairW(i, j) − ddof) when that denominator is positive. The diagonal PairM2(i, i) is variable i's own Σw(xi−mi)² (Pairwise: over the rows where x[i] is present). It is symmetric, 0 for a pair with no mass, and as merge- and worker-stable as Cov (the same stored figure). It panics when i or j is out of range.

func (*CoMoment) PairN

func (c *CoMoment) PairN(i, j int) int64

PairN returns the rows counted for the pair (i, j) — N for every pair in Listwise mode. It panics when i or j is out of range.

func (*CoMoment) PairNEff

func (c *CoMoment) PairNEff(i, j int) float64

PairNEff returns Kish's effective sample size (Σw)²/Σw² for the pair (i, j), over the rows the pair admits — NEff for every pair in Listwise mode — and 0 when the pair carries no mass. It is exactly PairW(i, j)²/PairW2(i, j), the arithmetic the probability-weight inference basis uses for N*. It panics when i or j is out of range.

func (*CoMoment) PairNMass

func (c *CoMoment) PairNMass(i, j int) int64

PairNMass returns the rows counted for the pair (i, j) that carry mass (w > 0) — PairN less the pair's zero-weight rows, and PairN itself under unit weights. A weighted test that skips a zero-weight row (TEST_PEARSON_R's raw n) reads this count. In Listwise mode it is the global count of complete rows with w > 0 for every pair. It panics when i or j is out of range.

func (*CoMoment) PairW

func (c *CoMoment) PairW(i, j int) float64

PairW returns Σw for the pair (i, j) — W for every pair in Listwise mode. It panics when i or j is out of range.

func (*CoMoment) PairW2

func (c *CoMoment) PairW2(i, j int) float64

PairW2 returns Σw² for the pair (i, j) — the squared-weight sum over the rows the pair admits, accumulated (and merged) beside PairW, so the pair's Kish n_eff is PairW²/PairW2 (see PairNEff). In Listwise mode it is the global Σw² for every pair. Zero-weight rows add nothing; it is 0 when the pair carries no mass. It panics when i or j is out of range.

func (*CoMoment) W

func (c *CoMoment) W() float64

W returns Σw over the rows N counts.

type CoMomentMode

type CoMomentMode int

CoMomentMode selects how a CoMoment treats a row with a missing value (a NaN in x).

const (
	// Listwise (the zero value, so the default) drops a row with ANY
	// missing value; every figure is over the same complete rows.
	Listwise CoMomentMode = iota
	// Pairwise keeps, for every pair (i, j), its own count, weight sum,
	// means and co-moments over the rows where BOTH x[i] and x[j] are
	// present. Variable i's mean and variance read the pair (i, i).
	Pairwise
)

func (CoMomentMode) String

func (m CoMomentMode) String() string

String returns "listwise", "pairwise" or "unknown".

type LUFactor

type LUFactor struct {
	// contains filtered or unexported fields
}

LUFactor is the reference LU factorisation P·A = L·U of a general (not necessarily symmetric or positive-definite) square matrix, built by FactorLU with partial pivoting and reused across any number of Solve, SolveMatrix and Inverse calls. L is unit lower triangular, U upper triangular, P the row permutation chosen by the pivoting. It owns its storage: mutating the Matrix it was built from does not affect it.

The factor and every solve against it are a REFERENCE kernel: pure Go, FMA-free, bit-identical on every architecture (see FactorLU for the operation order).

func FactorLU

func FactorLU(a *Matrix) (*LUFactor, error)

FactorLU returns the reference LU factorisation of the square matrix a, P·a = L·U, by Doolittle elimination with partial pivoting.

Operation order (the bit contract): for each column k ascending, the pivot is the row i ≥ k with the largest |a[i][k]| (a magnitude tie goes to the lowest row); that row is swapped into row k; then for every row i > k ascending the multiplier l = a[i][k] / a[k][k] is a true division, and a[i][j] -= l·a[k][j] runs over j > k ascending, one FMA-free term at a time.

Singularity policy — the ONE rank tolerance: a is refused when it is rank-deficient under RankTolerance(n, n, σ_max), i.e. exactly when ConditionNumber(a) is +Inf, and (defensively) when a pivot is exactly zero. The decision rides the gonum-backed singular values, so — like every Rank decision in this package — it is portable except within an ulp of the tolerance itself; the factor and solution VALUES are FMA-free and bit-identical everywhere. A full-rank but ill-conditioned matrix still factors: a solve against it loses roughly log10(κ) digits, and LUFactor.ConditionNumber reports κ for a caller that wants a stricter reliability gate.

Errors: a nil a, or a non-square a, is PULSE_MATRIX_SHAPE_MISMATCH. A NaN or infinite element is PULSE_MATRIX_SINGULAR with details "reason" = "non_finite". A rank-deficient a is PULSE_MATRIX_SINGULAR with details "reason" = "singular", "n", "rank" (Rank(a, 0)) and "condition_number" (+Inf, null on the wire); a zero pivot adds "pivot" (the failing column). A refused matrix returns an error, never a NaN-bearing factor or solution. An empty a factors to an empty factor.

func (*LUFactor) ConditionNumber

func (f *LUFactor) ConditionNumber() float64

ConditionNumber returns the 2-norm condition number σ_max / σ_min of the factored matrix, as ConditionNumber computed it at factor time (gonum-backed, so it may differ in the last ulps across architectures; 1 for an empty factor). It is always finite: a matrix whose condition number is +Inf is refused by FactorLU.

func (*LUFactor) Inverse

func (f *LUFactor) Inverse() (*Matrix, error)

Inverse returns A⁻¹ for the factored A, column j being the solve against the j-th unit vector (SolveMatrix on the identity).

func (*LUFactor) N

func (f *LUFactor) N() int

N returns the order of the factored matrix.

func (*LUFactor) Solve

func (f *LUFactor) Solve(b *Vec) (*Vec, error)

Solve returns x with A·x = b for the factored A: forward substitution L·y = P·b, then back substitution U·x = y, each a running subtraction over k ascending (back substitution: k ascending from i+1) and, for U, a true division by the diagonal. FMA-free. A nil b or len(b) != N() is PULSE_MATRIX_SHAPE_MISMATCH.

func (*LUFactor) SolveMatrix

func (f *LUFactor) SolveMatrix(b *Matrix) (*Matrix, error)

SolveMatrix returns X with A·X = B for the factored A, solving each column of B exactly as Solve does (so column j of X equals Solve on column j of B, bit for bit). A nil B or B.Rows() != N() is PULSE_MATRIX_SHAPE_MISMATCH; a B with zero columns yields an N()×0 matrix.

type ListwiseState

type ListwiseState struct {
	// N, NMass and NWeightInvalid are CoMoment.N, the rows counted with
	// mass (w > 0) and CoMoment.NWeightInvalid.
	N, NMass, NWeightInvalid int64
	// W and W2 are Σw and Σw² over the rows N counts.
	W, W2 float64
	// Mean is the stored running mean, axis order — the raw figure, so
	// an accumulator with no mass keeps its zeros (Mean() reports NaN
	// there).
	Mean []float64
	// M2 is the raw weighted co-moment Σw(xi−mi)(xj−mj) — PairM2 at
	// every (i, j).
	M2 *Sym
}

ListwiseState is the complete stored state of a Listwise CoMoment: every figure its outputs (Mean, Cov, Corr, NEff, …) derive from, at full precision. A fitted state is captured once (CoMoment.ListwiseState) and rebuilt later (NewCoMomentFromListwise) — on another request, or after a round trip through storage — so a consumer that scores against a prebuilt fit derives every block (covariance, correlation, a decomposition) by the same arithmetic, to the same bits, as one that folded the rows itself.

type Matrix

type Matrix struct {
	// contains filtered or unexported fields
}

Matrix is a dense rows×cols matrix of float64 in row-major order. The zero value is an empty 0×0 matrix.

func Cholesky

func Cholesky(s *Sym) (*Matrix, error)

Cholesky returns the lower-triangular factor L of the symmetric positive-definite s, with s = L·Lᵀ and zeros above the diagonal.

It is a REFERENCE kernel: pure Go, FMA-free, bit-identical on every architecture. The operation order is the contract — Banachiewicz row order (i ascending, then j = 0..i), a running subtraction sum -= L[i][k]·L[j][k] over k ascending, a square root on the diagonal and a true division sum / L[j][j] off it. Only the lower triangle of s is read.

A diagonal pivot ≤ 0 is PULSE_MATRIX_SINGULAR (details "pivot", the failing row, plus the singular diagnostics — "rank", "condition_number" and, where identifiable, "dependent_indices"; see singularDiagnostics). A NaN pivot is not ≤ 0 and propagates into L, exactly as the comparison reads. A nil s is PULSE_MATRIX_SHAPE_MISMATCH; an empty s factors to an empty L.

The diagnostics run only on the failure path, through the gonum-backed SymEigen; they never touch L, so the factor's bit contract is unchanged.

func CholeskyRidge

func CholeskyRidge(s *Sym, schedule RidgeSchedule) (*Matrix, float64, error)

CholeskyRidge factors s like Cholesky, retrying a failed attempt with the next schedule increment added to every diagonal element (a nil or empty schedule means DefaultRidgeSchedule). It returns the factor and the ACCUMULATED ridge — the sum of the increments added before the successful attempt, 0 when the first attempt succeeds. s itself is never modified.

When every attempt fails the error is PULSE_MATRIX_SINGULAR (details "attempts" and "ridge") and the returned ridge is the total of the whole schedule, including the increment added after the last attempt.

func Mul

func Mul(a, b *Matrix) (*Matrix, error)

Mul returns the dense product a·b, gonum-backed (no cross-architecture bit contract). A nil operand or a.Cols() != b.Rows() is PULSE_MATRIX_SHAPE_MISMATCH; a zero dimension yields the zero matrix of shape a.Rows()×b.Cols().

func NewMatrix

func NewMatrix(rows, cols int, data []float64) (*Matrix, error)

NewMatrix returns a rows×cols matrix over a COPY of data, which must hold exactly rows*cols values in row-major order. A nil data allocates a zero matrix. Negative dimensions or a length mismatch are PULSE_MATRIX_SHAPE_MISMATCH.

func NewMatrixFromRows

func NewMatrixFromRows(rows [][]float64) (*Matrix, error)

NewMatrixFromRows returns a matrix over a copy of a rectangular [][]float64. A ragged input is PULSE_MATRIX_SHAPE_MISMATCH; an empty input is a 0×0 matrix.

func (*Matrix) At

func (m *Matrix) At(i, j int) float64

At returns element (i, j). It panics when (i, j) is out of range.

func (*Matrix) Cols

func (m *Matrix) Cols() int

Cols returns the column count.

func (*Matrix) RowMajor

func (m *Matrix) RowMajor() []float64

RowMajor returns a copy of the elements in row-major order.

func (*Matrix) Rows

func (m *Matrix) Rows() int

Rows returns the row count.

func (*Matrix) Set

func (m *Matrix) Set(i, j int, v float64)

Set assigns element (i, j). It panics when (i, j) is out of range.

func (*Matrix) ToRows

func (m *Matrix) ToRows() [][]float64

ToRows returns a copy of the matrix as a [][]float64.

type QRResult

type QRResult struct {
	Q *Matrix
	R *Matrix
}

QRResult is the THIN QR decomposition m = Q·R of a rows×cols matrix with rows ≥ cols: Q is rows×cols with orthonormal columns and R is cols×cols upper triangular.

Sign: each diagonal element of R is non-negative (a negative one has its R row and Q column flipped together), which makes the decomposition of a full-column-rank matrix unique.

func QR

func QR(m *Matrix) (*QRResult, error)

QR returns the thin QR decomposition of m under the sign policy documented on QRResult. A nil m or rows < cols is PULSE_MATRIX_SHAPE_MISMATCH; a non-finite element is PULSE_MATRIX_SINGULAR (details "reason"). A matrix with zero columns decomposes to a rows×0 Q and a 0×0 R.

type RidgeSchedule

type RidgeSchedule []float64

RidgeSchedule is the sequence of diagonal increments CholeskyRidge adds between attempts: attempt t runs with the first t increments already added, so a schedule of length m makes m attempts. Increments are the caller's responsibility; they are added as given.

func DefaultRidgeSchedule

func DefaultRidgeSchedule() RidgeSchedule

DefaultRidgeSchedule returns a fresh copy of the default schedule: eight attempts, increments 10^(t-6) for t = 0..7 (1e-6 … 1e1), so a matrix that fails every attempt reports a total ridge of 11.111111.

type SPDFactor

type SPDFactor struct {
	// contains filtered or unexported fields
}

SPDFactor is the gonum-backed Cholesky factorisation of a symmetric positive-definite matrix, built by FactorSPD and reused across any number of Solve and Inverse calls. It owns its storage: mutating the Sym it was built from does not affect it.

func FactorSPD

func FactorSPD(s *Sym) (*SPDFactor, error)

FactorSPD returns the gonum-backed Cholesky factorisation of the symmetric positive-definite s. Only the upper triangle is consulted, which is all a Sym stores.

A matrix that is not positive definite — a pivot that is ≤ 0 or NaN, which covers a singular, indefinite or NaN-bearing matrix — is PULSE_MATRIX_SINGULAR (details "reason" = "not_positive_definite", "n", plus the singular diagnostics "rank", "condition_number" and, where identifiable, "dependent_indices" — see singularDiagnostics; none of them for a NaN-bearing matrix). There is no separate finiteness pre-check: an infinite element is handed to the factorisation as-is. A nil s is PULSE_MATRIX_SHAPE_MISMATCH; an empty s factors to an empty factor.

func (*SPDFactor) ConditionNumber

func (f *SPDFactor) ConditionNumber() float64

ConditionNumber returns gonum's 1-norm estimate of the factored matrix's condition number (1 for an empty factor). Solve and Inverse fail when it exceeds ConditionTolerance.

func (*SPDFactor) Inverse

func (f *SPDFactor) Inverse() (*Sym, error)

Inverse returns A⁻¹ for the factored A. A factor whose condition number exceeds ConditionTolerance, or whose triangular inverse breaks down, is PULSE_MATRIX_SINGULAR with the condition estimate in details "condition_number" (+Inf for a breakdown).

func (*SPDFactor) N

func (f *SPDFactor) N() int

N returns the order of the factored matrix.

func (*SPDFactor) Solve

func (f *SPDFactor) Solve(b *Vec) (*Vec, error)

Solve returns x with A·x = b for the factored A. A nil b or len(b) != N() is PULSE_MATRIX_SHAPE_MISMATCH. A factor whose condition number exceeds ConditionTolerance is PULSE_MATRIX_SINGULAR, with the estimate in details "condition_number".

type SVDResult

type SVDResult struct {
	U      *Matrix
	Values *Vec
	V      *Matrix
}

SVDResult is the THIN singular value decomposition m = U·diag(Values)·Vᵀ of a rows×cols matrix, with k = min(rows, cols): U is rows×k, V is cols×k, both with orthonormal columns.

Order: Values are descending; values within RankTolerance(rows, cols, σ_max) of their neighbour form one tie cluster ordered by the dominant component index of the V column (original variable order).

Sign: each pair (U column k, V column k) is flipped TOGETHER so the dominant component of the V column is positive (DominanceTolerance; a tie goes to the lowest index). Flipping both keeps U·Σ·Vᵀ unchanged.

func SVD

func SVD(m *Matrix) (*SVDResult, error)

SVD returns the thin singular value decomposition of m under the order and sign policy documented on SVDResult. A nil m is PULSE_MATRIX_SHAPE_MISMATCH; a non-finite element or a decomposition that does not converge is PULSE_MATRIX_SINGULAR (details "reason"). A matrix with a zero dimension decomposes to empty factors.

type Sym

type Sym struct {
	// contains filtered or unexported fields
}

Sym is an n×n symmetric matrix stored as its packed upper triangle in row order: (0,0) (0,1) … (0,n-1) (1,1) … (n-1,n-1). At(i, j) and At(j, i) read the same element. The zero value is an empty 0×0 matrix.

func InverseSPD

func InverseSPD(s *Sym) (*Sym, error)

InverseSPD returns s⁻¹ for the symmetric positive-definite s, computed as L⁻ᵀ·L⁻¹ from the reference Cholesky factor so the result is symmetric by construction. L⁻¹ is built column by column with a running subtraction over k ascending and a true division by the diagonal; each (i, j ≥ i) element of the inverse is then a running sum of L⁻¹[k][i]·L⁻¹[k][j] over k ascending from j. FMA-free.

A nil s is PULSE_MATRIX_SHAPE_MISMATCH; a factorisation failure is PULSE_MATRIX_SINGULAR.

func NewSym

func NewSym(n int, packed []float64) (*Sym, error)

NewSym returns an n×n symmetric matrix over a COPY of packed, the n(n+1)/2 upper-triangle elements in row order. A nil packed allocates a zero matrix. A negative n or a length mismatch is PULSE_MATRIX_SHAPE_MISMATCH.

func NewSymFromRows

func NewSymFromRows(rows [][]float64) (*Sym, error)

NewSymFromRows returns the symmetric matrix whose LOWER triangle is read from a square [][]float64: element (i, j) for j ≤ i is rows[i][j], and the strict upper triangle of the input is ignored. Reading the lower triangle is deliberate — it is exactly what a row-wise Cholesky consumes. A non-square or ragged input is PULSE_MATRIX_SHAPE_MISMATCH; an empty input is a 0×0 matrix.

func (*Sym) At

func (s *Sym) At(i, j int) float64

At returns element (i, j) == (j, i). It panics when out of range.

func (*Sym) N

func (s *Sym) N() int

N returns the order of the matrix.

func (*Sym) Set

func (s *Sym) Set(i, j int, v float64)

Set assigns element (i, j) and therefore (j, i). It panics when out of range.

func (*Sym) ToRows

func (s *Sym) ToRows() [][]float64

ToRows returns a copy of the full (both triangles) n×n matrix.

type SymEigenResult

type SymEigenResult struct {
	Values  *Vec
	Vectors *Matrix
}

SymEigenResult is the eigen-decomposition s = V·diag(Values)·Vᵀ of a symmetric matrix. Column k of Vectors is the unit eigenvector of Values[k].

Order: Values are descending. Eigenvalues within the rank tolerance RankTolerance(n, n, max|λ|) of their neighbour form one tie cluster; inside a cluster the vectors are ordered by the index of their dominant component (original variable order), lowest first.

Sign: every eigenvector is flipped so its dominant component (see DominanceTolerance; a tie goes to the lowest index) is positive.

func SymEigen

func SymEigen(s *Sym) (*SymEigenResult, error)

SymEigen returns the eigen-decomposition of the symmetric s under the order and sign policy documented on SymEigenResult. A nil s is PULSE_MATRIX_SHAPE_MISMATCH; a non-finite element, or a decomposition that does not converge, is PULSE_MATRIX_SINGULAR (details "reason"). An empty s decomposes to empty Values and Vectors.

type Vec

type Vec struct {
	// contains filtered or unexported fields
}

Vec is a dense float64 vector. The zero value is empty.

func NewVec

func NewVec(data []float64) *Vec

NewVec returns a vector over a COPY of data.

func Solve

func Solve(a *Matrix, b *Vec) (*Vec, error)

Solve solves the general square system a·x = b through the reference LU factorisation (FactorLU, then LUFactor.Solve): pure Go, FMA-free, bit-identical on every architecture. Use it for a non-symmetric or indefinite a; for a symmetric positive-definite a, SolveSPD is the cheaper reference route. To reuse one factorisation across several right-hand sides, call FactorLU once and Solve / SolveMatrix on the factor.

A nil operand, a non-square a or len(b) != a.Rows() is PULSE_MATRIX_SHAPE_MISMATCH; a singular a (see FactorLU's policy) or a NaN / infinite element is PULSE_MATRIX_SINGULAR — never a NaN result.

func SolveSPD

func SolveSPD(s *Sym, b *Vec) (*Vec, error)

SolveSPD solves s·x = b for the symmetric positive-definite s through its reference Cholesky factor: forward substitution L·y = b, then back substitution Lᵀ·x = y, each a running subtraction over k ascending followed by a true division by the diagonal. FMA-free.

A nil operand or len(b) != order(s) is PULSE_MATRIX_SHAPE_MISMATCH; a factorisation failure is PULSE_MATRIX_SINGULAR.

func (*Vec) At

func (v *Vec) At(i int) float64

At returns element i. It panics when out of range.

func (*Vec) Len

func (v *Vec) Len() int

Len returns the vector length.

func (*Vec) Set

func (v *Vec) Set(i int, x float64)

Set assigns element i. It panics when out of range.

func (*Vec) Slice

func (v *Vec) Slice() []float64

Slice returns a copy of the elements.

Directories

Path Synopsis
Package metric is Pulse's one registry of distance and similarity metrics: every metric is defined here once — its kernel, the vector kinds it suits and fits poorly, the direction(s) it can answer in, whether it needs fitted statistics, and its domain rules — and every consumer (the expression metric functions, MAT_DISTANCE, GROUP_KMEANS, the reference-comparison attributes and OVERLAY_CORRESPONDENCE's chi_square profile distances) resolves a metric name through it.
Package metric is Pulse's one registry of distance and similarity metrics: every metric is defined here once — its kernel, the vector kinds it suits and fits poorly, the direction(s) it can answer in, whether it needs fitted statistics, and its domain rules — and every consumer (the expression metric functions, MAT_DISTANCE, GROUP_KMEANS, the reference-comparison attributes and OVERLAY_CORRESPONDENCE's chi_square profile distances) resolves a metric name through it.

Jump to

Keyboard shortcuts

? : This menu
/ : Search site
f or F : Jump to
y or Y : Canonical URL