By Islam BarakaFor analysts and developers who know aligned return vectors, matrix multiplication, and covariance versus correlation. The goal is to reproduce the mechanism, inspect its failure states, and decide what the output can legitimately tell you—not to fit or endorse a trading strategy.
Two assets can move together because both share a relationship with a third. Graphical lasso looks for a sparse structure in the precision matrix, the inverse of covariance, to describe conditional associations under a Gaussian model.
The distinction matters. Penalizing covariance entries directly is a different problem. Drawing an edge because two raw returns are correlated is a different graph. This tutorial makes the optimization objective, matrix roles and convergence evidence visible.
Open this figure at full size.
Fix the objective—and say which entries are penalized
Let C be the centered empirical covariance with denominator n. The precision estimate Θ solves
Only off-diagonal precision entries are penalized. The graphical-lasso API reference documents the sparse inverse-covariance problem and numerical solver controls. Penalizing the diagonal as well would change the solution; the norm shorthand must not obscure this convention.
The covariance estimate is W=Θ⁻¹. A zero Θ_ij corresponds to conditional independence for a jointly Gaussian distribution, conditional on the other included variables. Outside that model, it should not be promoted to a universal independence statement. It never establishes causation.
Solve a two-asset example exactly
Take C=[[1,0.5],[0.5,1]] and α=0.1. For this two-asset off-diagonal-penalty problem, the covariance diagonals stay at one and the off-diagonal is soft-thresholded to 0.4:
Thus Θ_11≈1.19047619 and Θ_12≈−0.476190476. The partial-correlation normalization −Θ_12/√(Θ_11Θ_22) is 0.4. The precision off-diagonal has the opposite sign from this partial correlation; a graph legend must say which quantity it displays.
At α≥0.5, the off-diagonal covariance and precision entries are zero in this special two-variable example. That elementwise soft-threshold rule is not a general solution for larger matrices.
Open this figure at full size.
Download the exact worked input and expected values.
The diagonal is a powerful regression check
Differentiate the smooth objective with respect to an unpenalized diagonal precision entry. The optimality condition is C_ii−W_ii=0. Therefore the fitted covariance diagonal equals the empirical covariance diagonal for this objective.
A solver that adds α or an extra quadratic term to W_ii is solving something else or is wrong. An earlier implementation in this package violated this condition while still reporting convergence. The repaired version keeps the diagonal fixed and checks optimality independently of iteration-to-iteration movement.
For nonzero off-diagonal precision entries, the KKT condition is C_ij−W_ij+α sign(Θ_ij)=0. At a zero entry, |C_ij−W_ij|≤α. These conditions explain what “converged” should mean numerically, not merely that a displayed matrix stopped changing.
Read the diagnostics before reading the graph
The teaching solver uses block coordinate descent with inner lasso updates. It checks positive definiteness, a scale-normalized KKT residual, and the dual gap. It reports iteration history and rejects an exhausted iteration budget instead of returning an unverified matrix as a successful fit.
| Diagnostic | What it checks | What it cannot prove |
|---|---|---|
| Positive-definite covariance | Valid invertible numerical state | Correct optimum by itself |
| KKT residual | Objective's optimality conditions | Scientific suitability of the model |
| Dual gap | Primal/dual consistency near optimum | Stability of selected edges across samples |
| Edge count | Sparsity at the selected penalty | Causal structure |
The solver is restricted to 2–12 assets for transparent teaching. It is not a replacement for a production-scale, well-tested optimization library. Very ill-conditioned inputs can fail numerical checks; that failure should remain visible.
Open this figure at full size.
Penalty scale is not dimensionless
α has covariance-scale units in this objective. If all residuals are multiplied by c, their covariance scales by c²; scaling α by c² preserves the corresponding solution structure after rescaling. Moving from decimals to percentages without adjusting α dramatically changes regularization.
Standardizing each asset before fitting is another choice, not a cosmetic chart operation. It changes the effective penalty relative to asset scales. Record whether the model used raw covariance or standardized variables and how any result was transformed back.
Use the playground as a penalty experiment
Hold the synthetic four-asset sample fixed and vary α separately from the observation count. Inspect C, W, Θ and the partial-correlation graph. Some edges may vanish as the penalty increases; that does not imply that the corresponding marginal covariance is necessarily zero in a larger system.
Step adds one complete observation row and refits at the same penalty. This reveals sample sensitivity. Solver history is displayed as optimization evidence, not confused with market time. The zero-variance-asset failure shows why a degenerate input is rejected by this particular teaching solver.
Open this figure at full size.
Open the standalone guided playground. The embedded playground and runnable code are available on this page. Download the 64-observation teaching input.
What makes the graph useful—and what would make it misleading
A sparse precision graph can summarize conditional association structure and support a regularized covariance estimate. Its meaning depends on the included variables, distributional assumptions, sample and selected penalty. An omitted common driver can change the apparent graph.
Choosing α, assessing edge stability and evaluating held-out likelihood or risk behavior require an explicit research protocol. The package supplies coefficients and synthetic data; it does not claim an economically optimal penalty or discovered market network.
Independent tests compare the fitted covariance and precision with scikit-learn across several dimensions, plus the exact two-variable solution above. That evidence establishes a numerical implementation within its stated scope. The reader's final task is to keep that numerical claim separate from the much stronger scientific interpretation of an edge.
Reproduce and inspect the calculation
The Python and TypeScript tabs contain standalone implementations, not imports into an unseen runtime. Both expose calculate(input_data). Feed the worked JSON's input object into that entry point. For the longer experiment, use the teaching-path JSON directly.
import json
from pathlib import Path
from graphical_lasso_covariance import calculate
data = json.loads(Path("teaching-path.json").read_text())
result = calculate(data)
print(result["latest"])
import {calculate} from './graphical_lasso_covariance.ts';
const result = calculate(inputData); // inputData is the downloaded JSON object
console.log(result.latest);
Place the downloaded input beside your script and the standalone source on its import path. The Python reference uses the standard library; the TypeScript reference has no external runtime dependency. Shared tests include independent numeric anchors, valid boundaries, rejected inputs and cross-language output comparisons. They establish arithmetic, not forecasting performance.
Evidence and scope
This article uses authored synthetic calculations and primary technical references, reviewed 2026-09-10. Historical market examples are deferred until identity, adjustment basis, chronology and redistribution rights can be verified. No personal trading history or search-ranking superiority is asserted.
ML covariance; off-diagonal precision penalty only. Fixed covariance diagonal C_ii, positive alpha, <=12 assets teaching solver. Proximal KKT residual, SPD and iteration budget checked; non-convergence is an error.
- scikit-learn — graphical_lasso API and implementation notes
- Friedman, Hastie and Tibshirani — Sparse Inverse Covariance Estimation with the Graphical Lasso
Continue the investigation
- Sample Covariance: compare its assumptions and information boundary before comparing the numbers.
- Factor-Model Covariance: compare its assumptions and information boundary before comparing the numbers.
Rendered from the canonical Mermaid sources linked by this article.
Graphical-Lasso Covariance — calculation-flow
Graphical-Lasso Covariance — decision-boundary
ReferencesPrimary sources and evidence notesExpand the source trail, evidence role, and limitations behind the engineering choices.
Expand the source trail, evidence role, and limitations behind the engineering choices.
Reviewed 2026-09-10. Primary technical documentation and papers; synthetic arithmetic is author-derived. This is a targeted source review, not a verified review of Google's top ten results and not a claim of ranking superiority.
- S1: scikit-learn — graphical_lasso API and implementation notes — accessed 2026-09-10. Rolling official documentation snapshot; exact package versions used for numerical comparisons are recorded in the repair numeric-evidence.json. Supports the definition and declared convention, not investment performance. Jurisdiction: not applicable to this mathematical reference.
- S2: Friedman, Hastie and Tibshirani — Sparse Inverse Covariance Estimation with the Graphical Lasso — accessed 2026-09-10. Original research publication; not a current market observation. Supports the definition and declared convention, not investment performance. Jurisdiction: not applicable to this mathematical reference.
Convention boundary: the original paper uses an all-entry penalty in its displayed objective. This package follows the off-diagonal-only software convention and verifies the fixed covariance diagonal; do not silently equate the two objectives.
Scope of evidence
ML covariance; off-diagonal precision penalty only. Fixed covariance diagonal C_ii, positive alpha, <=12 assets teaching solver. Proximal KKT residual, SPD and iteration budget checked; non-convergence is an error.
Historical case: deferred. No public provider dataset, historical performance claim, or personal trading anecdote is used. Synthetic examples demonstrate arithmetic, not market efficacy. Sources are not copied as article prose.
Accessed: 2026-09-10.
Supports: estimator definition and the explicitly declared variants.
Limitations: technical documentation does not verify a real market feed, author experience, forecast efficacy or search-result superiority. Original-paper access limitations are recorded in the repair report.
Full dependency-light reference implementations in both supported languages.
/** Standalone D10-F04-A06 reference. Generated from validated D10 v2 source. */
export class ContractError extends Error {}
type RecordValue=Record<string, any>;
type Matrix=number[][];
const sum=(x:number[]):number=>x.reduce((a,b)=>a+b,0);
function requireValue(ok: unknown, code: string, message: string): asserts ok {
if (!ok) throw new ContractError(`${code}: ${message}`);
}
function finite(x: unknown, name: string): number {
requireValue(typeof x === 'number' && Number.isFinite(x), 'NUMBER', `${name} must be a finite number`);
return x;
}
function integer(x: unknown, name: string, minimum = 0, maximum = 10000): number {
const v = finite(x, name);
requireValue(Number.isInteger(v) && v >= minimum && v <= maximum, 'INTEGER', `${name} must be an integer in [${minimum}, ${maximum}]`);
return v;
}
function param(p: RecordValue, key: string, fallback: number): number {
return finite(Object.hasOwn(p, key) ? p[key] : fallback, key);
}
function option(p: RecordValue, key: string, fallback: any): any {
return Object.hasOwn(p, key) ? p[key] : fallback;
}
function positive(p: RecordValue, key: string, fallback: number): number {
const v = param(p, key, fallback);
requireValue(v > 0, 'RANGE', `${key} must be positive`);
return v;
}
function vector(value: unknown, name: string, minimum = 1): number[] {
requireValue(Array.isArray(value) && value.length >= minimum, 'SHAPE', `${name} needs ${minimum} or more values`);
return value.map((v, i) => finite(v, `${name}[${i}]`));
}
function matrix(value: unknown, name: string, minRows = 1, minCols = 1): Matrix {
requireValue(Array.isArray(value) && value.length >= minRows, 'SHAPE', `${name}: too few rows`);
const rows = value.map(row => vector(row, name, minCols));
requireValue(rows.every(row => row.length === rows[0].length), 'SHAPE', `${name} must be rectangular`);
return rows;
}
function eye(n: number, value = 1): Matrix {
return Array.from({length: n}, (_, i) => Array.from({length: n}, (_, j) => i === j ? value : 0));
}
function psd(a: Matrix, name: string, strict = false): number | null {
const n = a.length;
requireValue(a.every(row => row.length === n), 'SHAPE', `${name} must be square`);
const tolerance = 1e-12 * Math.max(...a.flat().map(Math.abs), 1e-300);
requireValue(a.every((row, i) => row.every((v, j) => Math.abs(v - a[j][i]) <= tolerance)), 'PSD', `${name} must be symmetric`);
const lower = eye(n), pivots: number[] = [];
for (let j = 0; j < n; j++) {
const pivot = a[j][j] - sum(pivots.map((v, k) => lower[j][k] ** 2 * v));
requireValue(strict ? pivot > tolerance : pivot >= -tolerance, 'PSD', `${name} must be positive semidefinite (strict when requested)`);
pivots.push(pivot > tolerance ? pivot : 0);
for (let i = j + 1; i < n; i++) {
const residual = a[i][j] - sum(pivots.slice(0, j).map((v, k) => lower[i][k] * lower[j][k] * v));
requireValue(pivots[j] !== 0 || Math.abs(residual) <= tolerance, 'PSD', `${name} has invalid zero pivot`);
lower[i][j] = pivots[j] ? residual / pivots[j] : 0;
}
}
return strict ? sum(pivots.map(Math.log)) : null;
}
function inverse(a: Matrix): Matrix {
const n = a.length, identity = eye(n), work = a.map((row, i) => [...row, ...identity[i]]);
for (let j = 0; j < n; j++) {
let k = j;
for (let i = j + 1; i < n; i++) if (Math.abs(work[i][j]) > Math.abs(work[k][j])) k = i;
[work[j], work[k]] = [work[k], work[j]];
requireValue(work[j][j] !== 0, 'NUMERIC', 'singular inverse');
const pivot = work[j][j];
work[j] = work[j].map(v => v / pivot);
for (let i = 0; i < n; i++) if (i !== j) {
const factor = work[i][j];
work[i] = work[i].map((v, l) => v - factor * work[j][l]);
}
}
return work.map(row => row.slice(n));
}
function glassoCovariance(S: Matrix, alpha: number, maxIter = 200, tol = 1e-8): [Matrix, RecordValue] {
const n = S.length;
requireValue(n >= 2 && n <= 12 && alpha > 0 && Math.min(...S.map((row, i) => row[i])) > 0, 'RANGE', 'glasso needs 2..12 nonconstant assets and positive alpha');
const scale = Math.max(...S.map((row, i) => row[i])), C = S.map(row => row.map(v => v / scale)), penalty = alpha / scale;
// Dual-feasible SPD initialization; a diagonal start is not generally feasible.
const maximumOff = Math.max(...C.flatMap((row, i) => row.filter((_, j) => i !== j).map(Math.abs)));
const blend = maximumOff ? Math.min(.1, penalty / maximumOff) : .1;
const W = C.map((row, i) => row.map((v, j) => i === j ? v : (1 - blend) * v)), history: RecordValue[] = [];
for (let iteration = 0; iteration < maxIter; iteration++) {
for (let j = 0; j < n; j++) {
const indices = Array.from({length: n}, (_, i) => i).filter(i => i !== j), beta = indices.map(() => 0);
for (let inner = 0; inner < 2000; inner++) {
let change = 0;
indices.forEach((i, k) => {
const partial = C[i][j] - sum(indices.map((q, l) => l === k ? 0 : W[i][q] * beta[l]));
const updated = Math.sign(partial) * Math.max(Math.abs(partial) - penalty, 0) / W[i][i];
change = Math.max(change, Math.abs(updated - beta[k])); beta[k] = updated;
});
if (change < tol * 1e-5) break;
}
const values = indices.map(i => sum(indices.map((q, l) => W[i][q] * beta[l])));
indices.forEach((i, k) => { W[i][j] = W[j][i] = values[k]; });
}
const logdet = psd(W, 'glasso covariance', true)!, precision = inverse(W);
let residual = 0, trace = 0, l1 = 0;
for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) {
const gradient = C[i][j] - W[i][j];
const trial = precision[i][j] - gradient;
const prox = Math.sign(trial) * Math.max(Math.abs(trial) - penalty, 0);
const error = i === j ? Math.abs(gradient) : Math.abs(precision[i][j] - prox);
residual = Math.max(residual, error); trace += C[i][j] * precision[j][i];
if (i !== j) l1 += penalty * Math.abs(precision[i][j]);
}
const gap = trace - n + l1;
history.push({iteration: iteration + 1, kkt_residual: residual, dual_gap: gap, objective: trace + logdet + l1 + n * Math.log(scale)});
if (residual <= tol && Math.abs(gap) <= tol * Math.max(1, n)) return [W.map(row => row.map(v => v * scale)),
{iterations: iteration + 1, converged: true, alpha, precision: precision.map(row => row.map(v => v / scale)), kkt_residual: residual, dual_gap: gap, history}];
}
throw new ContractError('CONVERGENCE: graphical lasso did not satisfy KKT and dual-gap tolerances; iteration budget exhausted');
}
function covariance(data: RecordValue, p: RecordValue, kind: string): RecordValue {
if (kind === 'factor_covariance') {
const B = matrix(data.loadings, 'loadings', 2), F = matrix(data.factor_covariance, 'factor_covariance'), D = vector(data.specific_variances, 'specific_variances');
const assets = B.length, factors = B[0].length;
requireValue(F.length === factors && D.length === assets, 'SHAPE', 'factor dimensions do not match');
psd(F, 'factor_covariance'); requireValue(Math.min(...D) >= 0, 'RANGE', 'specific variances must be nonnegative');
const common = B.map(row => B.map(other => sum(row.map((v, k) => sum(other.map((u, l) => v * F[k][l] * u))))));
const result = common.map((row, i) => row.map((v, j) => v + (i === j ? D[i] : 0)));
return {matrix: result, latest: result, ready: true, ready_at: 0, diagnostics: {common, specific_variances: D, assets, factors, causal: true}};
}
const X = matrix(data.returns, 'returns', 2, 2), n = X.length, dim = X[0].length;
const means = X[0].map((_, j) => sum(X.map(row => row[j])) / n), Z = X.map(row => row.map((v, j) => v - means[j]));
const C = X[0].map((_, i) => X[0].map((_, j) => sum(Z.map(z => z[i] * z[j])) / n));
const diagnostics: RecordValue = {observations: n, assets: dim, means, causal: true, ml_covariance: C};
let result: Matrix;
if (kind === 'sample_covariance') { result = C.map(row => row.map(v => v * n / (n - 1))); diagnostics.denominator = n - 1; }
else if (kind === 'ewma_covariance') {
const decay = param(p, 'decay', .94);
requireValue(decay > 0 && decay < 1, 'RANGE', 'decay must be in (0,1)');
result = Object.hasOwn(data, 'initial_covariance') ? matrix(data.initial_covariance, 'initial_covariance') : eye(dim, 0);
requireValue(result.length === dim, 'SHAPE', 'initial covariance dimensions mismatch'); psd(result, 'initial_covariance');
const states: Matrix[] = [];
for (const row of X) { result = result.map((old, i) => old.map((v, j) => decay * v + (1 - decay) * row[i] * row[j])); states.push(result); }
Object.assign(diagnostics, {states, mean_convention: 'supplied zero-mean residuals', weight_mass: 1 - decay ** n, seed_weight: decay ** n});
} else if (kind === 'graphical_lasso') {
const solver = glassoCovariance(C, positive(p, 'alpha', .00002), integer(option(p, 'max_iter', 200), 'max_iter', 1, 2000), positive(p, 'tolerance', 1e-8));
result = solver[0]; diagnostics.graphical_lasso = solver[1];
} else {
const trace = sum(C.map((row, i) => row[i])), mu = trace / dim, tr2 = sum(C.flat().map(v => v * v));
const delta = sum(C.flatMap((row, i) => row.map((v, j) => (v - (i === j ? mu : 0)) ** 2)));
let rho: number;
if (kind === 'ledoit_wolf') {
const noise = sum(Z.map(z => sum(C.flatMap((row, i) => row.map((v, j) => (z[i] * z[j] - v) ** 2))))) / n ** 2;
rho = delta > 0 ? Math.min(1, Math.max(0, noise / delta)) : 0;
Object.assign(diagnostics, {noise_estimate: noise, target_distance: delta});
} else {
const numerator = (1 - 2 / dim) * tr2 + trace * trace, denominator = (n + 1 - 2 / dim) * (tr2 - trace * trace / dim);
rho = tr2 === 0 ? 0 : denominator > 1e-14 * tr2 ? Math.min(1, Math.max(0, numerator / denominator)) : 1;
Object.assign(diagnostics, {numerator, denominator, variant: 'original finite-p OAS'});
}
result = C.map((row, i) => row.map((v, j) => (1 - rho) * v + (i === j ? rho * mu : 0)));
Object.assign(diagnostics, {shrinkage: rho, target: mu});
}
return {matrix: result, latest: result, ready: true, ready_at: 0, diagnostics};
}
export function calculate(data: RecordValue): RecordValue {
requireValue(data && typeof data === 'object' && !Array.isArray(data),'SHAPE','input must be an object');
const p=Object.hasOwn(data,'parameters')?data.parameters:{};
requireValue(p && typeof p === 'object' && !Array.isArray(p),'SHAPE','parameters must be an object');
const result=covariance(data,p,"graphical_lasso");
function check(v:any):void {
if(typeof v==='number')requireValue(Number.isFinite(v),'NUMERIC','nonfinite computed output');
else if(Array.isArray(v))v.forEach(check);
else if(v && typeof v==='object')Object.values(v).forEach(check);
}
check(result);
return {topic_id:"D10-F04-A06",title:"Graphical-Lasso Covariance",parameters:p,...result};
}
The embedded lab now expands to its full document height, keeping the article as the only scroll surface.