Markowitz Mean-Variance: When a Return Floor Adds Risk

How much extra risk does a return target force you to take? Start with the lowest-variance allocation. If it already clears the target, keep it. If it does not, move only far enough toward higher expected return to meet the floor. This is the decision the tutorial computes and audits.

This tutorial uses a small synthetic fixture. It does not turn an unverified holding or a vendor snapshot into a historical case.

Before you start

This tutorial is for developers comfortable with weighted averages, vectors, and basic covariance. By the end you can derive the return-floor boundary, distinguish slack from binding targets, and verify the minimum-variance result. 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 “Markowitz” means here

Markowitz's portfolio model is commonly taught through weighted expected return and covariance risk. The University of Washington teaching note writes the portfolio mean and variance in this form and introduces a target-return quadratic program (the note, pages 0–3). MOSEK's portfolio cookbook shows that target-return minimum variance, risk-bound maximum return, and mean-minus-variance utility are different formulations (MOSEK §2).

This article chooses one of them:

\begin{aligned} \min_w\quad &w^{\mathsf T}\Sigma w\\ \text{subject to}\quad &\mathbf{1}^{\mathsf T}w=1,\\ &\mu^{\mathsf T}w\ge\tau,\qquad w_i\ge0. \end{aligned}

w is the vector of portfolio weights, μ is the expected return vector for one declared horizon, Σ is the covariance matrix for that same return unit, and τ is the return floor. The non-negativity and budget constraints define a fully invested, long-only simplex. There is no shorting, leverage, cash position, fee model, turnover penalty, or risk-free asset in this variant.

The target is deliberately an inequality. If the global minimum-variance (GMV) portfolio already clears the floor, the target is slack. Replacing ≥ with = would silently change that decision problem.

The input is part of the algorithm

The core accepts supplied assetIds, μ, Σ, and τ. It does not infer a history from holdings. Before a caller estimates moments, it should be able to answer these questions:

Is each number a current valuation mark or an acquisition cost? A purchase price can be a valid accounting fact while being unsuitable as a current return observation.
If a split occurred, were quantity and per-share price transformed together? A split-adjusted price paired with raw quantity, or the reverse, changes economic exposure.
What are the local and base currencies, FX quote direction, rate date, local units, and portfolio cutoff? “Convert to USD” is not a complete contract.
Are dividends represented as cash/receivables, a price return, or a total-return/adjusted series? Adding a dividend to a series that already includes it counts the same wealth twice. Investor.gov's ex-dividend timing explainer is a useful primary reference for separating entitlement and payment dates (Investor.gov).
Do every return observation and covariance entry share the same horizon, calendar, revision rule, and point-in-time cutoff? A backward-filled future quote is not evidence available at an earlier decision time.

These are upstream lineage checks. Numeric vectors alone cannot reveal whether someone mislabeled returns as losses, inverted an FX quote, or used cost instead of mark. The solver reports what it was given; it does not certify economics it never received.

A reproducible two-asset fixture

The following inputs are synthetic and use the same declared return horizon:

Plain text

assets       = [A, B]
mu           = [0.06, 0.14]
Sigma        = [[0.04, 0.01],
                [0.01, 0.09]]
target tau   = 0.10

With two assets, the budget gives w_B = 1 − w_A. The target floor becomes:

Plain text

0.06 w_A + 0.14(1 − w_A) ≥ 0.10
w_A ≤ 0.50

The same covariance has a long-only GMV point at w = [8/11, 3/11]. Its expected return is 0.0818181818…, below the target, so the target is binding at [0.50, 0.50].

Plain text

expected return = 0.50·0.06 + 0.50·0.14 = 0.10
variance        = 0.25·0.04 + 2·0.25·0.01 + 0.25·0.09 = 0.0375
volatility      = sqrt(0.0375) = 0.1936491673…

The budget residual and target residual are zero in exact arithmetic. The implementation also reports a normalized Frank–Wolfe gap: with g = 2Σ̃w, it compares gᵀw with the minimum of the same linearization over the complete feasible set F (including the return floor). For a PSD covariance, a non-negative gap below the fixed tolerance bounds the global objective error for this convex program. If a tiny negative eigenvalue is accepted as roundoff, the certificate is tolerance-limited rather than an exact convex proof.

Synthetic Markowitz risk-return opportunity curve showing the GMV comparator and target-boundary portfolio

The full-size opportunity-curve SVG is a synthetic two-asset diagram. It shows the distinction between the lower-risk GMV comparator (which clears floors only up to 0.0818…) and the target-floor point (floor binding); the full curve is not an efficient frontier for every possible floor, and it is not a forecast or a claim about an investable market frontier.

Boundaries worth testing

Change only τ while keeping the moments fixed:

Target	Result	Why
`0.08`	GMV `[8/11,3/11]`; expected return `0.081818…`; floor slack	The unconstrained-within-simplex risk minimum already clears the floor.
`0.10`	`[0.5,0.5]`; variance `0.0375`; floor binds	The GMV point falls below the requested return.
`0.15`	`infeasible`	No long-only fully invested portfolio can exceed the maximum asset return `0.14`.

An infeasible target is not a numerical failure and should not be replaced by equal weights. Conversely, a numerical_issue means the solver did not earn its certificate; it is not permission to clip a negative variance or silently regularize a matrix.

What the covariance gate checks

The reference implementation scales the matrix by its actual maximum absolute entry and checks symmetry, eigenvalue, and eigensolver residuals in the dimensionless matrix. It accepts a symmetry residual and eigensolver residual up to 1e−12, and a minimum eigenvalue down to −1e−10 for roundoff. A materially indefinite covariance returns invalid_input. No eigenvalue clipping or ridge is applied.

Return values are scaled separately by the actual maximum absolute value among μ and τ; this preserves the constraint when returns are tiny. Weight feasibility, return-row tolerance, variance tolerance, and FW gap tolerance are separate controls. That distinction matters: using max(1, …) for a tiny return row can make the target disappear numerically.

For PSD singular matrices, the result may not be unique. The package returns a deterministic representative and labels non_unique only for an exact flat direction or a boundary with multiple feasible minimizers. Near-flat but strictly positive curvature is unknown, not a made-up uniqueness claim.

A prior-method lesson: stale observations can manufacture “riskless” weights

The companion stale-quote case is a controlled synthetic calculation, not a market anecdote. Four aligned returns were:

Plain text

A = [-0.02, -0.01,  0.01,  0.02]
B = [-0.01,  0.01, -0.01,  0.01]
Sigma = [[1/3000, 1/15000],
         [1/15000, 1/7500]]

The clean covariance gives GMV weights [0.2, 0.8] and variance 3/25000. Replacing the A returns with imputed zeros produces a diagonal covariance with zero A variance, so the displayed GMV becomes [1,0] with variance zero. That is a consequence of the altered input, not proof of a riskless investment. The correct diagnostic is to expose the missing-data evidence and timestamp before estimating the covariance. The prior D01-F02-A04 Stale-Quote Detector provides the data-quality method. This D14 fixture is a new, controlled synthetic consequence calculation, not a claim that a particular public holding or historical quote behaved this way.

The same accounting principle applies to corporate actions and cash. A split must reconcile price and quantity. A cash dividend must reconcile stock value, receivable, and cash. An adjusted total-return series must not receive the same dividend again. These checks improve reproducibility and diagnosis; they do not make the optimizer a data vendor.

Implementation and independent checks

The Python and TypeScript implementations are kept as separate, reviewable ports with a shared JSON fixture:

From the topic directory, this is a copy/paste Python entry example using the same canonical inputs:

Python

from implementations.python.markowitz import markowitz_mean_variance

result = markowitz_mean_variance(
    ["A", "B"],
    [0.06, 0.14],
    [[0.04, 0.01], [0.01, 0.09]],
    0.10,
)
print(result["status"], result["weights"], result["variance"])
# optimal [0.5, 0.5] 0.0375

The solver uses an active-set KKT path. If a singular face blocks that path, it uses exact support enumeration for small N, and an exact-line-search Frank–Wolfe fallback for larger N. The output exposes status, weights, expected return, variance, volatility, budget/target/bound residuals, matrix scale, eigen diagnostics, iterations, and both normalized/original FW gaps.

The tests include general multi-asset cases, scale stress, malformed source types, an indefinite matrix, singular unique/non-unique cases, a near-flat case, a true singleton return boundary, and a four-asset regression where clipping a shorting solution and renormalizing is not the direct long-only answer. An independent reviewer compares objective and feasibility against support enumeration; the two implementations are also compiled and run against the same fixtures.

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

Lower the floor from 10% to 8%. Does the solution need to earn exactly 8%?
Raise it to 14%, then 15%. Which boundary is feasible?
Multiply every return and the target by 100 and covariance by 10,000. Should weights change?

Answers. The 8% floor is slack: GMV earns about 8.1818%. At 14%, only B is feasible; 15% is infeasible. Consistent unit conversion preserves weights; variance scales by 10,000. Converting only one input changes the model.

The calculation path

Rendering system map…

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

Summary

Markowitz mean-variance is a precise convex program only after its feasible set, horizon, units, covariance basis, and target are declared. For the canonical fixture, [0.5,0.5] is the target-active minimum-variance answer, while [8/11,3/11] is the slack-target GMV comparator. The useful engineering habit is to show the inputs and certificate that support the answer—and to stop at the lineage gate when marks, splits, FX, dividends, or stale observations do not reconcile.

Decision support: what should change the answer?

Observation	Correct response	Tempting but unsafe shortcut
Target is below the GMV return	Report that the constraint is slack and return GMV	Force a point onto the active frontier anyway
Target exceeds the largest long-only asset mean	Return `infeasible`	Extrapolate weights or introduce hidden leverage
Covariance is materially indefinite	Reject the model input	Clip eigenvalues without declaring a new estimator
Unconstrained solution has a negative weight	Solve the long-only problem directly	Clip the weight and renormalize
Price and quantity use different split bases	Stop upstream and reconcile history	Let covariance “average out” the error
Assets use different currencies or horizons	Translate and align before estimating moments	Mix local returns and a USD target

Small glossary

Term	Meaning in this article
Efficient frontier	Portfolios with minimum variance for each feasible expected return.
Active target	A return floor met at equality because relaxing it would reduce variance.
Slack target	A return floor already exceeded by the GMV solution.
PSD covariance	A covariance matrix whose quadratic variance is never negative.
KKT certificate	Feasibility and stationarity evidence for the declared constrained optimum.
FW gap	A scale-aware first-order gap over the simplex with the return floor included.

Continue through Family 01

Global Minimum Variance removes the expected-return target and isolates covariance risk.
Maximum Sharpe Ratio replaces the return floor with an excess-return-to-volatility ratio.
Mean-CVaR Optimization uses an explicit loss tail rather than squared dispersion.
Mean-Absolute-Deviation Optimization uses probability-weighted absolute scenario deviations.

References

Technical and input-lineage sources are listed with access extent in REFERENCES.md; the package’s research record preserves the discovery and destination-access observations separately from this reader-facing tutorial.

Production handoff checklist

Freeze the ordered asset universe, moment horizon, base currency, target, and feasible set.
Reconcile split-adjusted prices with quantities and document dividend treatment before estimating returns.
Record covariance scale, PSD diagnostics, target feasibility, primal residuals, and optimality gap with the weights.
Re-run the fixture and at least one infeasible, singular, and clip-renormalize regression before release.
Treat a later change to estimator, bounds, cash, shorting, leverage, or costs as a new model version.

ReferencesPrimary sources and evidence notes

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

The source-access observations below were recorded during the 2026-09-14 research pass. The local claim ledger separates source theory, synthetic calculations and implementation choices. The dated independent recheck at the end records what was inspected again; it does not claim to repeat every earlier discovery search.

Primary technical sources

R01 — University of Washington, “Markowitz Mean-Variance Portfolio Theory”

Organization: University of Washington, Department of Mathematics.
Source type: academic course note (PDF), primary teaching formulation.
URL: https://sites.math.washington.edu/~burke/crs/408/fin-proj/mark1.pdf
Accessed: 2026-09-14; full 17-page PDF opened; pages 0–3/readable extracted lines 98–151 inspected.
Supports: weighted portfolio return, expected return μᵀw, covariance risk wᵀΣw, and target-return quadratic-program setup.
Limitations: classroom note; it does not define this package's input provenance, tolerance policy, or cross-language certificate.

R02 — MOSEK, “Markowitz portfolio optimization,” Portfolio Optimization Cookbook §2

Organization: MOSEK ApS.
Source type: official modeling documentation.
URL: https://docs.mosek.com/portfolio-cookbook/markowitz.html
Version/date observed: page states last updated 2025-11-05.
Accessed: 2026-09-14; full loaded accessibility-tree page inspected, with §§2.1–2.3 and §2.4.1/equations (2.1)–(2.4) read.
Supports: target-return minimum variance, risk-bound and utility variants, full-investment and long-only constraints, PSD/Cholesky considerations, and practical infeasibility warnings.
Limitations: a modeling cookbook, not a claim that any particular market moments or optimizer weights are economically correct.

R03 — Sharma et al., “Decision-Focused Learning for Mean-Variance Portfolio Optimization”

Authors: Sharma et al.
Source type: primary research preprint (arXiv HTML).
URL: https://arxiv.org/html/2409.09684v1
Accessed: 2026-09-14; full 488-line HTML opened; lines 34–97 and 282–323 inspected for MVO framing, estimation uncertainty, and objective context.
Supports: the warning that estimated expected returns/covariances can affect portfolio decisions; it motivates visible data-window/provenance diagnostics.
Limitations: it is not the package's definition, proof of performance, or license for a historical fixture.

R04 — Harry Markowitz, “Portfolio Selection”

Source type: foundational bibliographic record.
DOI: https://doi.org/10.1111/j.1540-6261.1952.tb01525.x
Accessed: 2026-09-14; title, journal, volume, issue, and pages 77–91 verified in DOI metadata; full article text was not accessible in this pass.
Supports: attribution only. No uninspected full-text claim is used.

Input-lineage sources

R05 — U.S. SEC, “Stock Splits”

URL: https://www.sec.gov/answers/stocksplit.htm
Accessed: 2026-09-14; full 17-line investor page read.
Supports: a stock split changes shares and per-share price mechanically; the package uses it only to motivate a paired quantity/price reconciliation.
Limitations: provider-specific adjustment conventions and dates still require an auditable source record.

R06 — Investor.gov, “Ex-Dividend Dates: When Are You Entitled to Stock and Its Dividend?”

URL: https://www.investor.gov/introduction-investing/investing-basics/glossary/ex-dividend-dates-when-are-you-entitled-stock-and
Accessed: 2026-09-14; relevant page lines 151–185 read.
Supports: entitlement/payment timing belongs in a cash/receivable ledger and must be kept distinct from a price-only series.
Limitations: it does not define a universal adjusted-price vendor formula.

Contextual Google destinations inspected

These pages were opened to identify reader gaps, not to support technical claims or rank sources:

Medium — Portfolio Optimization: The Markowitz Mean-Variance Model — full rendered article, AX nodes 32–267; package-driven max_sharpe example, no target-return/PSD/certificate contract.
QuestDB — Mean-Variance Optimization — full page, AX nodes 13–167; glossary and operational context, no reproducible fixture.
GitHub — Portfolio optimization using PyPortfolioOpt — public README/repo, AX nodes 82–176; source notebook/PDF listed but not separately opened.
ResearchGate/MDPI — Predictive stock selection with MVO — publisher panel and visible article beginning, AX text 90–753; combined prediction/preselection/MV, not a pure solver contract.
Wikipedia — Modern portfolio theory — full article accessibility tree, AX nodes 80–4,000; terminology map only and tertiary.
MBrenndoerfer — Modern Portfolio Theory and Mean-Variance Optimization — full interactive chapter, AX nodes 19–2,740; explanatory variants, no strict schema/certificate.
Springer — Markowitz Mean-Variance Optimization — full chapter §§7.1–7.5 through Egyptian Knowledge Bank, AX nodes 28–588; useful context, not this package's general-N numerical contract.

Citation boundaries

No historical market-data example is included. The synthetic fixtures are independently calculated and clearly labeled. The author's unverified 1010.SR lead is intentionally absent. The core does not claim to detect semantic errors in unlabeled numeric vectors; it requires explicit provenance metadata upstream.

Independent recheck — 2026-09-16

MOSEK 1.6.0 §2: inspected model and constraint definitions.

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.

markowitz.ts

/**
 * Long-only target-return Markowitz optimizer.
 *
 * This core consumes supplied expected returns and covariance moments only.
 * It does not infer moments from holdings, repair corporate actions, convert
 * currencies, or silently regularize a covariance matrix.
 */
export const DEFAULT_MAX_ITERATIONS = 10000;
export const TOL_SYMMETRY = 1e-12;
export const TOL_PSD = 1e-10;
export const TOL_EIGEN_RESIDUAL = 1e-12;
export const TOL_FEASIBILITY = 1e-12;
export const TOL_RETURN_SCALED = 1e-12;
export const TOL_VARIANCE = 1e-12;
export const TOL_LINEAR_SOLVE = 1e-14;
export const TOL_ACTIVE = 1e-12;
export const TOL_GAP_RELATIVE = 1e-10;
export const MAX_EXACT_SUPPORT_ASSETS = 12;
export type MarkowitzStatus = "optimal" | "infeasible" | "invalid_input" | "numerical_issue";
export type MarkowitzSolutionClass = "unique" | "non_unique" | "unknown";
export interface MarkowitzOptions {
    maxIterations?: unknown;
}
export interface MarkowitzResult {
    assetIds: string[] | null;
    mu: number[] | null;
    targetReturn: number | null;
    weights: number[] | null;
    expectedReturn: number | null;
    variance: number | null;
    volatility: number | null;
    budgetResidual: number | null;
    targetResidual: number | null;
    lowerBoundResidual: number | null;
    fwGapS: number | null;
    fwGap: number | null;
    matrixScale: number | null;
    status: MarkowitzStatus;
    solutionClass: MarkowitzSolutionClass;
    variant: "markowitz-target-return-min-variance-long-only";
    method: "active-set-kkt-face-enumeration";
    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 lexicographicLess(left: readonly number[], right: readonly number[]): boolean {
    for (let index = 0; index < left.length; index += 1) {
        if (left[index] !== right[index]) {
            return left[index] < right[index];
        }
    }
    return false;
}
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) {
                if (Math.abs(work[row][column]) > largest) {
                    largest = Math.abs(work[row][column]);
                    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;
            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 solutionClass(matrix: readonly (readonly number[])[]): MarkowitzSolutionClass {
    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 column = (index: number): number[] => basis.map((row) => row[index]);
    const restricted = Array.from({ length: size - 1 }, (_, row) => Array.from({ length: size - 1 }, (_, other) => dot(column(row), matVec(matrix, column(other)))));
    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 unresolved.
    if (maxAbsMatrix(symmetric) === 0) {
        return "non_unique";
    }
    const eigensystem = symmetricEigenvalues(symmetric);
    if (eigensystem.residual <= TOL_EIGEN_RESIDUAL && Math.min(...eigensystem.values) > TOL_PSD) {
        return "unique";
    }
    return "unknown";
}
function exactDuplicatePair(matrix: readonly (readonly number[])[]): [
    number,
    number
] | null {
    const size = matrix.length;
    for (let left = 0; left < size; left += 1) {
        for (let right = left + 1; right < size; right += 1) {
            let equal = true;
            for (let index = 0; index < size; index += 1) {
                if (matrix[left][index] !== matrix[right][index] || matrix[index][left] !== matrix[index][right]) {
                    equal = false;
                    break;
                }
            }
            if (equal) {
                return [left, right];
            }
        }
    }
    return null;
}
function classifySolution(matrix: readonly (readonly number[])[], muScaled: readonly number[], targetScaled: number, weights: readonly number[]): MarkowitzSolutionClass {
    const size = matrix.length;
    if (size <= 1) {
        return "unique";
    }
    const maximum = Math.max(...muScaled);
    const maximumAssets = muScaled
        .map((value, index) => Math.abs(value - maximum) <= TOL_RETURN_SCALED ? index : -1)
        .filter((index) => index >= 0);
    const expected = dot(muScaled, weights);
    const targetBinding = Math.abs(expected - targetScaled) <= TOL_RETURN_SCALED;
    if (targetBinding && Math.abs(targetScaled - maximum) <= TOL_RETURN_SCALED && maximumAssets.length === 1) {
        return "unique";
    }
    const duplicate = exactDuplicatePair(matrix);
    if (duplicate !== null) {
        const [left, right] = duplicate;
        if (weights[left] + weights[right] > TOL_ACTIVE && (!targetBinding || Math.abs(muScaled[left] - muScaled[right]) <= TOL_RETURN_SCALED)) {
            return "non_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 column = (index: number): number[] => basis.map((row) => row[index]);
    const restricted = Array.from({ length: size - 1 }, (_, row) => Array.from({ length: size - 1 }, (_, other) => dot(column(row), matVec(matrix, column(other)))));
    const symmetric = restricted.map((row, rowIndex) => row.map((value, columnIndex) => (value + restricted[columnIndex][rowIndex]) / 2));
    if (maxAbsMatrix(symmetric) === 0) {
        return "non_unique";
    }
    const eigensystem = symmetricEigenvalues(symmetric);
    if (eigensystem.residual <= TOL_EIGEN_RESIDUAL && Math.min(...eigensystem.values) > TOL_PSD) {
        return "unique";
    }
    return "unknown";
}
function feasibleLinearVertices(mu: readonly number[], target: number, tolReturn: number): number[][] {
    const size = mu.length;
    const candidates: number[][] = [];
    for (let index = 0; index < size; index += 1) {
        if (mu[index] >= target - tolReturn) {
            const candidate = Array(size).fill(0) as number[];
            candidate[index] = 1;
            candidates.push(candidate);
        }
    }
    for (let left = 0; left < size; left += 1) {
        for (let right = left + 1; right < size; right += 1) {
            const denominator = mu[right] - mu[left];
            if (denominator === 0) {
                continue;
            }
            const weightLeft = (mu[right] - target) / denominator;
            const weightRight = (target - mu[left]) / denominator;
            if (weightLeft < -tolReturn || weightRight < -tolReturn) {
                continue;
            }
            const candidate = Array(size).fill(0) as number[];
            candidate[left] = Math.max(0, weightLeft);
            candidate[right] = Math.max(0, weightRight);
            const total = candidate[left] + candidate[right];
            if (total > 0 && Math.abs(dot(mu, candidate) - target) <= 10 * tolReturn) {
                candidates.push(candidate.map((value) => value / total));
            }
        }
    }
    return candidates;
}
function linearMinimum(gradient: readonly number[], mu: readonly number[], target: number, tolReturn: number): {
    value: number;
    point: number[];
} | null {
    const candidates = feasibleLinearVertices(mu, target, tolReturn);
    if (candidates.length === 0) {
        return null;
    }
    let best = candidates[0];
    let bestValue = dot(gradient, best);
    for (const candidate of candidates.slice(1)) {
        const value = dot(gradient, candidate);
        if (value < bestValue || (value === bestValue && lexicographicLess(candidate, best))) {
            best = candidate;
            bestValue = value;
        }
    }
    return { value: bestValue, point: best };
}
function gmvFace(matrix: readonly (readonly number[])[], maxIterations: number, tolGap: number): {
    weights: number[] | null;
    iterations: number;
} {
    const size = matrix.length;
    const active = Array.from({ length: size }, (_, index) => index);
    for (let iteration = 1; iteration <= maxIterations; iteration += 1) {
        const face = active.map((row) => active.map((column) => matrix[row][column]));
        const solved = solveLinearSystem(face, Array(active.length).fill(1));
        if (solved === null) {
            break;
        }
        const denominator = solved.reduce((sum, value) => sum + value, 0);
        if (!Number.isFinite(denominator) || denominator <= TOL_LINEAR_SOLVE) {
            break;
        }
        const faceWeights = solved.map((value) => value / denominator);
        const negatives = faceWeights
            .map((value, index) => value < -TOL_FEASIBILITY ? index : -1)
            .filter((index) => index >= 0);
        if (negatives.length > 0) {
            if (active.length === 1) {
                break;
            }
            const removePosition = negatives.reduce((best, index) => faceWeights[index] < faceWeights[best] ? index : best, negatives[0]);
            active.splice(removePosition, 1);
            continue;
        }
        const weights = Array(size).fill(0) as number[];
        active.forEach((index, position) => { weights[index] = Math.max(0, faceWeights[position]); });
        const total = weights.reduce((sum, value) => sum + value, 0);
        if (!Number.isFinite(total) || total <= 0) {
            break;
        }
        for (let index = 0; index < size; index += 1) {
            weights[index] /= total;
        }
        const gradient = matVec(matrix, weights).map((value) => 2 * value);
        const positive = weights
            .map((value, index) => value > TOL_ACTIVE ? index : -1)
            .filter((index) => index >= 0);
        if (positive.length === 0) {
            break;
        }
        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, iterations: iteration };
    }
    return { weights: null, iterations: Math.min(maxIterations, size + 1) };
}
function targetFace(matrix: readonly (readonly number[])[], muScaled: readonly number[], targetScaled: number, maxIterations: number, tolGap: number): {
    weights: number[] | null;
    iterations: number;
    eta: number | null;
} {
    const size = matrix.length;
    const feasible = feasibleLinearVertices(muScaled, targetScaled, TOL_RETURN_SCALED);
    if (feasible.length === 0) {
        return { weights: null, iterations: 0, eta: null };
    }
    const exactVertex = feasible.find((candidate) => candidate.filter((value) => value > TOL_ACTIVE).length === 1
        && Math.abs(dot(muScaled, candidate) - targetScaled) <= TOL_RETURN_SCALED);
    const initial = exactVertex ?? feasible.find((candidate) => candidate.filter((value) => value > TOL_ACTIVE).length >= 2) ?? feasible[0];
    let active = initial
        .map((value, index) => value > TOL_ACTIVE ? index : -1)
        .filter((index) => index >= 0);
    if (active.length === 0) {
        return { weights: null, iterations: 0, eta: null };
    }
    for (let iteration = 1; iteration <= maxIterations; iteration += 1) {
        if (active.length === 1) {
            const index = active[0];
            if (Math.abs(muScaled[index] - targetScaled) <= TOL_FEASIBILITY) {
                const weights = Array(size).fill(0) as number[];
                weights[index] = 1;
                return { weights, iterations: iteration, eta: 0 };
            }
            break;
        }
        const count = active.length;
        const kkt = Array.from({ length: count + 2 }, () => Array(count + 2).fill(0));
        for (let row = 0; row < count; row += 1) {
            const asset = active[row];
            for (let column = 0; column < count; column += 1) {
                kkt[row][column] = 2 * matrix[asset][active[column]];
            }
            kkt[row][count] = 1;
            kkt[row][count + 1] = -muScaled[asset];
        }
        for (let column = 0; column < count; column += 1) {
            const asset = active[column];
            kkt[count][column] = 1;
            kkt[count + 1][column] = muScaled[asset];
        }
        const solved = solveLinearSystem(kkt, [...Array(count).fill(0), 1, targetScaled]);
        if (solved === null) {
            break;
        }
        const faceWeights = solved.slice(0, count);
        const negatives = faceWeights
            .map((value, index) => value < -TOL_FEASIBILITY ? index : -1)
            .filter((index) => index >= 0);
        if (negatives.length > 0) {
            let removed = false;
            const ordered = [...negatives].sort((left, right) => faceWeights[left] - faceWeights[right]);
            for (const removePosition of ordered) {
                const remaining = active.filter((_, position) => position !== removePosition);
                const remainingValues = remaining.map((index) => muScaled[index]);
                if (remainingValues.length > 0 && Math.min(...remainingValues) <= targetScaled + TOL_RETURN_SCALED && Math.max(...remainingValues) >= targetScaled - TOL_RETURN_SCALED) {
                    active = remaining;
                    removed = true;
                    break;
                }
            }
            if (removed) {
                continue;
            }
            break;
        }
        const weights = Array(size).fill(0) as number[];
        active.forEach((index, position) => { weights[index] = Math.max(0, faceWeights[position]); });
        const total = weights.reduce((sum, value) => sum + value, 0);
        if (!Number.isFinite(total) || total <= 0) {
            break;
        }
        for (let index = 0; index < size; index += 1) {
            weights[index] /= total;
        }
        const gradient = matVec(matrix, weights).map((value) => 2 * value);
        const eta = solved[count + 1];
        const stationarity = Array.from({ length: size }, (_, index) => gradient[index] + solved[count] - eta * muScaled[index]);
        const positive = weights
            .map((value, index) => value > TOL_ACTIVE ? index : -1)
            .filter((index) => index >= 0);
        const violating = Array.from({ length: size }, (_, index) => index)
            .filter((index) => !positive.includes(index) && stationarity[index] < -tolGap);
        if (eta < -tolGap) {
            break;
        }
        if (violating.length > 0) {
            const addIndex = violating.reduce((best, index) => stationarity[index] < stationarity[best] ? index : best, violating[0]);
            if (!active.includes(addIndex)) {
                active.push(addIndex);
                active.sort((left, right) => left - right);
                continue;
            }
        }
        return { weights, iterations: iteration, eta };
    }
    return { weights: null, iterations: Math.min(maxIterations, size + 1), eta: null };
}
function enumerateGmvSupports(matrix: readonly (readonly number[])[], tolGap: number): {
    weights: number[] | null;
    iterations: number;
} {
    const size = matrix.length;
    let bestWeights: number[] | null = null;
    let bestObjective = Number.POSITIVE_INFINITY;
    for (let mask = 1; mask < (1 << size); mask += 1) {
        const support = Array.from({ length: size }, (_, index) => index).filter((index) => (mask & (1 << index)) !== 0);
        const count = support.length;
        let weights: number[];
        let multiplier: number;
        if (count === 1) {
            weights = Array(size).fill(0) as number[];
            weights[support[0]] = 1;
            multiplier = -2 * matrix[support[0]][support[0]];
        }
        else {
            const kkt = Array.from({ length: count + 1 }, () => Array(count + 1).fill(0));
            for (let row = 0; row < count; row += 1) {
                for (let column = 0; column < count; column += 1) {
                    kkt[row][column] = 2 * matrix[support[row]][support[column]];
                }
                kkt[row][count] = 1;
            }
            for (let column = 0; column < count; column += 1) {
                kkt[count][column] = 1;
            }
            const solved = solveLinearSystem(kkt, [...Array(count).fill(0), 1]);
            if (solved === null) {
                continue;
            }
            weights = Array(size).fill(0) as number[];
            support.forEach((asset, position) => { weights[asset] = solved[position]; });
            multiplier = solved[count];
        }
        if (weights.some((value) => value < -TOL_FEASIBILITY)) {
            continue;
        }
        weights = weights.map((value) => Math.max(0, value));
        const total = weights.reduce((sum, value) => sum + value, 0);
        if (!Number.isFinite(total) || total <= 0) {
            continue;
        }
        weights = weights.map((value) => value / total);
        const gradient = matVec(matrix, weights).map((value) => 2 * value);
        if (Array.from({ length: size }, (_, index) => index).some((index) => !support.includes(index) && gradient[index] + multiplier < -tolGap)) {
            continue;
        }
        const objective = quadratic(matrix, weights);
        if (objective < bestObjective || (objective === bestObjective && (bestWeights === null || lexicographicLess(weights, bestWeights)))) {
            bestObjective = objective;
            bestWeights = weights;
        }
    }
    return { weights: bestWeights, iterations: bestWeights === null ? 0 : 1 };
}
function enumerateTargetSupports(matrix: readonly (readonly number[])[], muScaled: readonly number[], targetScaled: number, tolGap: number): {
    weights: number[] | null;
    iterations: number;
    eta: number | null;
} {
    const size = matrix.length;
    let bestWeights: number[] | null = null;
    let bestObjective = Number.POSITIVE_INFINITY;
    let bestEta: number | null = null;
    for (let mask = 1; mask < (1 << size); mask += 1) {
        const support = Array.from({ length: size }, (_, index) => index).filter((index) => (mask & (1 << index)) !== 0);
        const count = support.length;
        let weights: number[];
        let multiplier: number;
        let eta: number;
        if (count === 1) {
            const asset = support[0];
            if (Math.abs(muScaled[asset] - targetScaled) > TOL_RETURN_SCALED) {
                continue;
            }
            weights = Array(size).fill(0) as number[];
            weights[asset] = 1;
            multiplier = -2 * matrix[asset][asset];
            eta = 0;
        }
        else {
            const kkt = Array.from({ length: count + 2 }, () => Array(count + 2).fill(0));
            for (let row = 0; row < count; row += 1) {
                const asset = support[row];
                for (let column = 0; column < count; column += 1) {
                    kkt[row][column] = 2 * matrix[asset][support[column]];
                }
                kkt[row][count] = 1;
                kkt[row][count + 1] = -muScaled[asset];
            }
            for (let column = 0; column < count; column += 1) {
                const asset = support[column];
                kkt[count][column] = 1;
                kkt[count + 1][column] = muScaled[asset];
            }
            const solved = solveLinearSystem(kkt, [...Array(count).fill(0), 1, targetScaled]);
            if (solved === null) {
                continue;
            }
            weights = Array(size).fill(0) as number[];
            support.forEach((asset, position) => { weights[asset] = solved[position]; });
            multiplier = solved[count];
            eta = solved[count + 1];
        }
        if (eta < -tolGap || weights.some((value) => value < -TOL_FEASIBILITY)) {
            continue;
        }
        weights = weights.map((value) => Math.max(0, value));
        const total = weights.reduce((sum, value) => sum + value, 0);
        if (!Number.isFinite(total) || total <= 0) {
            continue;
        }
        weights = weights.map((value) => value / total);
        if (Math.abs(dot(muScaled, weights) - targetScaled) > 10 * TOL_RETURN_SCALED) {
            continue;
        }
        const gradient = matVec(matrix, weights).map((value) => 2 * value);
        if (Array.from({ length: size }, (_, index) => index).some((index) => !support.includes(index) && gradient[index] + multiplier - eta * muScaled[index] < -tolGap)) {
            continue;
        }
        const objective = quadratic(matrix, weights);
        if (objective < bestObjective || (objective === bestObjective && (bestWeights === null || lexicographicLess(weights, bestWeights)))) {
            bestObjective = objective;
            bestWeights = weights;
            bestEta = eta;
        }
    }
    return { weights: bestWeights, iterations: bestWeights === null ? 0 : 1, eta: bestEta };
}
function frankWolfeSimplex(matrix: readonly (readonly number[])[], maxIterations: number, tolGap: number): {
    weights: number[] | null;
    iterations: number;
    gap: number | null;
} {
    const size = matrix.length;
    let weights = Array(size).fill(1 / size) as number[];
    for (let iteration = 1; iteration <= maxIterations; iteration += 1) {
        const gradient = matVec(matrix, weights).map((value) => 2 * value);
        const oracleIndex = Array.from({ length: size }, (_, index) => index).reduce((best, index) => gradient[index] < gradient[best] ? index : best, 0);
        const gap = dot(gradient, weights) - gradient[oracleIndex];
        if (gap <= tolGap) {
            return { weights, iterations: iteration, gap };
        }
        const oracle = Array(size).fill(0) as number[];
        oracle[oracleIndex] = 1;
        const direction = oracle.map((value, index) => value - weights[index]);
        const curvature = quadratic(matrix, direction);
        const linearChange = dot(gradient, direction);
        const denominator = 2 * curvature;
        let step: number;
        if (denominator > 0) {
            step = Math.min(1, Math.max(0, -linearChange / denominator));
        }
        else if (linearChange < 0) {
            step = 1;
        }
        else {
            return { weights, iterations: iteration, gap };
        }
        weights = weights.map((value, index) => value + step * direction[index]);
    }
    const gradient = matVec(matrix, weights).map((value) => 2 * value);
    const oracleIndex = Array.from({ length: size }, (_, index) => index).reduce((best, index) => gradient[index] < gradient[best] ? index : best, 0);
    return { weights, iterations: maxIterations, gap: dot(gradient, weights) - gradient[oracleIndex] };
}
function frankWolfeTarget(matrix: readonly (readonly number[])[], muScaled: readonly number[], targetScaled: number, maxIterations: number, tolGap: number): {
    weights: number[] | null;
    iterations: number;
    gap: number | null;
} {
    const candidates = feasibleLinearVertices(muScaled, targetScaled, TOL_RETURN_SCALED);
    if (candidates.length === 0) {
        return { weights: null, iterations: 0, gap: null };
    }
    let weights = candidates.reduce((best, candidate) => quadratic(matrix, candidate) < quadratic(matrix, best) ? candidate : best, candidates[0]);
    for (let iteration = 1; iteration <= maxIterations; iteration += 1) {
        const gradient = matVec(matrix, weights).map((value) => 2 * value);
        const oracle = linearMinimum(gradient, muScaled, targetScaled, TOL_RETURN_SCALED);
        if (oracle === null) {
            return { weights: null, iterations: iteration, gap: null };
        }
        const gap = dot(gradient, weights) - oracle.value;
        if (gap <= tolGap) {
            return { weights, iterations: iteration, gap };
        }
        const direction = oracle.point.map((value, index) => value - weights[index]);
        const curvature = quadratic(matrix, direction);
        const linearChange = dot(gradient, direction);
        const denominator = 2 * curvature;
        let step: number;
        if (denominator > 0) {
            step = Math.min(1, Math.max(0, -linearChange / denominator));
        }
        else if (linearChange < 0) {
            step = 1;
        }
        else {
            return { weights, iterations: iteration, gap };
        }
        weights = weights.map((value, index) => value + step * direction[index]);
    }
    const gradient = matVec(matrix, weights).map((value) => 2 * value);
    const oracle = linearMinimum(gradient, muScaled, targetScaled, TOL_RETURN_SCALED);
    return { weights, iterations: maxIterations, gap: oracle === null ? null : dot(gradient, weights) - oracle.value };
}
function failure(status: MarkowitzStatus, message: string, options: {
    assetIds?: string[] | null;
    mu?: number[] | null;
    targetReturn?: number | null;
    maxIterations?: number;
    diagnostics?: Record<string, unknown>;
} = {}): MarkowitzResult {
    return {
        assetIds: options.assetIds ?? null,
        mu: options.mu ?? null,
        targetReturn: options.targetReturn ?? null,
        weights: null,
        expectedReturn: null,
        variance: null,
        volatility: null,
        budgetResidual: null,
        targetResidual: null,
        lowerBoundResidual: null,
        fwGapS: null,
        fwGap: null,
        matrixScale: null,
        status,
        solutionClass: "unknown",
        variant: "markowitz-target-return-min-variance-long-only",
        method: "active-set-kkt-face-enumeration",
        iterations: 0,
        maxIterations: options.maxIterations ?? DEFAULT_MAX_ITERATIONS,
        diagnostics: { message, ...(options.diagnostics ?? {}) },
        warnings: [],
    };
}
function result(assetIds: string[], mu: number[], targetReturn: number, weights: number[], variance: number, iterations: number, maxIterations: number, scale: number, fwGapS: number | null, fwGap: number | null, status: MarkowitzStatus, className: MarkowitzSolutionClass, warnings: string[], diagnostics: Record<string, unknown>): MarkowitzResult {
    const expectedReturn = dot(mu, weights);
    const varianceTolerance = scale > 0 ? TOL_VARIANCE * scale : TOL_VARIANCE;
    const cleanVariance = variance < 0 && variance >= -varianceTolerance ? 0 : variance;
    return {
        assetIds,
        mu,
        targetReturn,
        weights,
        expectedReturn,
        variance: cleanVariance >= 0 ? cleanVariance : null,
        volatility: cleanVariance >= 0 ? Math.sqrt(cleanVariance) : null,
        budgetResidual: Math.abs(weights.reduce((sum, value) => sum + value, 0) - 1),
        targetResidual: Math.max(0, targetReturn - expectedReturn),
        lowerBoundResidual: Math.max(0, -Math.min(...weights)),
        fwGapS,
        fwGap,
        matrixScale: scale,
        status,
        solutionClass: className,
        variant: "markowitz-target-return-min-variance-long-only",
        method: "active-set-kkt-face-enumeration",
        iterations,
        maxIterations,
        diagnostics,
        warnings,
    };
}
export function markowitzMeanVariance(assetIdsInput: unknown, muInput: unknown, covarianceInput: unknown, targetReturnInput: unknown, options: MarkowitzOptions = {}): MarkowitzResult {
    if (!Array.isArray(assetIdsInput)) {
        return failure("invalid_input", "assetIds must be an ordered array.");
    }
    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 (assetIds.length === 0) {
        return failure("invalid_input", "assetIds must not be empty.");
    }
    if (new Set(assetIds).size !== assetIds.length) {
        return failure("invalid_input", "assetIds must be unique.");
    }
    if (!Array.isArray(muInput) || muInput.length !== assetIds.length || muInput.some((value) => !isFiniteReal(value))) {
        return failure("invalid_input", "mu must be a finite vector aligned to assetIds.", { assetIds });
    }
    const mu = muInput.map((value) => value as number);
    if (!isFiniteReal(targetReturnInput)) {
        return failure("invalid_input", "targetReturn must be a finite real number.", { assetIds, mu });
    }
    const targetReturn = targetReturnInput;
    const maxIterationsInput = options.maxIterations ?? DEFAULT_MAX_ITERATIONS;
    if (!Number.isInteger(maxIterationsInput) || (maxIterationsInput as number) < 0) {
        return failure("invalid_input", "maxIterations must be a non-negative integer.", { assetIds, mu, targetReturn });
    }
    const maxIterations = maxIterationsInput as number;
    if (maxIterations > DEFAULT_MAX_ITERATIONS) {
        return failure("invalid_input", "maxIterations may not exceed the fixed 10,000-iteration budget.", { assetIds, mu, targetReturn, maxIterations });
    }
    if (!Array.isArray(covarianceInput)) {
        return failure("invalid_input", "covariance must be a square numeric matrix.", { assetIds, mu, targetReturn, maxIterations });
    }
    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, mu, targetReturn, maxIterations });
    }
    const matrix: number[][] = [];
    for (const row of covarianceInput) {
        if (!Array.isArray(row) || row.some((value) => !isFiniteReal(value))) {
            return failure("invalid_input", "covariance entries must be finite real numbers.", { assetIds, mu, targetReturn, maxIterations });
        }
        matrix.push(row.map((value) => value as number));
    }
    const scale = maxAbsMatrix(matrix);
    let returnScale = Math.max(...mu.map((value) => Math.abs(value)), Math.abs(targetReturn));
    if (returnScale === 0) {
        returnScale = 1;
    }
    const muScaled = mu.map((value) => value / returnScale);
    const targetScaled = targetReturn / returnScale;
    if (scale === 0) {
        if (targetScaled > Math.max(...muScaled) + TOL_RETURN_SCALED) {
            return failure("infeasible", "targetReturn exceeds the maximum feasible simplex expected return.", { assetIds, mu, targetReturn, maxIterations, diagnostics: { matrixScale: 0, returnScale, maxAssetReturn: Math.max(...mu) } });
        }
        const candidates = feasibleLinearVertices(muScaled, targetScaled, TOL_RETURN_SCALED);
        let weights: number[] = candidates.length > 0 ? candidates[0] : Array(size).fill(1 / size) as number[];
        for (const candidate of candidates.slice(1)) {
            if (lexicographicLess(candidate, weights)) {
                weights = candidate;
            }
        }
        const expected = dot(mu, weights);
        return result(assetIds, mu, targetReturn, weights, 0, 0, maxIterations, 0, 0, 0, "optimal", classifySolution(Array.from({ length: size }, () => Array(size).fill(0)), muScaled, targetScaled, weights), ["all_zero_covariance_flat_objective"], {
            matrixScale: 0,
            returnScale,
            normalizedObjective: 0,
            targetBinding: Math.abs(expected / returnScale - targetScaled) <= TOL_RETURN_SCALED,
            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, targetReturn, maxIterations, diagnostics: { 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 baseDiagnostics: Record<string, unknown> = {
        matrixScale: scale,
        symmetryResidualS: symmetryResidual,
        minEigenvalueS: minimumEigenvalue,
        maxEigenvalueS: maximumEigenvalue,
        eigenResidualS: eigensystem.residual,
        tolPsd: TOL_PSD,
        tolGapRelative: TOL_GAP_RELATIVE,
    };
    if (eigensystem.residual > TOL_EIGEN_RESIDUAL) {
        return failure("numerical_issue", "symmetric eigensolver residual exceeds the fixed normalized tolerance.", { assetIds, mu, targetReturn, maxIterations, diagnostics: baseDiagnostics });
    }
    if (minimumEigenvalue < -TOL_PSD) {
        return failure("invalid_input", "normalized covariance is materially indefinite.", { assetIds, mu, targetReturn, maxIterations, diagnostics: baseDiagnostics });
    }
    const tolReturn = TOL_RETURN_SCALED;
    const tolGap = TOL_GAP_RELATIVE * Math.max(1, rowSumNorm(normalized));
    if (targetScaled > Math.max(...muScaled) + TOL_RETURN_SCALED) {
        return failure("infeasible", "targetReturn exceeds the maximum feasible simplex expected return.", { assetIds, mu, targetReturn, maxIterations, diagnostics: { ...baseDiagnostics, returnScale, maxAssetReturn: Math.max(...mu), tolGap } });
    }
    const equalWeights = Array(size).fill(1 / size) as number[];
    if (maxIterations === 0) {
        const gradient = matVec(normalized, equalWeights).map((value) => 2 * value);
        const oracle = linearMinimum(gradient, muScaled, targetScaled, tolReturn);
        const gap = oracle === null ? null : dot(gradient, equalWeights) - oracle.value;
        return result(assetIds, mu, targetReturn, equalWeights, quadratic(normalized, equalWeights) * scale, 0, maxIterations, scale, gap, gap === null ? null : gap * scale, "numerical_issue", classifySolution(normalized, muScaled, targetScaled, equalWeights), [], { ...baseDiagnostics, returnScale, tolGap, optimalityCertificate: "not-certified-max-iterations-zero" });
    }
    const warnings: string[] = [];
    if (symmetryResidual > 0) {
        warnings.push("symmetrized_within_tolerance");
    }
    if (minimumEigenvalue < 0) {
        warnings.push("psd_within_roundoff");
    }
    let iterations = 0;
    const gmv = gmvFace(normalized, maxIterations, tolGap);
    let gmvWeights = gmv.weights;
    let gmvIterations = gmv.iterations;
    if (gmvWeights === null && size <= MAX_EXACT_SUPPORT_ASSETS) {
        const enumerated = enumerateGmvSupports(normalized, tolGap);
        gmvWeights = enumerated.weights;
        gmvIterations = enumerated.iterations;
    }
    if (gmvWeights === null) {
        const fallback = frankWolfeSimplex(normalized, maxIterations, tolGap);
        gmvWeights = fallback.weights;
        gmvIterations = fallback.iterations;
    }
    if (gmvWeights === null) {
        return failure("numerical_issue", "simplex GMV solve did not produce a feasible candidate.", { assetIds, mu, targetReturn, maxIterations, diagnostics: { ...baseDiagnostics, returnScale, tolGap } });
    }
    iterations += gmvIterations;
    const gmvExpected = dot(muScaled, gmvWeights);
    let selectedWeights: number[] | null = null;
    let targetBinding = false;
    let eta: number | null = null;
    if (gmvExpected >= targetScaled - TOL_RETURN_SCALED) {
        selectedWeights = gmvWeights;
    }
    else {
        targetBinding = true;
        if (Math.abs(targetScaled - Math.max(...muScaled)) <= TOL_RETURN_SCALED) {
            const maximum = Math.max(...mu);
            const eligible = muScaled.map((value, index) => Math.abs(value - Math.max(...muScaled)) <= TOL_RETURN_SCALED ? index : -1).filter((index) => index >= 0);
            const submatrix = eligible.map((row) => eligible.map((column) => normalized[row][column]));
            const remaining = Math.max(0, maxIterations - iterations);
            const sub = gmvFace(submatrix, remaining, tolGap);
            let subWeights = sub.weights;
            let subIterations = sub.iterations;
            if (subWeights === null && remaining > 0 && eligible.length <= MAX_EXACT_SUPPORT_ASSETS) {
                const enumerated = enumerateGmvSupports(submatrix, tolGap);
                subWeights = enumerated.weights;
                subIterations = enumerated.iterations;
            }
            if (subWeights !== null) {
                selectedWeights = Array(size).fill(0) as number[];
                eligible.forEach((index, position) => { (selectedWeights as number[])[index] = subWeights?.[position] ?? 0; });
                iterations += Math.min(subIterations, remaining);
            }
        }
        if (selectedWeights === null) {
            const remaining = Math.max(0, maxIterations - iterations);
            if (remaining > 0) {
                const target = targetFace(normalized, muScaled, targetScaled, remaining, tolGap);
                selectedWeights = target.weights;
                eta = target.eta;
                iterations += Math.min(target.iterations, remaining);
                if (selectedWeights === null && size <= MAX_EXACT_SUPPORT_ASSETS) {
                    const enumerated = enumerateTargetSupports(normalized, muScaled, targetScaled, tolGap);
                    selectedWeights = enumerated.weights;
                    eta = enumerated.eta;
                    iterations += Math.min(enumerated.iterations, Math.max(0, maxIterations - iterations));
                }
                if (selectedWeights === null && size > MAX_EXACT_SUPPORT_ASSETS) {
                    const fallback = frankWolfeTarget(normalized, muScaled, targetScaled, Math.max(0, maxIterations - iterations), tolGap);
                    selectedWeights = fallback.weights;
                    iterations += Math.min(fallback.iterations, Math.max(0, maxIterations - iterations));
                }
            }
        }
    }
    if (selectedWeights === null) {
        return failure("numerical_issue", "active-set KKT solve did not produce a certified feasible point.", { assetIds, mu, targetReturn, maxIterations, diagnostics: { ...baseDiagnostics, returnScale, tolGap, targetBinding } });
    }
    const expectedReturn = dot(mu, selectedWeights);
    const normalizedObjective = quadratic(normalized, selectedWeights);
    const gradient = matVec(normalized, selectedWeights).map((value) => 2 * value);
    const oracle = linearMinimum(gradient, muScaled, targetScaled, tolReturn);
    const gap = oracle === null ? null : dot(gradient, selectedWeights) - oracle.value;
    const budgetResidual = Math.abs(selectedWeights.reduce((sum, value) => sum + value, 0) - 1);
    const targetResidual = Math.max(0, targetReturn - expectedReturn);
    const targetResidualScaled = Math.max(0, targetScaled - expectedReturn / returnScale);
    const lowerResidual = Math.max(0, -Math.min(...selectedWeights));
    const status: MarkowitzStatus = gap !== null && gap <= tolGap && budgetResidual <= TOL_FEASIBILITY && targetResidualScaled <= TOL_RETURN_SCALED && lowerResidual <= TOL_FEASIBILITY && normalizedObjective >= -TOL_VARIANCE ? "optimal" : "numerical_issue";
    const certificate = status === "optimal" && minimumEigenvalue >= 0 ? "simplex-target-return-frank-wolfe-convex-certificate" : status === "optimal" ? "tolerance-limited-near-psd" : "not-certified";
    return result(assetIds, mu, targetReturn, selectedWeights, scale * normalizedObjective, Math.min(iterations, maxIterations), maxIterations, scale, gap, gap === null ? null : gap * scale, status, classifySolution(normalized, muScaled, targetScaled, selectedWeights), warnings, { ...baseDiagnostics, returnScale, tolReturnScaled: TOL_RETURN_SCALED, tolGap, normalizedObjective, varianceRaw: scale * normalizedObjective, targetBinding: Math.abs(expectedReturn / returnScale - targetScaled) <= TOL_RETURN_SCALED, targetMultiplierScaled: eta, optimalityCertificate: certificate });
}
export const globalMinimumVarianceTarget = markowitzMeanVariance;

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

Islam Baraka

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

Markowitz Mean-Variance: When a Return Floor Adds Risk

Before you start

What “Markowitz” means here

The input is part of the algorithm

A reproducible two-asset fixture

Boundaries worth testing

What the covariance gate checks

A prior-method lesson: stale observations can manufacture “riskless” weights

Implementation and independent checks

Reproduce, then challenge the answer

Check your understanding

The calculation path

Summary

Decision support: what should change the answer?

Small glossary

Continue through Family 01

References

Production handoff checklist

Primary technical sources

R01 — University of Washington, “Markowitz Mean-Variance Portfolio Theory”

R02 — MOSEK, “Markowitz portfolio optimization,” Portfolio Optimization Cookbook §2

R03 — Sharma et al., “Decision-Focused Learning for Mean-Variance Portfolio Optimization”

R04 — Harry Markowitz, “Portfolio Selection”

Input-lineage sources

R05 — U.S. SEC, “Stock Splits”

R06 — Investor.gov, “Ex-Dividend Dates: When Are You Entitled to Stock and Its Dividend?”

Contextual Google destinations inspected

Citation boundaries

Independent recheck — 2026-09-16

Islam Baraka

Related released builds.

Global Minimum Variance

Maximum Sharpe Ratio

Mean-CVaR Optimization