D14-F01-A02 / Released engineering topic

Global Minimum Variance: Why Clipping Misses the Optimum

A production-minded guide to Global Minimum Variance.

9 min read Reviewed Sep 16, 2026 Updated Sep 16, 2026 Python + TypeScript
Global Minimum Variance: Why Clipping Misses the OptimumD14 / D14-F01

If expected returns are too uncertain to use as a target, you can ask a narrower question: which feasible combination has the smallest variance under the covariance estimate? Diversification makes the answer depend on how assets move together, so choosing the least volatile asset alone need not solve the problem.

This tutorial freezes one narrow core so those boundaries are visible. Given a finite covariance matrix in one declared currency and horizon, find the fully-invested portfolio with the smallest variance, subject to non-negative weights. The core does not reconstruct history from a holdings snapshot and it does not repair corporate actions silently.

Before you start

This tutorial is for developers comfortable with weighted averages, vectors, and basic covariance. By the end you can derive a two-asset GMV allocation, explain why clipping fails, and inspect an optimality gap. All percentages use one common horizon; .05 means 5%, and variance has squared-return units. No live market feed or numerical library is required to run the package examples.

What “global” means here

The feasible set is the long-only simplex:

Δ={w∈RN:wi≥0,  1Tw=1}.\Delta = \{w \in \mathbb{R}^N : w_i \ge 0,\;\mathbf{1}^{\mathsf T}w = 1\}.

The frozen objective is:

min⁡w∈Δq(w),q(w)=wTΣw.\min_{w \in \Delta} q(w), \qquad q(w)=w^{\mathsf T}\Sigma w.

This objective and its budget/long-only constraint are the declared Markowitz variant used here; compare the formulation in the official MOSEK Portfolio Optimization Cookbook.

“Global” means global over this exact feasible set. It does not mean best, safest, or guaranteed to outperform another allocation. A target return, shorting, leverage, turnover, transaction costs and a risk-free asset belong to different variants and are not smuggled into this one.

mu may be supplied for a reported expected-return diagnostic, but it does not change the GMV answer. Return horizon, annualization and base currency must be declared by the caller; the core does not infer them.

A small calculation we can audit by hand

Take two assets and the synthetic covariance matrix:

Plain text
Sigma = [[0.04, 0.01],
         [0.01, 0.09]]

With w_B = 1 - w_A, the stationary solution is:

wA=0.09−0.010.04+0.09−2(0.01)=811,wB=311.w_A = \frac{0.09-0.01}{0.04+0.09-2(0.01)} = \frac{8}{11}, \qquad w_B=\frac{3}{11}.

Both weights satisfy the long-only constraints. Direct multiplication gives:

q(w)=(811)2(0.04)+2(811)(311)(0.01)+(311)2(0.09)=7220.q(w)=\left(\frac{8}{11}\right)^2(0.04)+2\left(\frac{8}{11}\right)\left(\frac{3}{11}\right)(0.01)+\left(\frac{3}{11}\right)^2(0.09)=\frac{7}{220}.

Thus variance is 0.03181818181818182 and volatility is 0.17837651700316892. Equal weights have variance 0.0375; that is an arithmetic contrast for this fixture, not an investment-performance claim.

The production result also carries primal residuals and a certificate. For g = 2 Sigma w, the simplex Frank–Wolfe gap is:

gTw−min⁡igi.g^{\mathsf T}w - \min_i g_i.

For a positive-semidefinite covariance, a feasible result with a sufficiently small gap has a convex global-optimality certificate. A small change in the last few weights is not such a certificate.

Why clipping a shorting answer is not optimization

The unconstrained stationary solution is useful as a diagnostic, but it is not the answer after the feasible set changes. Consider:

Plain text
Sigma = [[1, 0.85, 0.5],
         [0.85, 1, 0.1],
         [0.5, 0.1, 1]]

The unconstrained fully-invested vector is [-18/19, 22/19, 15/19]. A tempting shortcut clips its first coordinate and renormalizes the rest, producing [0,22/37,15/37] with variance 775/1369 ~= 0.5661066472.

The constrained solver instead returns [0,0.5,0.5], whose variance is 11/20 = 0.55 and whose FW gap is zero. The active constraint changes the problem; it is not a post-processing detail.

Variance comparison for direct long-only optimization and clipping

Figure: both displayed portfolios are feasible. The direct solution has lower variance. Bars share a common zero-based scale; the unconstrained shorting vector is explained above.

Open the full-size variance comparison

Validate the economics before the optimizer

The covariance core accepts moments only after an upstream audit. That audit should preserve the following distinctions:

CheckFailure that must remain visible
Corporate actionA split-adjusted per-share mark paired with unadjusted quantity changes apparent wealth. Adjust both price and shares under the event convention; the SEC stock-split explanation describes the continuity problem.
Valuation basisAcquisition cost is not a current valuation mark. Do not build return history from a cost field.
CurrencyRecord local units, FX quote direction, FX date, cutoff/staleness rule, and one conversion to the declared base currency.
DividendCash/receivable belongs in the wealth ledger. An adjusted total-return series already incorporates the price/cash treatment; adding the cash again double counts it. See Investor.gov's ex-dividend explanation.
HistoryA point-in-time holdings snapshot is not a return series. Preserve observation and availability clocks, revisions and missing data.

An isolated ticker/quantity/price lead without an observation date, price basis and quantity history is not enough to reconstruct a return. It is excluded from this tutorial; the examples below are synthetic and reproducible.

A stale observation changes the covariance problem

The upstream Stale-Quote Detector keeps source-event age, transport time, unchanged duration and session heartbeat separate. A flat return sequence alone does not prove that a feed is healthy or broken. The optimizer should receive a covariance only after the missing/stale decision is visible in the input audit.

Here is a small synthetic consequence that can be recomputed without a market claim. Aligned returns are

Plain text
A = [-.02, -.01, .01, .02]
B = [-.01,  .01, -.01, .01]

Using the sample (n - 1) covariance gives [[1/3000, 1/15000], [1/15000, 1/7500]]. Its GMV weights are [1/5, 4/5] and its variance is 3/25000 = 0.00012. If the missing A observations are instead replaced by zero returns, the displayed covariance becomes diag(0, 1/7500). The optimizer then returns [1, 0] with variance zero. That is a data-quality diagnostic caused by the substitution, not evidence of a riskless investment. The clock/session evidence and the imputation policy must therefore travel with the covariance.

Where the optimizer fits in the input chain

The related Backward Split Adjustment and Cash-Dividend Total-Return Adjustment topics are input transformations and ledger checks. They do not replace the GMV objective, and this tutorial does not rank one method above another. Keep the transformations' conventions explicit, then pass one ordered, same-horizon moment matrix to the core.

Numerical policy in the reference implementation

The implementation normalizes by the actual positive scale s = max(abs(Sigma[i][j])). It does not use max(1,s), because that would hide a tiny but materially indefinite daily covariance. Symmetry is checked in normalized units and only a within-tolerance matrix is symmetrized. Materially indefinite or non-finite input is rejected; eigenvalues are not clipped and no ridge is added.

The solver starts at equal weights, takes projected-gradient steps, projects exactly onto the simplex, and line-searches the quadratic. If floating-point cancellation stalls the direction, an active-set KKT polish solves each candidate face without regularization. The final result still has to pass the scale-aware FW gap and feasibility tolerances. Singular does not automatically mean non-unique: diag(0,1) has a unique zero-risk simplex solution, whereas an all-ones covariance is flat on the simplex. A strictly positive but near-flat representable curvature is reported as unknown when uniqueness cannot be proved.

The implementation emits optimal, invalid_input or numerical_issue and includes raw scale, eigenvalue/residual, objective, gap and tolerance diagnostics. A covariance accepted only in the small negative-eigenvalue roundoff band receives a tolerance-limited certificate, not an exact convex proof.

Run the same fixture in both languages

Python:

Python
from implementations.python.gmv import global_minimum_variance

result = global_minimum_variance(
    ["A", "B"],
    [[0.04, 0.01], [0.01, 0.09]],
)
print(result["weights"], result["variance"], result["fwGapS"])

TypeScript:

TypeScript
import {globalMinimumVariance} from "../implementations/typescript/gmv.js";

const result = globalMinimumVariance(
  ["A", "B"],
  [[0.04, 0.01], [0.01, 0.09]],
);
console.log(result.weights, result.variance, result.fwGapS);

Independent tests calculate the canonical result, the three-asset clipping failure, boundary negative-weight comparators, singular unique and non-unique cases, all-zero N=1 uniqueness, near-flat classification, malformed inputs, and zero-iteration non-certification. Python and TypeScript share only the input fixture and expected arithmetic; neither calls the other language as its oracle.

Reproduce, then challenge the answer

From this topic's root folder, run the dependency-free Python example and tests:

Shell
python -B examples/run.py
python -B -m unittest discover -s tests -p "test_*.py"

The earlier example and the packaged runnable example use the same canonical inputs.

For TypeScript, the runnable entry point imports the canonical core. The commands below use Node.js 22 and this repository’s TypeScript 7 compiler. Compile from the topic folder; if an older compiler rejects --ignoreConfig, omit that flag:

Shell
tsc --ignoreConfig --target ES2022 --module NodeNext --moduleResolution NodeNext --strict --skipLibCheck --outDir .topic-build examples/run.ts
node .topic-build/examples/run.js

The shared input fixture is byte-identical to the test fixture. The guided playground runs a bundled copy of the TypeScript core; preset names identify synthetic examples and deliberate failures.

Check your understanding

  1. Double every covariance entry. Does GMV change?
  2. Supply a different expected-return vector while keeping covariance fixed.
  3. Compare diag(0,1) with an all-ones covariance.

Answers. Weights stay fixed and variance doubles. Expected returns change only the reported return, not the objective. The first singular matrix has the unique optimum [1,0]; the second gives every feasible portfolio the same variance. Singular does not imply non-unique.

The calculation path

Stop when validation fails. A successful numerical certificate verifies the declared model, not the quality of the market estimates supplied to it.

Limitations and next variants

GMV is sensitive to the covariance estimator, window, missing-data policy, regime changes and outliers. It does not forecast returns and it does not validate semantic labels inside numeric arrays. The next family topics cover a target-return Markowitz variant, maximum Sharpe, mean-CVaR and MAD; each needs its own objective, feasible set, data contract and independent tests.

Evidence reviewed

Technical model review included the official MOSEK Portfolio Optimization Cookbook, the CFR GMV working paper, the SEC stock-split explanation, and Investor.gov's ex-dividend guidance. The covariance and stale-observation numbers in this article are synthetic fixtures with their arithmetic shown explicitly.

Decision support: diagnose before changing the solver

SymptomFirst questionCorrect next action
One asset receives nearly all weightIs its estimated variance genuinely low or created by stale/missing rows?Inspect timestamps and pairwise sample support.
Closed-form unconstrained weights include negativesIs shorting allowed in the declared mandate?If no, solve the simplex problem directly.
Covariance is singularIs the simplex objective flat, uniquely minimized, or unsupported numerically?Classify the solution set; do not assume singular means non-unique.
Daily variances are tinyWas PSD checked relative to matrix scale?Use scale-aware tolerances.
Portfolio spans currenciesWere returns translated consistently before covariance estimation?Rebuild the panel in one declared base-currency convention.
A split or dividend date dominates covarianceAre price and quantity/return adjustment bases reconciled?Repair the input history and re-estimate.

Small glossary

TermMeaning here
SimplexNonnegative weights that sum to one.
GMVThe feasible portfolio with the smallest variance.
Scale-aware toleranceA numerical threshold measured relative to the covariance magnitude.
Singular covarianceA matrix with at least one zero eigenvalue; uniqueness still needs analysis.
FW gapA first-order certificate for optimality over the simplex.
Clip-and-renormalizeAn unsafe shortcut that generally does not solve the constrained quadratic program.

Continue through Family 01

Production handoff checklist

  • Archive asset order, covariance estimator, return window, missing-row policy, horizon, and base currency.
  • Compare direct long-only optimization with the documented unconstrained and clip-renormalize counterexample.
  • Preserve matrix scale, eigen diagnostics, budget/bound residuals, FW gap, and solution-class output.
  • Investigate extreme concentration as a data symptom before adding an arbitrary weight cap.
  • Re-estimate after any change to identity, splits, dividends, FX conversion, or stale-price handling.
References7 primary sources and evidence notes

Expand the source trail, evidence role, and limitations behind the engineering choices.

The source-access observations below come from the 2026-09-14 research pass. The claim ledger distinguishes theory, calculations and implementation choices. The dated independent recheck at the end records the narrower set of sources inspected again. Search snippets do not establish the numerical contract.

R01 — MOSEK Portfolio Optimization Cookbook — Markowitz portfolio

  • Organization: MOSEK ApS
  • Source type: official solver/modeling documentation
  • URL: https://docs.mosek.com/portfolio-cookbook/markowitz.html
  • Accessed: 2026-09-14; relevant sections 2.1–2.3 and 2.4.1 read
  • Supports: quadratic portfolio objective/constraints, covariance-based model
  • Limitations: documentation is not a historical market-data source and does not establish this reference solver's tolerances or adapter policy.

R02 — On the Estimation of the Global Minimum Variance Portfolio

  • Authors/organization: Kempf & Memmel; Centre for Financial Research Cologne
  • Source type: primary working-paper PDF
  • URL: https://www.cfr-cologne.de/download/workingpaper/cfr-05-02.pdf
  • Accessed: 2026-09-14; full 35-page PDF opened, P1 and P5–P12 read
  • Supports: GMV definition and estimation-risk context
  • Limitations: estimation discussion is not a specification for a production input adapter or numerical status schema.

R03 — Portfolio Optimization / Markowitz model

R04 — SEC: Stock Splits

  • Organization: U.S. Securities and Exchange Commission
  • Source type: primary investor education
  • URL: https://www.sec.gov/answers/stocksplit.htm
  • Accessed: 2026-09-14; full 17-line page read
  • Supports: split price/share-count continuity guard used by the upstream input-integrity contract
  • Limitations: does not validate any private/provider ticker observation.

R05 — Investor.gov: Ex-Dividend Dates

R06 — UW Markowitz project notes

  • Authors/organization: University of Washington course material
  • Source type: academic lecture/project PDF
  • URL: https://sites.math.washington.edu/~burke/crs/408/fin-proj/mark1.pdf
  • Accessed: 2026-09-14; pages 1–6 read
  • Supports: classical Markowitz formulation context
  • Limitations: pedagogical notes; not used as the sole authority for any production numerical policy.

R07 — Historical/case boundary

  • The author-mentioned 1010.SR quantity/price lead is not a source. Date, acquisition/valuation basis, split history, quantity history and FX basis are missing, so it remains unverified and is excluded from examples.
  • Synthetic acceptance values are independently calculated and labeled in the staging README, article, lab and fixture JSON.

Independent recheck — 2026-09-16

MOSEK 1.6.0 §2: inspected covariance objective and constraints.

Recheck scope: theory and attribution only. Fixtures and solver limits are author-derived calculations and implementation choices. No new historical-market or personal-experience evidence is asserted.

Implementation

Full dependency-light reference implementations in both supported languages.

gmv.ts
/**
 * Pure general-N long-only global minimum-variance optimizer.
 *
 * This research-stage implementation accepts supplied moments only. It does
 * not infer returns from a holdings snapshot, fetch market data, repair
 * corporate actions, or convert currencies.
 */
export const DEFAULT_MAX_ITERATIONS = 10000;
export const TOL_SYMMETRY = 1e-12;
export const TOL_PSD = 1e-10;
export const TOL_FEASIBILITY = 1e-12;
export const TOL_VARIANCE = 1e-12;
export const TOL_CURVATURE = 1e-24;
export const TOL_LINEAR_SOLVE = 1e-14;
export const TOL_ACTIVE = 1e-12;
export const TOL_EIGEN_RESIDUAL = 1e-12;
export const TOL_STAGNATION_DIRECTION = 1e-8;
export const TOL_GAP_RELATIVE = 1e-10;
export type GMVStatus = "optimal" | "invalid_input" | "numerical_issue";
export type GMVSolutionClass = "unique" | "non_unique" | "unknown";
export interface GMVOptions {
    mu?: unknown;
    returnHorizon?: string;
    annualizationFactor?: unknown;
    maxIterations?: unknown;
}
export interface GMVResult {
    assetIds: string[] | null;
    mu: number[] | null;
    weights: number[] | null;
    variance: number | null;
    volatility: number | null;
    expectedReturn: number | null;
    budgetResidual: number | null;
    lowerBoundResidual: number | null;
    fwGapS: number | null;
    fwGap: number | null;
    matrixScale: number | null;
    status: GMVStatus;
    solutionClass: GMVSolutionClass;
    variant: "gmv-long-only-fully-invested";
    method: "projected-gradient-simplex-line-search";
    iterations: number;
    maxIterations: number;
    diagnostics: Record<string, unknown>;
    warnings: string[];
}
function isFiniteReal(value: unknown): value is number {
    return typeof value === "number" && Number.isFinite(value);
}
function dot(left: readonly number[], right: readonly number[]): number {
    let total = 0;
    for (let index = 0; index < left.length; index += 1) {
        total += left[index] * right[index];
    }
    return total;
}
function matVec(matrix: readonly (readonly number[])[], vector: readonly number[]): number[] {
    return matrix.map((row) => dot(row, vector));
}
function quadratic(matrix: readonly (readonly number[])[], vector: readonly number[]): number {
    return dot(vector, matVec(matrix, vector));
}
function maxAbsMatrix(matrix: readonly (readonly number[])[]): number {
    let maximum = 0;
    for (const row of matrix) {
        for (const value of row) {
            maximum = Math.max(maximum, Math.abs(value));
        }
    }
    return maximum;
}
function rowSumNorm(matrix: readonly (readonly number[])[]): number {
    let maximum = 0;
    for (const row of matrix) {
        let sum = 0;
        for (const value of row) {
            sum += Math.abs(value);
        }
        maximum = Math.max(maximum, sum);
    }
    return maximum;
}
function simplexProjection(values: readonly number[]): number[] {
    const ordered = values
        .map((value, index) => ({ value, index }))
        .sort((left, right) => right.value - left.value || left.index - right.index);
    let cumulative = 0;
    let rho = 0;
    let theta = 0;
    for (let position = 1; position <= ordered.length; position += 1) {
        cumulative += ordered[position - 1].value;
        const candidate = ordered[position - 1].value - (cumulative - 1) / position;
        if (candidate > 0) {
            rho = position;
            theta = (cumulative - 1) / position;
        }
    }
    if (rho === 0) {
        throw new Error("simplex projection found no positive threshold set");
    }
    return values.map((value) => Math.max(value - theta, 0));
}
function symmetricEigenvalues(matrix: readonly (readonly number[])[]): {
    values: number[];
    residual: number;
} {
    const size = matrix.length;
    if (size === 1) {
        return { values: [matrix[0][0]], residual: 0 };
    }
    const work = matrix.map((row) => [...row]);
    const limit = Math.max(32, 100 * size * size);
    for (let sweep = 0; sweep < limit; sweep += 1) {
        let p = 0;
        let q = 1;
        let largest = Math.abs(work[p][q]);
        for (let row = 0; row < size; row += 1) {
            for (let column = row + 1; column < size; column += 1) {
                const candidate = Math.abs(work[row][column]);
                if (candidate > largest) {
                    largest = candidate;
                    p = row;
                    q = column;
                }
            }
        }
        if (largest <= 1e-15) {
            break;
        }
        const app = work[p][p];
        const aqq = work[q][q];
        const apq = work[p][q];
        const angle = 0.5 * Math.atan2(2 * apq, aqq - app);
        const cosine = Math.cos(angle);
        const sine = Math.sin(angle);
        for (let index = 0; index < size; index += 1) {
            if (index === p || index === q) {
                continue;
            }
            const aip = work[index][p];
            const aiq = work[index][q];
            const newIp = cosine * aip - sine * aiq;
            const newIq = sine * aip + cosine * aiq;
            work[index][p] = newIp;
            work[p][index] = newIp;
            work[index][q] = newIq;
            work[q][index] = newIq;
        }
        work[p][p] = cosine * cosine * app - 2 * sine * cosine * apq + sine * sine * aqq;
        work[q][q] = sine * sine * app + 2 * sine * cosine * apq + cosine * cosine * aqq;
        work[p][q] = 0;
        work[q][p] = 0;
    }
    let residual = 0;
    for (let row = 0; row < size; row += 1) {
        for (let column = row + 1; column < size; column += 1) {
            residual = Math.max(residual, Math.abs(work[row][column]));
        }
    }
    return { values: work.map((row, index) => row[index]), residual };
}
function solveLinearSystem(matrix: readonly (readonly number[])[], rightHandSide: readonly number[]): number[] | null {
    const size = matrix.length;
    if (size === 0 || rightHandSide.length !== size) {
        return null;
    }
    const augmented = matrix.map((row, rowIndex) => [...row, rightHandSide[rowIndex]]);
    for (let column = 0; column < size; column += 1) {
        let pivotRow = column;
        for (let row = column + 1; row < size; row += 1) {
            if (Math.abs(augmented[row][column]) > Math.abs(augmented[pivotRow][column])) {
                pivotRow = row;
            }
        }
        const pivot = Math.abs(augmented[pivotRow][column]);
        if (!Number.isFinite(pivot) || pivot <= TOL_LINEAR_SOLVE) {
            return null;
        }
        if (pivotRow !== column) {
            [augmented[column], augmented[pivotRow]] = [augmented[pivotRow], augmented[column]];
        }
        const pivotValue = augmented[column][column];
        for (let row = column + 1; row < size; row += 1) {
            const factor = augmented[row][column] / pivotValue;
            if (factor === 0) {
                continue;
            }
            for (let entry = column; entry <= size; entry += 1) {
                augmented[row][entry] -= factor * augmented[column][entry];
            }
        }
    }
    const solution = Array(size).fill(0) as number[];
    for (let row = size - 1; row >= 0; row -= 1) {
        const diagonal = augmented[row][row];
        if (Math.abs(diagonal) <= TOL_LINEAR_SOLVE) {
            return null;
        }
        let remainder = 0;
        for (let column = row + 1; column < size; column += 1) {
            remainder += augmented[row][column] * solution[column];
        }
        solution[row] = (augmented[row][size] - remainder) / diagonal;
        if (!Number.isFinite(solution[row])) {
            return null;
        }
    }
    return solution;
}
function activeSetPolish(matrix: readonly (readonly number[])[], startingWeights: readonly number[], tolGap: number): {
    weights: number[];
    iterations: number;
} | null {
    const size = matrix.length;
    // A singular face can still have an exact vertex optimum (for example
    // diag(0, 1)).  Check vertices using their FW/KKT gap before solving a
    // nonsingular face; this is not a smallest-diagonal heuristic.
    for (let vertex = 0; vertex < size; vertex += 1) {
        const candidate = Array(size).fill(0) as number[];
        candidate[vertex] = 1;
        const gradient = matVec(matrix, candidate).map((value) => 2 * value);
        const vertexGap = gradient[vertex] - Math.min(...gradient);
        if (vertexGap <= tolGap) {
            return { weights: candidate, iterations: 0 };
        }
    }
    let active = startingWeights
        .map((value, index) => (value > TOL_ACTIVE ? index : -1))
        .filter((index) => index >= 0);
    if (active.length === 0) {
        let largest = 0;
        for (let index = 1; index < size; index += 1) {
            if (startingWeights[index] > startingWeights[largest]) {
                largest = index;
            }
        }
        active = [largest];
    }
    for (let polishIteration = 0; polishIteration <= size; polishIteration += 1) {
        const face = active.map((row) => active.map((column) => matrix[row][column]));
        const solved = solveLinearSystem(face, Array(active.length).fill(1));
        if (solved === null) {
            return null;
        }
        const denominator = solved.reduce((sum, value) => sum + value, 0);
        if (!Number.isFinite(denominator) || denominator <= TOL_LINEAR_SOLVE) {
            return null;
        }
        const faceWeights = solved.map((value) => value / denominator);
        const materiallyNegative = faceWeights
            .map((value, position) => (value < -TOL_FEASIBILITY ? position : -1))
            .filter((position) => position >= 0);
        if (materiallyNegative.length > 0) {
            if (active.length === 1) {
                return null;
            }
            const removePosition = materiallyNegative.reduce((best, position) => faceWeights[position] < faceWeights[best] ? position : best, materiallyNegative[0]);
            active.splice(removePosition, 1);
            continue;
        }
        const candidate = Array(size).fill(0) as number[];
        active.forEach((index, position) => {
            candidate[index] = Math.max(0, faceWeights[position]);
        });
        const total = candidate.reduce((sum, value) => sum + value, 0);
        if (!Number.isFinite(total) || total <= 0) {
            return null;
        }
        for (let index = 0; index < size; index += 1) {
            candidate[index] /= total;
        }
        const gradient = matVec(matrix, candidate).map((value) => 2 * value);
        const positive = candidate
            .map((value, index) => (value > TOL_ACTIVE ? index : -1))
            .filter((index) => index >= 0);
        if (positive.length === 0) {
            return null;
        }
        const multiplier = positive.reduce((sum, index) => sum + gradient[index], 0) / positive.length;
        const violating = Array.from({ length: size }, (_, index) => index)
            .filter((index) => !positive.includes(index) && gradient[index] < multiplier - tolGap);
        if (violating.length > 0) {
            const addIndex = violating.reduce((best, index) => gradient[index] < gradient[best] ? index : best, violating[0]);
            if (!active.includes(addIndex)) {
                active.push(addIndex);
                active.sort((left, right) => left - right);
                continue;
            }
        }
        return { weights: candidate, iterations: polishIteration + 1 };
    }
    return null;
}
function failure(status: GMVStatus, message: string, assetIds: string[] | null = null, mu: number[] | null = null, diagnostics: Record<string, unknown> = {}): GMVResult {
    return {
        assetIds,
        mu,
        weights: null,
        variance: null,
        volatility: null,
        expectedReturn: null,
        budgetResidual: null,
        lowerBoundResidual: null,
        fwGapS: null,
        fwGap: null,
        matrixScale: null,
        status,
        solutionClass: "unknown",
        variant: "gmv-long-only-fully-invested",
        method: "projected-gradient-simplex-line-search",
        iterations: 0,
        maxIterations: DEFAULT_MAX_ITERATIONS,
        diagnostics: { message, ...diagnostics },
        warnings: [],
    };
}
function basisColumn(index: number, basis: readonly (readonly number[])[]): number[] {
    return basis.map((row) => row[index]);
}
function solutionClass(matrix: readonly (readonly number[])[], tolerance: number): GMVSolutionClass {
    const size = matrix.length;
    if (size <= 1) {
        return "unique";
    }
    const basis = Array.from({ length: size }, () => Array(size - 1).fill(0));
    for (let index = 0; index < size - 1; index += 1) {
        basis[index][index] = 1;
        basis[size - 1][index] = -1;
    }
    const restricted = Array.from({ length: size - 1 }, (_, row) => Array.from({ length: size - 1 }, (_, column) => dot(basisColumn(row, basis), matVec(matrix, basisColumn(column, basis)))));
    const symmetric = restricted.map((row, rowIndex) => row.map((value, columnIndex) => (value + restricted[columnIndex][rowIndex]) / 2));
    // Only an exactly zero restricted matrix proves a flat objective. A
    // representable but tiny positive curvature remains numerically unresolved.
    if (maxAbsMatrix(symmetric) === 0) {
        return "non_unique";
    }
    const eigensystem = symmetricEigenvalues(symmetric);
    if (eigensystem.residual > TOL_EIGEN_RESIDUAL) {
        return "unknown";
    }
    if (Math.min(...eigensystem.values) > tolerance) {
        return "unique";
    }
    return "unknown";
}
export function globalMinimumVariance(assetIdsInput: unknown, covarianceInput: unknown, options: GMVOptions = {}): GMVResult {
    const del = options.returnHorizon;
    void del;
    if (!Array.isArray(assetIdsInput)) {
        return failure("invalid_input", "assetIds must be an ordered array.");
    }
    if (!Array.isArray(covarianceInput)) {
        return failure("invalid_input", "covariance must be a square numeric matrix.");
    }
    if (assetIdsInput.length === 0) {
        return failure("invalid_input", "assetIds must not be empty.");
    }
    const assetIds: string[] = [];
    for (const value of assetIdsInput) {
        if (typeof value !== "string" || value.trim() === "") {
            return failure("invalid_input", "assetIds must contain non-empty strings.");
        }
        assetIds.push(value);
    }
    if (new Set(assetIds).size !== assetIds.length) {
        return failure("invalid_input", "assetIds must be unique.");
    }
    const size = assetIds.length;
    if (covarianceInput.length !== size
        || covarianceInput.some((row) => !Array.isArray(row) || row.length !== size)) {
        return failure("invalid_input", "covariance must align to assetIds as an N by N matrix.", assetIds);
    }
    const matrix: number[][] = [];
    for (const row of covarianceInput) {
        if (!Array.isArray(row)) {
            return failure("invalid_input", "covariance rows must be arrays.", assetIds);
        }
        const converted: number[] = [];
        for (const value of row) {
            if (!isFiniteReal(value)) {
                return failure("invalid_input", "covariance entries must be finite real numbers.", assetIds);
            }
            converted.push(value);
        }
        matrix.push(converted);
    }
    let mu: number[] | null = null;
    if (options.mu !== undefined) {
        if (!Array.isArray(options.mu) || options.mu.length !== size) {
            return failure("invalid_input", "mu must be a finite vector aligned to assetIds.", assetIds);
        }
        if (options.mu.some((value) => !isFiniteReal(value))) {
            return failure("invalid_input", "mu entries must be finite real numbers.", assetIds);
        }
        mu = options.mu.map((value) => value as number);
    }
    if (options.annualizationFactor !== undefined
        && (!isFiniteReal(options.annualizationFactor) || options.annualizationFactor <= 0)) {
        return failure("invalid_input", "annualizationFactor must be finite and positive.", assetIds, mu);
    }
    const maxIterations = options.maxIterations ?? DEFAULT_MAX_ITERATIONS;
    if (!Number.isInteger(maxIterations) || (maxIterations as number) < 0) {
        return failure("invalid_input", "maxIterations must be a non-negative integer.", assetIds, mu);
    }
    const iterationsBudget = maxIterations as number;
    const base: Omit<GMVResult, "weights" | "variance" | "volatility" | "expectedReturn" | "budgetResidual" | "lowerBoundResidual" | "fwGapS" | "fwGap" | "diagnostics" | "warnings" | "status" | "solutionClass" | "iterations"> = {
        assetIds,
        mu,
        matrixScale: null,
        variant: "gmv-long-only-fully-invested",
        method: "projected-gradient-simplex-line-search",
        maxIterations: iterationsBudget,
    };
    const scale = maxAbsMatrix(matrix);
    base.matrixScale = scale;
    if (scale === 0) {
        const weights = Array(size).fill(1 / size);
        return {
            ...base,
            weights,
            variance: 0,
            volatility: 0,
            expectedReturn: mu === null ? null : dot(mu, weights),
            budgetResidual: Math.abs(weights.reduce((sum, value) => sum + value, 0) - 1),
            lowerBoundResidual: 0,
            fwGapS: 0,
            fwGap: 0,
            status: "optimal",
            solutionClass: size === 1 ? "unique" : "non_unique",
            iterations: 0,
            warnings: ["all_zero_covariance_flat_objective"],
            diagnostics: {
                matrixScale: 0,
                normalizedObjective: 0,
                optimalityCertificate: "exact-flat-objective",
            },
        };
    }
    let normalized = matrix.map((row) => row.map((value) => value / scale));
    let symmetryResidual = 0;
    for (let row = 0; row < size; row += 1) {
        for (let column = 0; column < size; column += 1) {
            symmetryResidual = Math.max(symmetryResidual, Math.abs(normalized[row][column] - normalized[column][row]));
        }
    }
    if (symmetryResidual > TOL_SYMMETRY) {
        return failure("invalid_input", "covariance symmetry residual exceeds the normalized tolerance.", assetIds, mu, { matrixScale: scale, symmetryResidualS: symmetryResidual });
    }
    normalized = normalized.map((row, rowIndex) => row.map((value, columnIndex) => (value + normalized[columnIndex][rowIndex]) / 2));
    const eigensystem = symmetricEigenvalues(normalized);
    const minimumEigenvalue = Math.min(...eigensystem.values);
    const maximumEigenvalue = Math.max(...eigensystem.values);
    const diagnosticsBase: Record<string, unknown> = {
        matrixScale: scale,
        symmetryResidualS: symmetryResidual,
        minEigenvalueS: minimumEigenvalue,
        maxEigenvalueS: maximumEigenvalue,
        eigenResidualS: eigensystem.residual,
        tolSymmetry: TOL_SYMMETRY,
        tolPsd: TOL_PSD,
        tolEigenResidual: TOL_EIGEN_RESIDUAL,
        tolFeasibility: TOL_FEASIBILITY,
        tolVariance: TOL_VARIANCE,
        tolCurvature: TOL_CURVATURE,
        tolActive: TOL_ACTIVE,
        tolLinearSolve: TOL_LINEAR_SOLVE,
        tolStagnationDirection: TOL_STAGNATION_DIRECTION,
        tolGapRelative: TOL_GAP_RELATIVE,
    };
    if (eigensystem.residual > TOL_EIGEN_RESIDUAL) {
        return failure("numerical_issue", "symmetric eigensolver residual exceeds the fixed normalized tolerance.", assetIds, mu, diagnosticsBase);
    }
    if (minimumEigenvalue < -TOL_PSD) {
        return failure("invalid_input", "normalized covariance is materially indefinite.", assetIds, mu, diagnosticsBase);
    }
    const warnings: string[] = [];
    if (symmetryResidual > 0) {
        warnings.push("symmetrized_within_tolerance");
    }
    if (minimumEigenvalue < 0) {
        warnings.push("psd_within_roundoff");
    }
    const rowNorm = rowSumNorm(normalized);
    const tolGap = TOL_GAP_RELATIVE * Math.max(1, rowNorm);
    let weights = Array(size).fill(1 / size);
    let normalizedObjective = quadratic(normalized, weights);
    let gap = Number.POSITIVE_INFINITY;
    let budgetResidual = Math.abs(weights.reduce((sum, value) => sum + value, 0) - 1);
    let lowerBoundResidual = 0;
    const className = solutionClass(normalized, TOL_PSD);
    let iterations = 0;
    let status: GMVStatus = "numerical_issue";
    let polishedActiveSet = false;
    let activeSetPolishIterations = 0;
    if (iterationsBudget === 0) {
        const gradient = matVec(normalized, weights).map((value) => 2 * value);
        gap = dot(gradient, weights) - Math.min(...gradient);
    }
    else if (maximumEigenvalue <= TOL_VARIANCE) {
        warnings.push("flat_objective_from_eigenvalue_tolerance");
        gap = 0;
        const variance = Math.max(0, scale * normalizedObjective);
        return {
            ...base,
            weights,
            variance,
            volatility: Math.sqrt(variance),
            expectedReturn: mu === null ? null : dot(mu, weights),
            budgetResidual,
            lowerBoundResidual,
            fwGapS: 0,
            fwGap: 0,
            status: "optimal",
            solutionClass: className,
            iterations: 0,
            warnings,
            diagnostics: {
                ...diagnosticsBase,
                normalizedObjective,
                tolGap,
                optimalityCertificate: "simplex-frank-wolfe-convex-certificate",
                polishedActiveSet: false,
                activeSetPolishIterations: 0,
                projection: "euclidean-simplex-sort-threshold",
            },
        };
    }
    else {
        const step = 1 / (2 * maximumEigenvalue);
        for (let iteration = 0; iteration < iterationsBudget; iteration += 1) {
            const gradient = matVec(normalized, weights).map((value) => 2 * value);
            const trial = weights.map((value, index) => value - step * gradient[index]);
            let projected: number[];
            try {
                projected = simplexProjection(trial);
            }
            catch {
                break;
            }
            const projectionBudget = Math.abs(projected.reduce((sum, value) => sum + value, 0) - 1);
            const projectionLower = Math.max(0, -Math.min(...projected));
            if (projectionBudget > TOL_FEASIBILITY || projectionLower > TOL_FEASIBILITY) {
                break;
            }
            const direction = projected.map((value, index) => value - weights[index]);
            const linear = dot(gradient, direction);
            const curvature = quadratic(normalized, direction);
            let alpha: number;
            if (curvature > TOL_CURVATURE) {
                alpha = Math.min(1, Math.max(0, -linear / (2 * curvature)));
            }
            else if (linear < 0) {
                alpha = 1;
            }
            else {
                alpha = 0;
            }
            const stagnated = alpha === 0
                && Math.max(...direction.map((value) => Math.abs(value))) <= TOL_STAGNATION_DIRECTION;
            weights = weights.map((value, index) => value + alpha * direction[index]);
            const currentIterations = iteration + 1;
            iterations = currentIterations;
            normalizedObjective = quadratic(normalized, weights);
            const currentGradient = matVec(normalized, weights).map((value) => 2 * value);
            gap = dot(currentGradient, weights) - Math.min(...currentGradient);
            budgetResidual = Math.abs(weights.reduce((sum, value) => sum + value, 0) - 1);
            lowerBoundResidual = Math.max(0, -Math.min(...weights));
            if (gap <= tolGap
                && budgetResidual <= TOL_FEASIBILITY
                && lowerBoundResidual <= TOL_FEASIBILITY
                && normalizedObjective >= -TOL_VARIANCE) {
                status = "optimal";
                break;
            }
            if (stagnated) {
                // Stop repeating floating-point cancellation; the KKT polish below
                // must still earn the FW certificate before status is optimal.
                break;
            }
        }
    }
    if (iterationsBudget > 0 && Number.isFinite(normalizedObjective)) {
        const polished = activeSetPolish(normalized, weights, tolGap);
        if (polished !== null) {
            const candidate = polished.weights;
            const candidateObjective = quadratic(normalized, candidate);
            const candidateGradient = matVec(normalized, candidate).map((value) => 2 * value);
            const candidateGap = dot(candidateGradient, candidate) - Math.min(...candidateGradient);
            const candidateBudget = Math.abs(candidate.reduce((sum, value) => sum + value, 0) - 1);
            const candidateLower = Math.max(0, -Math.min(...candidate));
            if (candidateObjective <= normalizedObjective + TOL_VARIANCE
                && candidateGap <= tolGap
                && candidateBudget <= TOL_FEASIBILITY
                && candidateLower <= TOL_FEASIBILITY
                && candidateObjective >= -TOL_VARIANCE) {
                weights = candidate;
                normalizedObjective = candidateObjective;
                gap = candidateGap;
                budgetResidual = candidateBudget;
                lowerBoundResidual = candidateLower;
                status = "optimal";
                polishedActiveSet = true;
                activeSetPolishIterations = polished.iterations;
            }
        }
    }
    if (normalizedObjective < -TOL_VARIANCE) {
        status = "numerical_issue";
        const varianceRaw = scale * normalizedObjective;
        return {
            ...base,
            weights,
            variance: null,
            volatility: null,
            expectedReturn: mu === null ? null : dot(mu, weights),
            budgetResidual,
            lowerBoundResidual,
            fwGapS: Number.isFinite(gap) ? gap : null,
            fwGap: Number.isFinite(gap) ? scale * gap : null,
            status: "numerical_issue",
            solutionClass: className,
            iterations,
            warnings,
            diagnostics: {
                ...diagnosticsBase,
                normalizedObjective,
                tolGap,
                varianceRaw,
                polishedActiveSet,
                activeSetPolishIterations,
                optimalityCertificate: "not-certified",
            },
        };
    }
    const varianceRaw = scale * normalizedObjective;
    const variance = normalizedObjective < 0 ? 0 : varianceRaw;
    const finalWarnings = normalizedObjective < 0 ? [...warnings, "roundoff_variance_zeroed"] : warnings;
    return {
        ...base,
        weights,
        variance,
        volatility: Math.sqrt(Math.max(variance, 0)),
        expectedReturn: mu === null ? null : dot(mu, weights),
        budgetResidual,
        lowerBoundResidual,
        fwGapS: Number.isFinite(gap) ? gap : null,
        fwGap: Number.isFinite(gap) ? scale * gap : null,
        status,
        solutionClass: className,
        iterations,
        warnings: finalWarnings,
        diagnostics: {
            ...diagnosticsBase,
            normalizedObjective,
            tolGap,
            fwGapSRaw: gap,
            varianceRaw,
            projection: "euclidean-simplex-sort-threshold",
            polishedActiveSet,
            activeSetPolishIterations,
            optimalityCertificate: iterationsBudget === 0
                ? "not-certified-max-iterations-zero"
                : minimumEigenvalue < 0 && status === "optimal"
                    ? "tolerance-limited-near-psd"
                    : status === "optimal"
                        ? "simplex-frank-wolfe-convex-certificate"
                        : "not-certified",
        },
    };
}
Playground

The embedded lab now expands to its full document height, keeping the article as the only scroll surface.

Full-height labplaygroundOpen full screen
Written by

Fintech engineer building market-data and financial systems, and the author of every article, glossary record, and reference implementation on The Fintech Builder.

Continue the system

Related released builds.

Back to the released library
01

Portfolio Construction / Mean-Risk Optimization

Markowitz Mean-Variance

A production-minded guide to Markowitz Mean-Variance.
10 min Python + TypeScriptAdvanced
02

Portfolio Construction / Mean-Risk Optimization

Maximum Sharpe Ratio

A production-minded guide to Maximum Sharpe Ratio.
8 min Python + TypeScriptAdvanced
03

Portfolio Construction / Mean-Risk Optimization

Mean-CVaR Optimization

A production-minded guide to Mean-CVaR Optimization.
9 min Python + TypeScriptAdvanced