alsgls package

Lightweight low-rank-plus-diagonal GLS/SUR estimation via ALS.

class alsgls.ALSGLS(*, rank='auto', rank_candidates=None, cv_folds=5, cv_random_state=None, lam_B=0.001, max_sweeps=12, rel_tol=1e-06, d_floor=1e-08, cg_maxit=800, cg_tol=3e-07)[source]

Bases: object

Scikit-learn style estimator for low-rank GLS via ALS.

Parameters:
get_params(deep=True)[source]

Return the estimator hyperparameters (scikit-learn protocol).

Parameters:

deep (bool)

Return type:

dict[str, Any]

set_params(**params)[source]

Set estimator hyperparameters in place (scikit-learn protocol).

Parameters:

params (Any)

Return type:

ALSGLS

als_kwargs()[source]

The keyword arguments this estimator passes to als_gls().

Return type:

dict[str, Any]

fit(Xs, Y)[source]

Fit the GLS system, selecting the rank first when requested.

Parameters:
Return type:

ALSGLS

predict(Xs)[source]

Predict responses for new design matrices, one per equation.

Parameters:

Xs (Sequence[Any])

Return type:

ndarray

score(Xs, Y)[source]

Return the negative Gaussian NLL per row of Y under the fitted Sigma.

Higher is better, as scikit-learn’s score protocol requires. This scores the whole fitted covariance, not just the conditional mean, so it is not the negative MSE; use alsgls.mse for that.

Parameters:
Return type:

float

covariance(method='kh')[source]

Covariance of the coefficients by the named method.

"kh" (default) is the Kackar-Harville corrected covariance, which accounts for Sigma being estimated; "plugin" is the df-rescaled plug-in that linearmodels, systemfit and Stata report. Measured on the test suite’s 4-equation, rank-1 fixture, the reported standard error is 0.85 (plug-in) and 0.89 (kh) of the actual spread at n = 20, and 0.94 / 0.96 at n = 40. Neither is calibrated at n = 20; for that use bootstrap().

Parameters:

method (str) – "kh" or "plugin".

Returns:

The (p_total, p_total) covariance.

Raises:

RuntimeError – If the training designs were not stored.

Return type:

ndarray

property cov_params_: ndarray

The default coefficient covariance, covariance("kh").

bootstrap(B=999, method='parametric', seed=None)[source]

Bootstrap the whole fit and return studentised (bootstrap-t) inference.

Each replicate refits F, D and beta on resampled data, which is what captures the part of the sampling variance the plug-in cov_params_ cannot: that Sigma is estimated. The calibrated object is the percentile-t interval, not the bootstrap standard error; see BootstrapResults.

Parameters:

  • B (int) – Number of replicates. 999 makes 0.05 * (B + 1) an integer, which is what a 95% percentile interval wants; 199 is enough if only the standard errors are of interest.

  • method (str) – "parametric" draws errors from the fitted F F' + diag(D); "wild" multiplies each residual row by a Rademacher sign; "residual" resamples whole residual rows. The last two make no distributional assumption.

  • seed (int | None) – Seed for the replicate stream; the same seed gives the same replicates.

Returns:

The bootstrap distribution and the inference built on it.

Return type:

BootstrapResults

predict_interval(Xs, alpha=0.05, return_type='prediction')[source]

Compute prediction or confidence intervals for new observations.

Parameters:

  • Xs (Sequence[Any]) – List of design matrices [X_0, …, X_{K-1}] for new observations.

  • alpha (float) – Significance level (default 0.05 for 95% intervals).

  • return_type (str) – “prediction” for prediction intervals (includes residual variance), “confidence” for confidence intervals (variance of E[y|X] only).

Returns:

{“mean”: (N, K), “lower”: (N, K), “upper”: (N, K)}

Return type:

dict

Raises:

ValueError – If return_type is neither "prediction" nor "confidence", or the design matrices do not match the fitted model in count or in number of columns.

class alsgls.ALSGLSSystem(system, *, rank='auto', lam_B=0.001, max_sweeps=12, rel_tol=1e-06, d_floor=1e-08, cg_maxit=800, cg_tol=3e-07)[source]

Bases: object

Statsmodels-style system container for ALS GLS fitting.

Parameters:
property equations: list[_SystemEquation]

The parsed equations, in system order.

property nobs: int

Number of observations shared by every equation.

property keqs: int

Number of equations in the system.

as_arrays()[source]

Return the system as (list of X_j, stacked Y) arrays.

Return type:

tuple[list[ndarray], ndarray]

fit()[source]

Fit the system with ALS-GLS and return a results container.

Return type:

ALSGLSSystemResults

class alsgls.ALSGLSSystemResults(model, estimator)[source]

Bases: object

Lightweight results container mimicking statsmodels outputs.

Parameters:
params_as_series()[source]

Return coefficients as a pandas Series with (equation, variable) index.

predict(exog=None)[source]

Predict fitted values, for the training design or new exog.

Parameters:

exog (Mapping[Any, Any] | Sequence[Any] | None)

Return type:

ndarray

summary_dict()[source]

Return the headline fit statistics as a plain dict.

Return type:

dict[str, Any]

covariance(method='kh')[source]

Covariance of the coefficients by the named method.

"kh" (default) is the Kackar-Harville corrected covariance, which accounts for Sigma being estimated; "plugin" is the df-rescaled plug-in that linearmodels, systemfit and Stata report. Measured on the test suite’s 4-equation, rank-1 fixture, the reported standard error is 0.85 (plug-in) and 0.89 (kh) of the actual spread at n = 20, and 0.94 / 0.96 at n = 40. Neither is calibrated at n = 20; for that use bootstrap().

Parameters:

method (str) – "kh" or "plugin".

Returns:

The (p_total, p_total) covariance.

Return type:

ndarray

property cov_params: ndarray

The default coefficient covariance, covariance("kh").

property bse: ndarray

Standard errors of parameter estimates (sqrt of diagonal of cov_params).

property tvalues: ndarray

β = 0.

Type:

t-statistics for H₀

property df_resid: int

nobs * keqs - n_params.

Type:

Residual degrees of freedom

property pvalues: ndarray

β = 0.

Type:

Two-sided p-values for H₀

conf_int(alpha=0.05)[source]

Confidence intervals for parameters.

Parameters:

alpha (float) – Significance level (default 0.05 for 95% CI)

Returns:

(n_params, 2) array with lower and upper bounds

Return type:

ci

get_prediction(exog=None, alpha=0.05)[source]

Get prediction results with standard errors and intervals.

Parameters:

  • exog (Mapping[Any, Any] | Sequence[Any] | None) – Design matrices for new observations. If None, uses the training data. Can be a dict {eq_name: X_j} or a list [X_0, …, X_{K-1}].

  • alpha (float) – Default significance level for intervals (default 0.05).

Returns:

Object with predicted_mean, se_mean, se_obs, and interval methods.

Return type:

PredictionResults

Raises:

ValueError – If design matrices are not supplied for every equation, or one has the wrong number of columns.

bootstrap(B=999, method='parametric', seed=None)[source]

Bootstrap the whole fit and return studentised (bootstrap-t) inference.

Each replicate refits F, D and beta on resampled data, which is what captures the part of the sampling variance the plug-in cov_params cannot: that Sigma is estimated. The calibrated object is the percentile-t interval, not the bootstrap standard error; see BootstrapResults.

Parameters:

  • B (int) – Number of replicates. 999 makes 0.05 * (B + 1) an integer, which is what a 95% percentile interval wants; 199 is enough if only the standard errors are of interest.

  • method (str) – "parametric" draws errors from the fitted F F' + diag(D); "wild" multiplies each residual row by a Rademacher sign; "residual" resamples whole residual rows. The last two make no distributional assumption.

  • seed (int | None) – Seed for the replicate stream; the same seed gives the same replicates.

Returns:

The bootstrap distribution and the inference built on it.

Return type:

BootstrapResults

summary(alpha=0.05)[source]

Text summary of estimation results (statsmodels-style).

Parameters:

alpha (float) – Significance level for confidence intervals

Returns:

Formatted summary table

Return type:

str

class alsgls.BootstrapResults(params, se_plugin, estimates, tstats, method, seed, df)[source]

Bases: object

Bootstrap distribution of the coefficients and the inference built on it.

Parameters:

  • params (ndarray)

  • se_plugin (ndarray)

  • estimates (ndarray)

  • tstats (ndarray)

  • method (str)

  • seed (int | None)

  • df (int)

params: ndarray

The parent fit’s coefficients, length p_total.

se_plugin: ndarray

The parent fit’s plug-in standard errors.

estimates: ndarray

Coefficients from each replicate, (B, p_total).

tstats: ndarray

Studentised deviations (beta(b) - beta_hat) / se(b), (B, p_total).

method: str

Which resampling scheme produced them.

seed: int | None

The seed the replicate stream was spawned from.

df: int

Residual degrees of freedom of the parent fit.

property B: int

Number of replicates.

property bse: ndarray

the spread of beta across replicates.

Better than the plug-in, but still biased down (Freedman and Peters measured 20-30% in their design), because each replicate is generated from a Sigma_hat that is itself too small. The interval, not this number, is the calibrated object.

Type:

Bootstrap standard errors

property tvalues: ndarray

params / se_plugin, the statistic the bootstrap-t p-value refers to.

conf_int(alpha=0.05)[source]

Percentile-t (bootstrap-t) confidence intervals.

[beta_hat - t*_{1-alpha/2} se, beta_hat - t*_{alpha/2} se] with the quantiles taken from the studentised replicates and se the parent plug-in. Asymmetric when the studentised distribution is.

Parameters:

alpha (float) – Significance level.

Returns:

(p_total, 2) array of lower and upper bounds.

Return type:

ndarray

property pvalues: ndarray

Two-sided bootstrap-t p-values for H0: beta = 0.

The share of replicates whose studentised deviation is at least as large in absolute value as the observed t, with the +1 continuity adjustment so B replicates cannot report zero.

summary(alpha=0.05)[source]

Text summary in the same layout as ALSGLSSystemResults.summary.

Parameters:

alpha (float) – Significance level for the intervals.

Returns:

The formatted table.

Return type:

str

class alsgls.PredictionResults(predicted_mean, se_mean, se_obs, _df, _alpha_default=0.05)[source]

Bases: object

Container for prediction results with intervals.

Parameters:

  • predicted_mean (ndarray)

  • se_mean (ndarray)

  • se_obs (ndarray)

  • _df (int)

  • _alpha_default (float)

predicted_mean: ndarray
se_mean: ndarray
se_obs: ndarray
conf_int_mean(alpha=None)[source]

Confidence intervals for E[y|X].

Parameters:

alpha (float | None) – Significance level (default 0.05 for 95% CI)

Returns:

(N, K, 2) array with lower and upper bounds

Return type:

ci

conf_int_obs(alpha=None)[source]

Prediction intervals for y|X (includes residual variance).

Parameters:

alpha (float | None) – Significance level (default 0.05 for 95% CI)

Returns:

(N, K, 2) array with lower and upper bounds

Return type:

ci

alsgls.XB_from_Blist(Xs, B_list)[source]

Return N x K matrix of predictions.

Parameters:

  • Xs (list[ndarray])

  • B_list (list[ndarray])

Return type:

ndarray

alsgls.als_gls(Xs, Y, k, lam_B=0.001, sweeps=8, d_floor=1e-08, cg_maxit=800, cg_tol=3e-07, *, rel_tol=1e-06, init_D=None)[source]

Alternating GLS with a low-rank-plus-diagonal covariance.

Alternates two exact steps until the likelihood stops falling: beta by matrix-free conjugate gradients at the current Sigma, then Sigma by the closed-form factor-analysis update in _sigma_step(). Woodbury is used throughout and no K x K matrix is formed.

Parameters:

  • Xs (list[ndarray]) – One design matrix per equation, each (N, p_j). The equations may have different numbers of regressors.

  • Y (ndarray) – Responses, (N, K), one column per equation.

  • k (int) – Rank of the latent factor block. The covariance is modelled as F F^T + D with F of shape (K, k), so k controls how much cross-equation dependence is shared; k << K is the point.

  • lam_B (float) – Ridge penalty on the coefficients beta, expressed relative to the mean initial residual variance rather than in absolute units. An absolute penalty would make the fit depend on the units of Y: X' Sigma^-1 X scales as 1/s^2 under Y -> sY while a fixed lam_B does not, so at lam_B = 1e-3 scaling Y by 1e4 moved every coefficient by 100%. Dividing by the variance scale makes the penalty transform correctly and the whole fit equivariant.

  • sweeps (int) – Maximum number of alternating passes over beta and Sigma.

  • d_floor (float) – Floor on each diagonal variance, as a fraction of the mean initial residual variance rather than an absolute variance, for the same reason lam_B is relative: the true D scales as s^2 under Y -> sY, so an absolute floor binds on every entry once Y is small enough.

  • cg_maxit (int) – Iteration cap for the conjugate-gradient solve in the beta step.

  • cg_tol (float) – Relative residual tolerance for that solve.

  • rel_tol (float) – Relative decrease in the negative log-likelihood below which the sweeps stop early.

  • init_D (ndarray | None) – Starting diagonal variances for the first Sigma step, length K. Defaults to the residual variances of the initial fit. A bootstrap replicate passes the parent fit’s D here: the replicate’s answer is close to the parent’s, so starting from it skips most of the inner alternation. Only D is needed, because the Sigma step recomputes F from D in closed form.

Returns:

info includes p_list, cg (the final beta solve), nll_trace (per sweep, non-increasing, and equal to nll_per_row at the returned parameters), nll_beta_trace (post-beta, per sweep), sigma_iters (inner alternation counts), and var_ref/lam_B_eff (the variance scale the relative penalties were resolved against, and the resulting absolute ridge).

Return type:

B_list, F, D, mem_MB_est, info

Raises:

  • ValueError – If an argument is outside its domain, if Xs or Y holds a non-finite entry, or if their shapes disagree.

  • np.linalg.LinAlgError – If a Cholesky factorisation of the Woodbury core fails, which means the current Sigma is not positive definite.

alsgls.mse(Y, Yhat)[source]

Mean squared error between two matrices.

Parameters:

  • Y (ndarray)

  • Yhat (ndarray)

Return type:

float

alsgls.nll_per_row(R, F, D)[source]

Negative log-likelihood per row for a residual matrix.

Computed for R under Σ = F F^T + diag(D) with Gaussian errors.

Parameters:

  • R (ndarray) – Residual matrix, (N, K).

  • F (ndarray) – Factor loadings, (K, k).

  • D (ndarray) – Diagonal variances, length K.

Returns:

0.5 * [ tr(R Σ^{-1} R^T)/N + log det(Σ) + K log(2π) ] where N is the number of rows in R.

Return type:

float

alsgls.select_rank_bic(Xs, Y, k_candidates=None, **als_kwargs)[source]

Select rank by minimizing BIC (Bayesian Information Criterion).

BIC(k) = 2 * N * nll_per_row + n_params * log(N)

which is the usual -2 * loglik + n_params * log(N), since N * nll_per_row is the negative log-likelihood of the sample. The value used to be reported at half this, so it was not comparable with the bic attribute of statsmodels or any other package; the rank chosen is unchanged, because halving is monotone.

n_params counts the free parameters of the whole fitted model: the K*k factor loadings in F less the k*(k-1)/2 orthogonal rotations F -> F Q that leave F F^T unchanged and so are not identified, the K diagonal variances in D, and the sum(p_j) regression coefficients, since nll_per_row is evaluated at the fitted beta and a BIC has to charge for it. This is the standard factor-analysis count; R’s factanal reports the complementary df = ((K-k)^2 - K - k) / 2.

A rank whose fit raised carries the message in its error key.

The count used to be K*(k+1) + k, which neither subtracted the rotational redundancy nor charged for beta. Only the first term varies with k, so on the fixtures tested the selected rank is unchanged; the reported value was wrong either way.

Parameters:

  • Xs (list[ndarray]) – Design matrices for each equation.

  • Y (ndarray) – Response matrix.

  • k_candidates (list[int] | None) – Candidate ranks to evaluate. Defaults to range(1, min(K//2, 12)+1).

  • **als_kwargs (Any) – Additional arguments passed to als_gls().

Returns:

The rank with minimum BIC. results: Per-rank dicts of ‘k’, ‘nll’, ‘bic’, ‘n_params’, ‘converged’.

Return type:

best_k

Raises:

RuntimeError – If no candidate rank produced a usable fit.

alsgls.select_rank_cv(Xs, Y, k_candidates=None, n_folds=5, random_state=None, **als_kwargs)[source]

Select rank by k-fold cross-validation on validation NLL.

Parameters:

  • Xs (list[ndarray]) – Design matrices for each equation.

  • Y (ndarray) – Response matrix.

  • k_candidates (list[int] | None) – Candidate ranks to evaluate. Defaults to range(1, min(K//2, 12)+1).

  • n_folds (int) – Number of cross-validation folds.

  • random_state (int | Generator | None) – Random state for reproducible fold splits.

  • **als_kwargs (Any) – Additional arguments passed to als_gls().

Returns:

The rank with minimum mean CV NLL. results: Per-rank results containing ‘k’, ‘cv_nll’, ‘cv_std’, ‘fold_nlls’.

Return type:

best_k

Raises:

  • ValueError – If n_folds is below 2 or exceeds the number of rows.

  • RuntimeError – If no candidate rank produced a usable fit.

alsgls.simulate_gls(N_tr, N_te, p_list, k, seed=0)[source]

Simulate a generalized least squares (GLS) dataset.

This variant allows each response equation to have its own number of features as specified by p_list.

Parameters:

  • N_tr (int) – Number of training samples.

  • N_te (int) – Number of test samples.

  • p_list (list[int]) – Number of features for each equation.

  • k (int) – Latent factor dimension controlling correlated noise.

  • seed (int) – Seed for the NumPy random number generator. Defaults to 0.

Returns:

per-equation feature matrices with X_tr[j] of shape (N_tr, p_list[j]) and X_te[j] of shape (N_te, p_list[j]), and responses of shape (N_tr, K) and (N_te, K) where K = len(p_list).

Return type:

A tuple (X_tr, Y_tr, X_te, Y_te)

Notes

Randomness is controlled via numpy.random.default_rng(seed); pass a different seed for different simulations.

Examples

>>> from alsgls import simulate_gls
>>> p_list = [3, 5, 2]
>>> Xtr, Ytr, Xte, Yte = simulate_gls(100, 20, p_list, k=2, seed=0)
alsgls.simulate_sur(N_tr, N_te, K, p, k, seed=0)[source]

Simulate a Seemingly Unrelated Regression (SUR) dataset.

Parameters:

  • N_tr (int) – Number of training samples.

  • N_te (int) – Number of test samples.

  • K (int) – Number of response equations.

  • p (int) – Number of features per equation.

  • k (int) – Latent factor dimension controlling correlated noise.

  • seed (int) – Seed for the NumPy random number generator. Defaults to 0.

Returns:

lists of per-equation feature matrices of shape (N_tr, p) and (N_te, p), and response matrices of shape (N_tr, K) and (N_te, K).

Return type:

A tuple (X_tr, Y_tr, X_te, Y_te)

Notes

Randomness is controlled via numpy.random.default_rng(seed); pass a different seed for different simulations.

Examples

>>> from alsgls import simulate_sur
>>> Xtr, Ytr, Xte, Yte = simulate_sur(100, 20, K=3, p=5, k=2, seed=42)

Submodules

alsgls.als module

Alternating-least-squares solver for low-rank-plus-diagonal GLS.

alsgls.als.als_gls(Xs, Y, k, lam_B=0.001, sweeps=8, d_floor=1e-08, cg_maxit=800, cg_tol=3e-07, *, rel_tol=1e-06, init_D=None)[source]

Alternating GLS with a low-rank-plus-diagonal covariance.

Alternates two exact steps until the likelihood stops falling: beta by matrix-free conjugate gradients at the current Sigma, then Sigma by the closed-form factor-analysis update in _sigma_step(). Woodbury is used throughout and no K x K matrix is formed.

Parameters:

  • Xs (list[ndarray]) – One design matrix per equation, each (N, p_j). The equations may have different numbers of regressors.

  • Y (ndarray) – Responses, (N, K), one column per equation.

  • k (int) – Rank of the latent factor block. The covariance is modelled as F F^T + D with F of shape (K, k), so k controls how much cross-equation dependence is shared; k << K is the point.

  • lam_B (float) – Ridge penalty on the coefficients beta, expressed relative to the mean initial residual variance rather than in absolute units. An absolute penalty would make the fit depend on the units of Y: X' Sigma^-1 X scales as 1/s^2 under Y -> sY while a fixed lam_B does not, so at lam_B = 1e-3 scaling Y by 1e4 moved every coefficient by 100%. Dividing by the variance scale makes the penalty transform correctly and the whole fit equivariant.

  • sweeps (int) – Maximum number of alternating passes over beta and Sigma.

  • d_floor (float) – Floor on each diagonal variance, as a fraction of the mean initial residual variance rather than an absolute variance, for the same reason lam_B is relative: the true D scales as s^2 under Y -> sY, so an absolute floor binds on every entry once Y is small enough.

  • cg_maxit (int) – Iteration cap for the conjugate-gradient solve in the beta step.

  • cg_tol (float) – Relative residual tolerance for that solve.

  • rel_tol (float) – Relative decrease in the negative log-likelihood below which the sweeps stop early.

  • init_D (ndarray | None) – Starting diagonal variances for the first Sigma step, length K. Defaults to the residual variances of the initial fit. A bootstrap replicate passes the parent fit’s D here: the replicate’s answer is close to the parent’s, so starting from it skips most of the inner alternation. Only D is needed, because the Sigma step recomputes F from D in closed form.

Returns:

info includes p_list, cg (the final beta solve), nll_trace (per sweep, non-increasing, and equal to nll_per_row at the returned parameters), nll_beta_trace (post-beta, per sweep), sigma_iters (inner alternation counts), and var_ref/lam_B_eff (the variance scale the relative penalties were resolved against, and the resulting absolute ridge).

Return type:

B_list, F, D, mem_MB_est, info

Raises:

  • ValueError – If an argument is outside its domain, if Xs or Y holds a non-finite entry, or if their shapes disagree.

  • np.linalg.LinAlgError – If a Cholesky factorisation of the Woodbury core fails, which means the current Sigma is not positive definite.

alsgls.ops module

Woodbury-based linear algebra kernels for the ALS-GLS solver.

alsgls.ops.woodbury_chol(F, D)[source]

Return (Dinv, C_chol) with the Cholesky factor of C = I + F^T D^{-1} F.

Intended for numerically stable downstream solves that avoid forming C^{-1} explicitly.

Parameters:

  • F (ndarray)

  • D (ndarray)

Return type:

tuple[ndarray, ndarray]

alsgls.ops.apply_siginv_to_matrix(M, F, D, *, Dinv=None, C_chol)[source]

Right-multiply an (NxK) matrix M by Σ^{-1} using Woodbury.

Σ = F F^T + diag(D), which is never formed densely.

Uses numerically stable Cholesky factorization approach.

Parameters:

  • M (ndarray) – (NxK) matrix to right-multiply by Σ^{-1}

  • F (ndarray) – (Kxk) factor loadings matrix

  • D (ndarray) – (K,) diagonal noise variances

  • Dinv (ndarray | None) – Pre-computed 1/D. If None, computed from D.

  • C_chol (ndarray) – Cholesky factor of C = I + F^T D^{-1} F

Returns:

M @ Σ^{-1}

Return type:

np.ndarray

alsgls.ops.stack_B_list(B_list)[source]

Stack list of (p_jx1) blocks into a flat vector.

Parameters:

B_list (list[ndarray])

Return type:

ndarray

alsgls.ops.unstack_B_vec(bvec, p_list)[source]

Inverse of stack: vector -> list of (p_jx1).

Parameters:

  • bvec (ndarray)

  • p_list (list[int])

Return type:

list[ndarray]

alsgls.ops.XB_from_Blist(Xs, B_list)[source]

Return N x K matrix of predictions.

Parameters:

  • Xs (list[ndarray])

  • B_list (list[ndarray])

Return type:

ndarray

alsgls.ops.cg_solve(operator_mv, b, x0=None, maxit=500, tol=1e-07, M_pre=None)[source]

Conjugate gradient for SPD operator A (matrix-free).

Parameters:

  • operator_mv (Callable[[np.ndarray], np.ndarray]) – Function that returns A @ x for a given x.

  • b (np.ndarray) – Right-hand side.

  • x0 (np.ndarray | None) – Initial guess.

  • maxit (int) – Maximum CG iterations.

  • tol (float) – Relative residual tolerance.

  • M_pre (Callable[[np.ndarray], np.ndarray] | None) – Preconditioner application: returns M^{-1} @ r.

Returns:

Approximate solution. info: Iterations and final residual norm.

Return type:

x

Raises:

ValueError – If the operator or the preconditioner turns out not to be positive definite, which conjugate gradients requires.

alsgls.ops.siginv_diag(F, Dinv, C_chol)[source]

Compute the diagonal of Σ^{-1} without forming the inverse.

Σ^{-1} = D^{-1} - D^{-1} F C^{-1} F^T D^{-1}, evaluated from Dinv and the Cholesky factor of C.

Parameters:

  • F (ndarray) – Factor loadings, (K, k).

  • Dinv (ndarray) – Reciprocal of the diagonal variances, length K.

  • C_chol (ndarray) – Cholesky factor of the k x k Woodbury core.

Returns:

The diagonal entries of Σ^{-1}.

Return type:

diag_Sinv

alsgls.ops.apply_siginv_F(F, Dinv, C_chol)[source]

Compute Σ^{-1} @ F efficiently using Woodbury.

Σ^{-1} @ F = D^{-1} F - D^{-1} F C^{-1} F^T D^{-1} F

Parameters:

  • F (ndarray) – Factor loadings

  • Dinv (ndarray) – Inverse diagonal

  • C_chol (ndarray) – Cholesky factor of C = I + F^T D^{-1} F

Returns:

Σ^{-1} @ F

Return type:

SinvF

class alsgls.ops.GramBlocks(Xs)[source]

Bases: object

The cross-products X_j' X_l laid out for fast block assembly.

Every precision-weighted quantity in the SUR system is a block matrix whose (j, l) block is a scalar times X_j' X_l. Rather than keep the K^2 blocks as a list and loop over them per assembly – which at K = 100 and 500 assemblies is five million Python-level block writes, and was the whole cost of the Kackar-Harville correction – the blocks are tiled once into a (p_total, p_total) template, and an assembly is one gather and one elementwise product.

Parameters:

Xs (Sequence[np.ndarray])

template

(p_total, p_total) with block (j, l) equal to X_j' X_l.

eq_of_row

Equation index of each row of the stacked coefficient vector.

p_list

Number of regressors per equation.

alsgls.ops.gram_blocks(Xs)[source]

Prepare GramBlocks for the design.

Parameters:

Xs (Sequence[np.ndarray]) – One design matrix per equation.

Returns:

The tiled cross-products.

Return type:

GramBlocks

alsgls.ops.assemble_blocks(M, gram)[source]

X' (M (x) I_n) X for the block-diagonal SUR design and a K x K M.

Its (j, l) block is M[j, l] * X_j' X_l. With M = Sigma^-1 this is the GLS normal matrix; with M = Sigma^-1 (dSigma/dtheta) Sigma^-1 it is a Kackar-Harville P matrix, and so on.

Parameters:
Returns:

The (p_total, p_total) assembled matrix.

Return type:

ndarray

alsgls.ops.df_rescaled(F, D, n, p_list)[source]

Rescale a fitted low-rank covariance for the residual degrees of freedom.

Sigma_hat is estimated from residuals of a fitted model, so its entries are biased toward zero by the fitting; the correction every SUR package applies is Sigma_ij * n / sqrt((n - p_i)(n - p_j)) (linearmodels debiased=True, R systemfit methodResidCov="geomean", Stata sureg, dfk). That elementwise scaling is diag(sqrt(c)) Sigma diag(sqrt(c)) with c_i = n / (n - p_i), and since

diag(sqrt(c)) (F F’ + diag(D)) diag(sqrt(c))

= (diag(sqrt(c)) F)(diag(sqrt(c)) F)’ + diag(c * D),

it preserves the low-rank-plus-diagonal structure exactly. Scale row i of F by sqrt(c_i) and D_i by c_i.

This corrects the bias in Sigma_hat itself. It does nothing about the variance of Sigma_hat, which is the larger part of the finite-sample shortfall in the plug-in standard errors; see BootstrapResults.

Parameters:

  • F (np.ndarray) – Factor loadings, (K, k).

  • D (np.ndarray) – Diagonal variances, length K.

  • n (int) – Number of rows each equation was fitted on.

  • p_list (Sequence[int]) – Number of regressors in each equation, length K.

Returns:

The rescaled (F, D).

Raises:

ValueError – If any equation has no residual degrees of freedom.

Return type:

tuple[np.ndarray, np.ndarray]

alsgls.ops.compute_XtSigmaInvX(Xs, F, D, lam_B=0.0)[source]

Compute (X’Σ⁻¹X + λI) using the Woodbury identity.

For GLS with Σ = FF’ + diag(D), we need the precision-weighted design matrix cross-product for computing coefficient standard errors:

Var(β̂) = (X’Σ⁻¹X + λI)⁻¹

This is exact only when lam_B is 0. For a ridge estimator the variance is the sandwich A⁻¹ (X'Σ⁻¹X) A⁻¹ with A = X'Σ⁻¹X + λI, and the form above understates it; at the default lam_B = 1e-3 the difference is negligible, but it grows with the penalty. Σ is also treated as known rather than estimated, which is the usual feasible-GLS convention.

Using Woodbury: Σ⁻¹ = D⁻¹ - D⁻¹F C⁻¹ F’D⁻¹ where C = I + F’D⁻¹F

The block structure gives:

[X’Σ⁻¹X]_{jl} = X_j’ [Σ⁻¹]_{jl} X_l

Parameters:

  • Xs (list[ndarray]) – List of design matrices [X_0, …, X_{K-1}] where X_j is (N, p_j)

  • F (ndarray) – (K, k) factor loadings matrix

  • D (ndarray) – (K,) diagonal noise variances

  • lam_B (float) – Ridge penalty to add to diagonal (for regularization)

Returns:

(p_total, p_total) matrix where p_total = sum(p_j)

Return type:

XtSinvX

alsgls.ops.compute_prediction_variance(Xs, F, D, cov_params, include_residual=True)[source]

Compute prediction variances for new observations.

For each observation i and equation j: - Var(E[y_j|X]) = X_j[i,:] @ Cov(β̂_j) @ X_j[i,:] - Var(y_j|X) = Var(E[y_j|X]) + Σ_jj where Σ_jj = ||F[j,:]||² + D[j]

Parameters:

  • Xs (Sequence[np.ndarray]) – List of design matrices [X_0, …, X_{K-1}] where X_j is (N_new, p_j)

  • F (np.ndarray) – (K, k) factor loadings matrix

  • D (np.ndarray) – (K,) diagonal noise variances

  • cov_params (np.ndarray) – (p_total, p_total) covariance matrix of parameter estimates

  • include_residual (bool) – If True, add Σ_jj (residual variance) for prediction intervals. If False, return only the variance of the mean prediction (confidence intervals).

Returns:

(N_new, K) array of prediction variances

Return type:

var_pred

Raises:

ValueError – If Xs is empty.

alsgls.api module

High-level estimator APIs for ALS-based GLS fitting.

alsgls.api.COV_METHODS = ('kh', 'plugin')

Covariance estimators covariance() accepts.

class alsgls.api.ALSGLS(*, rank='auto', rank_candidates=None, cv_folds=5, cv_random_state=None, lam_B=0.001, max_sweeps=12, rel_tol=1e-06, d_floor=1e-08, cg_maxit=800, cg_tol=3e-07)[source]

Bases: object

Scikit-learn style estimator for low-rank GLS via ALS.

Parameters:
get_params(deep=True)[source]

Return the estimator hyperparameters (scikit-learn protocol).

Parameters:

deep (bool)

Return type:

dict[str, Any]

set_params(**params)[source]

Set estimator hyperparameters in place (scikit-learn protocol).

Parameters:

params (Any)

Return type:

ALSGLS

als_kwargs()[source]

The keyword arguments this estimator passes to als_gls().

Return type:

dict[str, Any]

fit(Xs, Y)[source]

Fit the GLS system, selecting the rank first when requested.

Parameters:
Return type:

ALSGLS

predict(Xs)[source]

Predict responses for new design matrices, one per equation.

Parameters:

Xs (Sequence[Any])

Return type:

ndarray

score(Xs, Y)[source]

Return the negative Gaussian NLL per row of Y under the fitted Sigma.

Higher is better, as scikit-learn’s score protocol requires. This scores the whole fitted covariance, not just the conditional mean, so it is not the negative MSE; use alsgls.mse for that.

Parameters:
Return type:

float

covariance(method='kh')[source]

Covariance of the coefficients by the named method.

"kh" (default) is the Kackar-Harville corrected covariance, which accounts for Sigma being estimated; "plugin" is the df-rescaled plug-in that linearmodels, systemfit and Stata report. Measured on the test suite’s 4-equation, rank-1 fixture, the reported standard error is 0.85 (plug-in) and 0.89 (kh) of the actual spread at n = 20, and 0.94 / 0.96 at n = 40. Neither is calibrated at n = 20; for that use bootstrap().

Parameters:

method (str) – "kh" or "plugin".

Returns:

The (p_total, p_total) covariance.

Raises:

RuntimeError – If the training designs were not stored.

Return type:

ndarray

property cov_params_: ndarray

The default coefficient covariance, covariance("kh").

bootstrap(B=999, method='parametric', seed=None)[source]

Bootstrap the whole fit and return studentised (bootstrap-t) inference.

Each replicate refits F, D and beta on resampled data, which is what captures the part of the sampling variance the plug-in cov_params_ cannot: that Sigma is estimated. The calibrated object is the percentile-t interval, not the bootstrap standard error; see BootstrapResults.

Parameters:

  • B (int) – Number of replicates. 999 makes 0.05 * (B + 1) an integer, which is what a 95% percentile interval wants; 199 is enough if only the standard errors are of interest.

  • method (str) – "parametric" draws errors from the fitted F F' + diag(D); "wild" multiplies each residual row by a Rademacher sign; "residual" resamples whole residual rows. The last two make no distributional assumption.

  • seed (int | None) – Seed for the replicate stream; the same seed gives the same replicates.

Returns:

The bootstrap distribution and the inference built on it.

Return type:

BootstrapResults

predict_interval(Xs, alpha=0.05, return_type='prediction')[source]

Compute prediction or confidence intervals for new observations.

Parameters:

  • Xs (Sequence[Any]) – List of design matrices [X_0, …, X_{K-1}] for new observations.

  • alpha (float) – Significance level (default 0.05 for 95% intervals).

  • return_type (str) – “prediction” for prediction intervals (includes residual variance), “confidence” for confidence intervals (variance of E[y|X] only).

Returns:

{“mean”: (N, K), “lower”: (N, K), “upper”: (N, K)}

Return type:

dict

Raises:

ValueError – If return_type is neither "prediction" nor "confidence", or the design matrices do not match the fitted model in count or in number of columns.

class alsgls.api.PredictionResults(predicted_mean, se_mean, se_obs, _df, _alpha_default=0.05)[source]

Bases: object

Container for prediction results with intervals.

Parameters:

  • predicted_mean (ndarray)

  • se_mean (ndarray)

  • se_obs (ndarray)

  • _df (int)

  • _alpha_default (float)

predicted_mean: ndarray
se_mean: ndarray
se_obs: ndarray
conf_int_mean(alpha=None)[source]

Confidence intervals for E[y|X].

Parameters:

alpha (float | None) – Significance level (default 0.05 for 95% CI)

Returns:

(N, K, 2) array with lower and upper bounds

Return type:

ci

conf_int_obs(alpha=None)[source]

Prediction intervals for y|X (includes residual variance).

Parameters:

alpha (float | None) – Significance level (default 0.05 for 95% CI)

Returns:

(N, K, 2) array with lower and upper bounds

Return type:

ci

class alsgls.api.ALSGLSSystem(system, *, rank='auto', lam_B=0.001, max_sweeps=12, rel_tol=1e-06, d_floor=1e-08, cg_maxit=800, cg_tol=3e-07)[source]

Bases: object

Statsmodels-style system container for ALS GLS fitting.

Parameters:
property equations: list[_SystemEquation]

The parsed equations, in system order.

property nobs: int

Number of observations shared by every equation.

property keqs: int

Number of equations in the system.

as_arrays()[source]

Return the system as (list of X_j, stacked Y) arrays.

Return type:

tuple[list[ndarray], ndarray]

fit()[source]

Fit the system with ALS-GLS and return a results container.

Return type:

ALSGLSSystemResults

class alsgls.api.ALSGLSSystemResults(model, estimator)[source]

Bases: object

Lightweight results container mimicking statsmodels outputs.

Parameters:
params_as_series()[source]

Return coefficients as a pandas Series with (equation, variable) index.

predict(exog=None)[source]

Predict fitted values, for the training design or new exog.

Parameters:

exog (Mapping[Any, Any] | Sequence[Any] | None)

Return type:

ndarray

summary_dict()[source]

Return the headline fit statistics as a plain dict.

Return type:

dict[str, Any]

covariance(method='kh')[source]

Covariance of the coefficients by the named method.

"kh" (default) is the Kackar-Harville corrected covariance, which accounts for Sigma being estimated; "plugin" is the df-rescaled plug-in that linearmodels, systemfit and Stata report. Measured on the test suite’s 4-equation, rank-1 fixture, the reported standard error is 0.85 (plug-in) and 0.89 (kh) of the actual spread at n = 20, and 0.94 / 0.96 at n = 40. Neither is calibrated at n = 20; for that use bootstrap().

Parameters:

method (str) – "kh" or "plugin".

Returns:

The (p_total, p_total) covariance.

Return type:

ndarray

property cov_params: ndarray

The default coefficient covariance, covariance("kh").

property bse: ndarray

Standard errors of parameter estimates (sqrt of diagonal of cov_params).

property tvalues: ndarray

β = 0.

Type:

t-statistics for H₀

property df_resid: int

nobs * keqs - n_params.

Type:

Residual degrees of freedom

property pvalues: ndarray

β = 0.

Type:

Two-sided p-values for H₀

conf_int(alpha=0.05)[source]

Confidence intervals for parameters.

Parameters:

alpha (float) – Significance level (default 0.05 for 95% CI)

Returns:

(n_params, 2) array with lower and upper bounds

Return type:

ci

get_prediction(exog=None, alpha=0.05)[source]

Get prediction results with standard errors and intervals.

Parameters:

  • exog (Mapping[Any, Any] | Sequence[Any] | None) – Design matrices for new observations. If None, uses the training data. Can be a dict {eq_name: X_j} or a list [X_0, …, X_{K-1}].

  • alpha (float) – Default significance level for intervals (default 0.05).

Returns:

Object with predicted_mean, se_mean, se_obs, and interval methods.

Return type:

PredictionResults

Raises:

ValueError – If design matrices are not supplied for every equation, or one has the wrong number of columns.

bootstrap(B=999, method='parametric', seed=None)[source]

Bootstrap the whole fit and return studentised (bootstrap-t) inference.

Each replicate refits F, D and beta on resampled data, which is what captures the part of the sampling variance the plug-in cov_params cannot: that Sigma is estimated. The calibrated object is the percentile-t interval, not the bootstrap standard error; see BootstrapResults.

Parameters:

  • B (int) – Number of replicates. 999 makes 0.05 * (B + 1) an integer, which is what a 95% percentile interval wants; 199 is enough if only the standard errors are of interest.

  • method (str) – "parametric" draws errors from the fitted F F' + diag(D); "wild" multiplies each residual row by a Rademacher sign; "residual" resamples whole residual rows. The last two make no distributional assumption.

  • seed (int | None) – Seed for the replicate stream; the same seed gives the same replicates.

Returns:

The bootstrap distribution and the inference built on it.

Return type:

BootstrapResults

summary(alpha=0.05)[source]

Text summary of estimation results (statsmodels-style).

Parameters:

alpha (float) – Significance level for confidence intervals

Returns:

Formatted summary table

Return type:

str

alsgls.sim module

Synthetic data generators for GLS/SUR experiments and tests.

alsgls.sim.simulate_sur(N_tr, N_te, K, p, k, seed=0)[source]

Simulate a Seemingly Unrelated Regression (SUR) dataset.

Parameters:

  • N_tr (int) – Number of training samples.

  • N_te (int) – Number of test samples.

  • K (int) – Number of response equations.

  • p (int) – Number of features per equation.

  • k (int) – Latent factor dimension controlling correlated noise.

  • seed (int) – Seed for the NumPy random number generator. Defaults to 0.

Returns:

lists of per-equation feature matrices of shape (N_tr, p) and (N_te, p), and response matrices of shape (N_tr, K) and (N_te, K).

Return type:

A tuple (X_tr, Y_tr, X_te, Y_te)

Notes

Randomness is controlled via numpy.random.default_rng(seed); pass a different seed for different simulations.

Examples

>>> from alsgls import simulate_sur
>>> Xtr, Ytr, Xte, Yte = simulate_sur(100, 20, K=3, p=5, k=2, seed=42)
alsgls.sim.simulate_gls(N_tr, N_te, p_list, k, seed=0)[source]

Simulate a generalized least squares (GLS) dataset.

This variant allows each response equation to have its own number of features as specified by p_list.

Parameters:

  • N_tr (int) – Number of training samples.

  • N_te (int) – Number of test samples.

  • p_list (list[int]) – Number of features for each equation.

  • k (int) – Latent factor dimension controlling correlated noise.

  • seed (int) – Seed for the NumPy random number generator. Defaults to 0.

Returns:

per-equation feature matrices with X_tr[j] of shape (N_tr, p_list[j]) and X_te[j] of shape (N_te, p_list[j]), and responses of shape (N_tr, K) and (N_te, K) where K = len(p_list).

Return type:

A tuple (X_tr, Y_tr, X_te, Y_te)

Notes

Randomness is controlled via numpy.random.default_rng(seed); pass a different seed for different simulations.

Examples

>>> from alsgls import simulate_gls
>>> p_list = [3, 5, 2]
>>> Xtr, Ytr, Xte, Yte = simulate_gls(100, 20, p_list, k=2, seed=0)

alsgls.lsqr_gls module

High-accuracy GLS solves via LSQR/LSMR without squaring the condition number.

This module exposes utilities for solving the weighted least-squares problem

min_beta (y - X beta)^T Σ^{-1} (y - X beta)

in the common “low-rank plus diagonal” setting used across the package. The implementation avoids the explicit normal equations that the in-package CG routine currently relies on and instead wraps the design matrix inside a scipy.sparse.linalg.LinearOperator so that the LSQR/LSMR Krylov solvers can be used directly. In ill-conditioned designs this provides noticeably better convergence and is less sensitive to round-off.

The implementation follows the write-up in the project documentation and is careful to avoid forming dense KxK matrices except for a skinny SVD of the Woodbury core.

class alsgls.lsqr_gls.GLSLinearOperator(*args, **kwargs)[source]

Bases: LinearOperator

Linear operator representing A = W X for LSQR/LSMR.

Parameters:

  • X_dot (Callable[[np.ndarray], np.ndarray])

  • X_Tdot (Callable[[np.ndarray], np.ndarray])

  • W (WoodburyWeight)

  • N (int)

  • K (int)

  • P (int)

class alsgls.lsqr_gls.WoodburyWeight(d, F, d_floor=1e-12, sv_tol=1e-12)[source]

Bases: object

Row-wise operator W satisfying W.T @ W = Σ^{-1}.

Parameters:

  • d (ndarray | Sequence[float]) – Diagonal of D in Σ = D + F F^T.

  • F (ndarray | Sequence[float] | None) – Optional factor loadings. If None or empty then Σ is purely diagonal and the action reduces to simple scaling by D^{-1/2}.

  • d_floor (float) – Lower bound applied element-wise to d to avoid singularities.

  • sv_tol (float) – Relative tolerance used to trim tiny singular values when computing the skinny SVD of U = D^{-1/2} F.

d: ndarray | Sequence[float]
F: ndarray | Sequence[float] | None
d_floor: float = 1e-12
sv_tol: float = 1e-12
W_apply(T)[source]

Apply W to an (N, K) array row-by-row.

Parameters:

T (ndarray)

Return type:

ndarray

WT_apply(T)[source]

Apply the adjoint W.T to an (N, K) array row-by-row.

Parameters:

T (ndarray)

Return type:

ndarray

alsgls.lsqr_gls.make_block_design_ops(X_blocks)[source]

Build X_dot/X_Tdot callbacks for SUR-style block designs.

Parameters:

X_blocks (Sequence[ndarray])

alsgls.lsqr_gls.solve_gls_weighted(X_dot, X_Tdot, y, d, F, *, method='lsmr', atol=1e-10, btol=1e-10, conlim=100000000.0, maxiter=None, verbose=False)[source]

Solve argmin_beta || W (X beta - y) ||_2 via LSQR or LSMR.

The design is provided through matrix-free callbacks X_dot and X_Tdot matching the interfaces used throughout the rest of the alsgls package. The solver works directly with the GLS geometry and therefore avoids squaring the condition number of X.

Parameters:

  • X_dot (Callable[[ndarray], ndarray]) – Callback applying the stacked design to a coefficient vector.

  • X_Tdot (Callable[[ndarray], ndarray]) – Callback applying its transpose to a residual matrix.

  • y (ndarray) – Responses, (N, K).

  • d (ndarray | Sequence[float]) – Diagonal of D in Sigma = D + F F^T.

  • F (ndarray | Sequence[float] | None) – Factor loadings, or None when Sigma is purely diagonal.

  • method (str) – Krylov solver to use, "lsmr" or "lsqr".

  • atol (float) – Absolute stopping tolerance passed to the solver.

  • btol (float) – Relative stopping tolerance passed to the solver.

  • conlim (float) – Condition-number limit at which the solver gives up.

  • maxiter (int | None) – Iteration cap, or None for the solver’s own default.

  • verbose (bool) – Print solver progress.

Returns:

The concatenated coefficient vector. info: Diagnostics returned by the underlying Krylov solver.

Return type:

beta

Raises:

ValueError – If method is neither "lsmr" nor "lsqr".

alsgls.rank_selection module

Rank selection methods for ALS-GLS: BIC and cross-validation.

alsgls.rank_selection.select_rank_bic(Xs, Y, k_candidates=None, **als_kwargs)[source]

Select rank by minimizing BIC (Bayesian Information Criterion).

BIC(k) = 2 * N * nll_per_row + n_params * log(N)

which is the usual -2 * loglik + n_params * log(N), since N * nll_per_row is the negative log-likelihood of the sample. The value used to be reported at half this, so it was not comparable with the bic attribute of statsmodels or any other package; the rank chosen is unchanged, because halving is monotone.

n_params counts the free parameters of the whole fitted model: the K*k factor loadings in F less the k*(k-1)/2 orthogonal rotations F -> F Q that leave F F^T unchanged and so are not identified, the K diagonal variances in D, and the sum(p_j) regression coefficients, since nll_per_row is evaluated at the fitted beta and a BIC has to charge for it. This is the standard factor-analysis count; R’s factanal reports the complementary df = ((K-k)^2 - K - k) / 2.

A rank whose fit raised carries the message in its error key.

The count used to be K*(k+1) + k, which neither subtracted the rotational redundancy nor charged for beta. Only the first term varies with k, so on the fixtures tested the selected rank is unchanged; the reported value was wrong either way.

Parameters:

  • Xs (list[ndarray]) – Design matrices for each equation.

  • Y (ndarray) – Response matrix.

  • k_candidates (list[int] | None) – Candidate ranks to evaluate. Defaults to range(1, min(K//2, 12)+1).

  • **als_kwargs (Any) – Additional arguments passed to als_gls().

Returns:

The rank with minimum BIC. results: Per-rank dicts of ‘k’, ‘nll’, ‘bic’, ‘n_params’, ‘converged’.

Return type:

best_k

Raises:

RuntimeError – If no candidate rank produced a usable fit.

alsgls.rank_selection.select_rank_cv(Xs, Y, k_candidates=None, n_folds=5, random_state=None, **als_kwargs)[source]

Select rank by k-fold cross-validation on validation NLL.

Parameters:

  • Xs (list[ndarray]) – Design matrices for each equation.

  • Y (ndarray) – Response matrix.

  • k_candidates (list[int] | None) – Candidate ranks to evaluate. Defaults to range(1, min(K//2, 12)+1).

  • n_folds (int) – Number of cross-validation folds.

  • random_state (int | Generator | None) – Random state for reproducible fold splits.

  • **als_kwargs (Any) – Additional arguments passed to als_gls().

Returns:

The rank with minimum mean CV NLL. results: Per-rank results containing ‘k’, ‘cv_nll’, ‘cv_std’, ‘fold_nlls’.

Return type:

best_k

Raises:

  • ValueError – If n_folds is below 2 or exceeds the number of rows.

  • RuntimeError – If no candidate rank produced a usable fit.

alsgls.metrics module

Fit metrics for low-rank-plus-diagonal GLS models.

alsgls.metrics.mse(Y, Yhat)[source]

Mean squared error between two matrices.

Parameters:

  • Y (ndarray)

  • Yhat (ndarray)

Return type:

float

alsgls.metrics.nll_per_row(R, F, D)[source]

Negative log-likelihood per row for a residual matrix.

Computed for R under Σ = F F^T + diag(D) with Gaussian errors.

Parameters:

  • R (ndarray) – Residual matrix, (N, K).

  • F (ndarray) – Factor loadings, (K, k).

  • D (ndarray) – Diagonal variances, length K.

Returns:

0.5 * [ tr(R Σ^{-1} R^T)/N + log det(Σ) + K log(2π) ] where N is the number of rows in R.

Return type:

float