A Lightweight ALS Solver for Iterative GLS

When a GLS problem involves hundreds of equations, the $K × K$ covariance matrix becomes the computational bottleneck. A simple statistical remedy is to assume that most of the cross‑equation dependence can be captured by a handful of latent factors plus equation‑specific noise. This "low‑rank + diagonal" assumption slashes the number of unknowns from roughly $K^²$ to about $K×k$ parameters, where k (the latent factor rank) is much smaller than $K$. This is essentially poor man's multi-task learning: the shared latent factors let information flow across tasks (equations), improving estimates when tasks are related—without the complexity of neural networks. The model alone, however, does not guarantee speed: we still have to fit the parameters.

Installation

Install the library from PyPI:

pip install alsgls

For local development, clone the repo and use an editable install:

pip install -e .

Usage

from alsgls import ALSGLS, ALSGLSSystem, simulate_sur

Xs_tr, Y_tr, Xs_te, Y_te = simulate_sur(N_tr=240, N_te=120, K=60, p=3, k=4)

# Scikit-learn style estimator
est = ALSGLS(rank="auto", max_sweeps=12)
est.fit(Xs_tr, Y_tr)
test_score = est.score(Xs_te, Y_te)  # negative test NLL per observation

# Statsmodels-style system interface
system = {f"eq{j}": (Y_tr[:, j], Xs_tr[j]) for j in range(Y_tr.shape[1])}
sys_model = ALSGLSSystem(system, rank="auto")
sys_results = sys_model.fit()
params = sys_results.params_as_series()  # pandas optional

Rank Selection

The package supports automatic rank selection via BIC or cross-validation:

from alsgls import ALSGLS

# BIC-based rank selection
est = ALSGLS(rank="bic", max_sweeps=15)
est.fit(Xs, Y)
print(f"Selected rank: {est.rank_}")

# Cross-validation rank selection
est = ALSGLS(rank="cv", cv_folds=5, cv_random_state=42)
est.fit(Xs, Y)
print(f"Selected rank: {est.rank_}")

Real Data Example

See examples/real_data_fama_french.py for a demonstration using Fama-French 49 industry portfolios.

The benchmarks/compare_sur.py script contrasts ALS-GLS with statsmodels and linearmodels SUR implementations on matched simulation grids while recording peak memory (via Memray, Fil, or the POSIX RSS high-water mark).

Documentation and notebooks

Background material and reproducible experiments are available in the notebooks under als_sim/, such as als_sim/als_comparison.ipynb and als_sim/als_sur.ipynb.

Type-Safe ALS Solver

This package provides a modern, type-safe implementation of Alternating-Least-Squares (ALS) for low-rank GLS problems. The Woodbury identity reduces the expensive inverse to a tiny k × k system, and the β-update can be written without explicitly forming dense matrices.

Inference:

• Standard errors (bse), t-statistics (tvalues), p-values (pvalues), confidence intervals (conf_int()) and summary tables (summary()), statsmodels-style. By default these carry the Kackar–Harville correction for Σ being estimated — (X'Σ̂⁻¹X)⁻¹ + Λ after the degrees-of-freedom rescale — which no other SUR implementation applies. results.covariance("plugin") gives the uncorrected plug-in that linearmodels, systemfit and Stata report. Measured on a 4-equation system, reported SE over actual spread: plug-in 0.85, corrected 0.89 at n = 20; 0.94 / 0.96 at n = 40; both ~1.0 by n = 200.
• Calibrated small-sample inference via results.bootstrap(B=999), which refits the whole model on each replicate and returns percentile-t intervals and bootstrap-t p-values. This is the object to report when n is small relative to the number of equations; see docs/formal_methods.md §7 for the measurements.

New in v1.1.0:

• Rank selection: BIC and cross-validation for automatic rank selection
• Gradient-based factor update: Cleaner theory, same convergence guarantees
• Real-world example: Fama-French 49 industry portfolios demonstration
• Formal methods documentation: Rigorous mathematical foundations

Core features:

• Full type safety with mypy compliance and comprehensive type hints
• Numerically stable implementation using Cholesky factorization throughout
• Clean API with single computational path and enhanced error messages
• Memory efficient with O(K k) complexity, converging in 5–6 sweeps

Rule of thumb: if your GLS routine keeps looping between $\beta$ and a fresh $\hat{\Sigma}$, the ALS approach yields the same statistical fit with an order‑of‑magnitude smaller memory footprint and better numerical stability.

Beyond SUR: where the idea travels

Random‑effects models, feasible GLS with estimated heteroskedastic weights, optimal‑weight GMM, and spatial autoregressive GLS all iterate β ↔ Σ̂. Each can adopt the same ALS trick: treat the weight matrix as low‑rank + diagonal, invert only the k × k core, and avoid the dense K × K algebra. Memory savings in published examples range from 5× to 20×, depending on k.

A concrete case‑study: Seemingly‑Unrelated Regressions

To demonstrate performance, we benchmark ALS against traditional methods with N = 300 observations, three regressors, rank‑3 factors, and K ranging from 50 to 120 equations. The largest array that traditional methods need is the dense Σ⁻¹ (K×K), whereas ALS's largest is the skinny factor matrix F (K×k).

K	β‑RMSE EM	β‑RMSE ALS	Peak MB EM	Peak MB ALS	Memory ratio
50	0.021	0.021	0.020	0.002	10×
80	0.020	0.020	0.051	0.003	17×
120	0.020	0.020	0.115	0.004	29×

The ALS implementation achieves the same statistical performance while using only a few megabytes of memory, providing substantial computational advantages for large systems.

Defaults, tuning knobs, and failure modes

• Rank (k) – By default the high-level APIs pick min(8, ceil(K / 10)), a conservative fraction of the number of equations. Increase rank if the cross-equation correlation matrix is slow to decay; decrease it when the diagonal dominates.
• Ridge term (lam_B) – Defaults to 1e-3 on the regression update, and is applied relative to the residual variance scale so the fit does not depend on the units of Y. Raise it (e.g. 1e-2) if CG struggles to converge. There is no penalty on the factor loadings: the Σ-step is the exact conditional solution, so there is nothing for one to regularise.
• Noise floor (d_floor) – Keeps the diagonal component positive; the default 1e-8 is a fraction of the mean residual variance, not an absolute variance, so it transforms correctly under a change of units. Increase it in highly ill-conditioned settings.
• Stopping criteria – ALS stops when the relative drop in NLL per sweep is below 1e-6 (configurable via rel_tol) or after max_sweeps. Inspect info["nll_trace"] to diagnose stagnation.
• Possible failures – Large condition numbers or nearly-collinear regressors can make the β-step CG solve slow; adjust cg_tol/cg_maxit, add stronger ridge, or re-scale predictors. info["sigma_iters"] reports how many inner iterations each Σ-step needed; counts that sit at the cap mean the alternation is crawling, which happens when some diagonal variances are near zero (a Heywood case) and usually indicates the factor rank is too large relative to the sample size.