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.
A covariance matrix is not just a heatmap. Every cell summarizes paired deviations from estimated means, and every portfolio-risk calculation relies on those cells fitting together coherently.
This tutorial builds the sample matrix from complete aligned observations, distinguishes the n−1 denominator from maximum-likelihood covariance, and explains why a singular matrix can be mathematically valid but operationally inconvenient.
Open this figure at full size.
Rows are observations; columns are assets
Let X contain n observations of p assets, with each row covering the same interval for every asset. Compute each column mean, then center the rows: z_t=x_t−x̄. The sample covariance is S=Σz_tz_tᵀ/(n−1).
The diagonal S_ii is an asset's sample variance. Off-diagonal S_ij measures the average product of centered movements, with the sample correction. Negative covariance is valid. It is not a negative variance and not a claim of causal influence.
Maximum-likelihood empirical covariance uses denominator n instead. The scikit-learn empirical-covariance reference documents that convention. This topic deliberately uses n−1. The shrinkage topics use ML covariance, so a comparison must account for the scale difference before attributing it to regularization.
Construct a matrix with a known answer
Take four centered two-asset observations. Let the first column be [1,1,−1,−1]. Let the second be0.5 times that column plus √0.75 times [1,−1,1,−1]. The two underlying sign vectors are orthogonal and each has squared length four.
Both column means are zero. The sums of squared values are four for both assets; the cross-product sum is two. Dividing by n−1=3 gives
This intentionally large synthetic scale keeps the arithmetic simple. Multiplying every observation by 0.01 gives decimal-return-sized inputs and multiplies the matrix by 0.0001. Dividing by n instead would yield [[1,0.5],[0.5,1]]. Neither denominator should be hidden.
Open this figure at full size.
Download the exact worked input and expected values.
A portfolio check ties the cells together
For portfolio weights a, sample variance is aᵀSa. With a=[0.5,0.5] in the worked example, the result is 1. Compute the four portfolio observations directly and take their sample variance: the answer must match.
This identity checks more than symmetry. It confirms that the matrix corresponds to the same centered observation panel. A covariance matrix built from pairwise samples with different missing rows may not retain that coherent outer-product representation.
Positive semidefiniteness follows algebraically: aᵀSa=Σ(aᵀz_t)²/(n−1)≥0. Zero is allowed for some nonzero a. That is why “PSD” and “invertible” are different requirements.
Why dimensions can defeat inversion
Centering n observations leaves rank at most n−1. If p≥n, the p×p sample covariance cannot have full rank. Perfectly collinear assets can make it singular even when n is much larger than p.
A singular matrix is not necessarily a data error. It means the sample cannot identify independent variation in every portfolio direction. Inverting it blindly can fail; nearly singular matrices can amplify small estimation errors into large precision-matrix or optimizer changes.
| Need | Appropriate response |
|---|---|
| Explain observed co-movement | Inspect sample covariance and its data panel |
| Compare strength across differently scaled assets | Normalize to correlation, with zero-variance handling |
| Stabilize a noisy high-dimensional estimate | Evaluate shrinkage |
| Attribute common versus asset-specific risk | Specify a factor model |
| Infer sparse conditional associations | Consider a precision-model approach with assumptions |
Open this figure at full size.
Use the playground to connect data and geometry
The lab builds matrices from prefixes of a 64-row synthetic four-asset panel. Step adds a complete observation row, not one asset in isolation. The matrix cells and displayed sample count update together.
Change the dependence structure through the comparison control, then inspect the off-diagonal cells. The collinear edge case demonstrates that a coherent sample can be singular. A malformed rectangular panel is rejected; missing entries are not silently replaced with zero.
The small worked fixture also supports the portfolio check above. Predict how multiplying one asset's returns by two affects its variance and covariances: its diagonal quadruples, its cross-covariances double, and other cells stay unchanged.
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 the reference guarantees
It accepts a finite rectangular n×p matrix with n≥2 and p≥2. Rows are assumed to have been aligned upstream; the array alone does not verify timestamps, identifiers or adjustment histories. Asset order must remain stable. Mixing percent and decimal-return columns changes the economic meaning of the matrix.
Means are estimated from the supplied endpoint sample. That is appropriate for an endpoint estimate. It is not permission to use full-dataset means for every earlier backtest date. Recompute on each available prefix or use a separately declared causal mean model.
The direct implementation takes O(np²) time and returns the estimated means and denominator. Tests verify exact matrix cells, scaling, valid zero matrices, malformed inputs and cross-language equality. A heatmap's appearance is never used as a substitute for numeric checks.
Start here before Ledoit–Wolf shrinkage. Understanding which matrix you are shrinking is the prerequisite for understanding what the target and intensity actually change.
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 sample_covariance import calculate
data = json.loads(Path("teaching-path.json").read_text())
result = calculate(data)
print(result["latest"])
import {calculate} from './sample_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.
Complete aligned n by p observations, sample denominator n-1. Singular PSD is legitimate; no implicit inverse.
Continue the investigation
- EWMA Covariance: compare its assumptions and information boundary before comparing the numbers.
- Graphical-Lasso Covariance: compare its assumptions and information boundary before comparing the numbers.
Rendered from the canonical Mermaid sources linked by this article.
Sample Covariance — calculation-flow
Sample 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 — empirical_covariance 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.
Scope of evidence
Complete aligned n by p observations, sample denominator n-1. Singular PSD is legitimate; no implicit inverse.
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-A01 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,"sample_covariance");
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-A01",title:"Sample Covariance",parameters:p,...result};
}
The embedded lab now expands to its full document height, keeping the article as the only scroll surface.