By Islam Baraka“Keep at most two assets” sounds precise. Start with four existing holdings, though, and reaching that two-asset portfolio can require four trades. There is a second trap: selecting the two highest estimated returns need not produce the best risk-adjusted pair.
Cardinality is a count. A portfolio's support is the set of assets with nonzero weights. This article limits the cardinality of the final support. You will learn to distinguish it from a trade count, solve one two-asset subproblem, and compare all ten candidates in a small exact search. Familiarity with mean-variance optimization is enough; the binary variables are introduced below.
One portfolio, two different counts
Consider this synthetic rebalance:
current = [.25, .25, .25, .25]
final = [.60, .40, .00, .00]
trade = [.35, .15, −.25, −.25].
Two final weights are positive, but every asset changes. A rule allowing at most two final holdings accepts this candidate. A rule allowing at most two trades rejects it. Current holdings affect the second count; the final vector alone determines the first.
Count the nonzero entries in the final row and then in the trade row: two versus four. Open the count comparison.
In my portfolio-data work, changing weights and suspect imported quantities have made the starting snapshot worth checking before reasoning about a rebalance. Here, I also check stable asset identity. Two imported aliases for one security must not become two economic holdings merely because they occupy two rows. This is a data-validation rule; it does not merge genuinely different share classes or instruments.
Turn “at most K holdings” into constraints
Let x_i be asset i's final portfolio weight and y_i be a binary switch: zero switches the holding off; one permits a holding within its bounds. K is the maximum number of switches we may use.
| Input | Meaning |
|---|---|
μ, Σ | Expected simple returns and return covariance on one common horizon |
λ≥0 | Chosen coefficient balancing expected return and variance |
l_i, u_i | Minimum and maximum weight when asset i is selected |
K | Integer maximum count of final holdings |
The model is
maximize μᵀx − λ xᵀΣx
subject to Σ_i x_i = 1
l_i y_i ≤ x_i ≤ u_i y_i
Σ_i y_i ≤ K
x_i ≥ 0, y_i ∈ {0,1}.
If y_i=0, both linked bounds force x_i=0. If y_i=1 and the selected bounds are [.05,.60], the holding must occupy between 5% and 60% of the portfolio. A positive minimum is often called a buy-in threshold here; it applies to the final holding, not automatically to its trade.
There is a subtlety when l_i=0: y_i=1 still permits x_i=0. The switch says “allowed,” not “definitely held.” The at-most-K model remains valid, but the reported holding count must come from the actual positive x_i. Likewise, replacing Σy_i≤K by Σy_i=K does not guarantee exactly K positive holdings unless activation is tied to a positive minimum. Rounding tiny positions out of the display is not a substitute.
MOSEK's cardinality example uses this binary-linking idea on trades, x−x⁰. Applying it to final x is the declared adaptation here. The reference checks that adaptation by enumerating final supports.
Why top-K is not the optimizer
For the rest of the article, use this synthetic fixture, with selected lower bounds zero and upper bounds one:
μ = [0.12, 0.10, 0.09, 0.08]
Σ = diag(0.30, 0.02, 0.01, 0.01)
λ = 1
K = 2.
A has the highest expected return, but its variance is also much higher. Compare three decisions:
| Decision | Weights [A,B,C,D] | Utility |
|---|---|---|
| Highest two means, equal weights | [.5,.5,0,0] | .0300000 |
| Highest two means, best weights within that pair | [.09375,.90625,0,0] | .0828125 |
| Best pair across the whole universe | [0,.5,.5,0] | .0875000 |
Even optimizing the weights after screening to A and B misses the best support. In this declared objective, the low variance of C is worth more than A's higher mean. The covariance matrix is diagonal in this teaching case, so there are no cross-covariance terms; a general problem must include them too.
The winning pair emerges from comparing optimized supports, not from ranking means. Open the support comparison.
The implementation visits all four singleton supports and all six pairs. For each support it solves the continuous bounded quadratic subproblem. A free coordinate lies strictly between its bounds. The Karush–Kuhn–Tucker (KKT) optimality conditions require its marginal utility to equal a common budget multiplier:
μ_i − 2λ(Σx)_i = ν,
where ν is the multiplier for the full-investment equality. Lower- and upper-bound coordinates are checked with the corresponding one-sided inequalities. The best certified candidate is selected only after all supports are visited.
Walk one support before trusting the search
For the support {B,C} in the synthetic fixture, set A and D to zero and write x_B=b, x_C=1−b. The restricted objective is
f(b) = .10b + .09(1−b) − .02b² − .01(1−b)²
= .08 + .03b − .03b².
Its derivative, .03−.06b, vanishes at b=.5; the concave quadratic therefore gives [0,.5,.5,0] and f=.0875. This is one continuous subproblem, not yet the cardinality proof. The proof for the small reference comes from comparing it with every allowed singleton and pair, including their bounds and KKT conditions. With thousands of assets, that enumeration is combinatorial; a mixed-integer solver or a labeled heuristic is a different engineering choice.
There is a useful feasibility precheck. If K=1 but every asset's upper bound is below one, no one-asset portfolio can satisfy Σx=1. Conversely, if a positive lower bound applies to each selected holding, choosing too many holdings can force their minimum total above one. A solver should report these conflicts explicitly. K must be an integer in [1,n]. With full investment, K=0 is mathematically impossible; this API rejects it as an invalid input before searching.
Compare every allowed support
For any pair {i,j} in this diagonal fixture, write its weights as [a,1−a]. If their variances are v_i and v_j, the free optimum is
a = (μ_i − μ_j + 2v_j) / [2(v_i + v_j)],
clipped to [0,1] for these zero-to-one bounds. The following table applies that formula independently; weights appear in the support's listed order. Values are rounded only for display.
| Allowed support | Best weights within it | Best utility |
|---|---|---|
| A | [1] | −.180000000 |
| B | [1] | .080000000 |
| C | [1] | .080000000 |
| D | [1] | .070000000 |
| A, B | [.093750,.906250] | .082812500 |
| A, C | [.080645,.919355] | .082016129 |
| A, D | [.096774,.903226] | .072903226 |
| B, C | [.500000,.500000] | .087500000 |
| B, D | [.666667,.333333] | .083333333 |
| C, D | [.750000,.250000] | .081250000 |
Now the small-instance global claim is inspectable: B/C beats every allowed singleton and pair. “At most two” required checking the singletons as well. If several candidates tie, an implementation should make its tie-breaking rule explicit rather than suggest a unique economic answer.
Honest solver statuses
The reference is intentionally small (n≤12, K≤6). The number of supports is Σ_{j=1..K} choose(n,j); every support can also require several bound configurations. Increasing the universe quickly makes exhaustive search expensive.
| Status | What the reader may conclude |
|---|---|
optimal | All allowed supports and active-set candidates were checked within stated tolerances. |
feasible | A candidate exists, but max_supports stopped the proof; no best bound is invented. |
incomplete | Search stopped before finding a candidate; this does not prove infeasibility. |
infeasible | The completed search found no support satisfying the declared bounds. |
An incumbent is simply the best feasible candidate found so far. Its existence and its optimality are different facts.
The code also rejects a nonsymmetric or indefinite covariance at the matrix’s actual scale, invalid bounds, unnormalized current holdings, and out-of-scope dimensions. A singular covariance with positive risk aversion is not silently regularized; use a separately certified mixed-integer/QP workflow for that case.
Where a sparse answer can mislead
- Support is defined on final holdings. Rounding a tiny positive position to zero for display must not quietly change the reported support or the budget. If an execution system imposes a minimum trade or share increment, add that constraint and solve again.
- A top-K screen is only a heuristic. It can discard a low-mean asset that diversifies a high-risk asset. The worked example shows exactly that: C enters the optimal pair even though A and B have higher means.
- A sparse result may require many trades. The
[.60,.40,0,0]candidate needs four trades from the stated current portfolio. The turnover tutorial limits trade size; the cost tutorial funds those trades. Neither is implied byK=2. - A universe is dated data. Preserve one stable asset ID per economic holding and the universe available at the decision time. Duplicate aliases and later-listed assets can make a historical sparse portfolio look feasible when it was not.
For a production mixed-integer solve, retain incumbent, bound, gap, tolerance, and termination reason. “Feasible” after a timeout means a candidate passed the checks; it does not establish that no better support exists. After solving, recompute Σx, each selected bound, realized support, and the stated utility from the returned vector. This direct audit is the small habit that keeps a polished sparse allocation from becoming an unsupported optimality claim.
Try the playground
Open the interactive cardinality lab. Start by predicting these three outcomes:
- For the canonical K=2 final-holdings model, the winner is B/C with
.0875utility. Count its final holdings and its changes from[.25,.25,.25,.25]: two holdings, four trades. - If only one holding is allowed, B and C tie at
.08. The smaller count changes the attainable objective, not merely the labels on the same bars. - If enumeration stops early, even a plausible allocation must remain an incumbent. Inspect the partial-search scenario's status before comparing its objective with the complete search.
Switching to trade-support mode changes the question being solved. Use the displayed current/final/trade rows to verify what the count means in each mode.
Reproduce the fixture
From the topic folder:
python -B implementations/python/cardinality_constrained.py
python -B -m unittest discover -s tests -p 'test_*.py' -v
npx --yes tsx tests/cardinality-constrained.test.ts
What I want the reader to keep
I would write the policy in one sentence before opening a solver: “at most K final economic holdings,” or “at most K nonzero trades.” Then I would attach the binary variable to that exact object, deduplicate the dated asset universe, and check the returned support directly. The four-trade/two-holding example shows why the noun after K matters as much as the number. For a different operational question, continue to the gross/net exposure tutorial; it counts signed size, not holdings.
All examples are synthetic. A named account case remains deferred pending a dated asset universe, reliable positions and prices, and publication rights. The canonical definition and source record document the numerical limits and source adaptation.
ReferencesPrimary sources and evidence notesExpand the source trail, evidence role, and limitations behind the engineering choices.
Expand the source trail, evidence role, and limitations behind the engineering choices.
MOSEK Portfolio Optimization Cookbook, §6.4 and §6.5
- Organization or authors: MOSEK ApS.
- Source type: Official solver documentation and portfolio-optimization cookbook.
- Publication or effective date: Online version accessed 2026-09-14.
- Version: Portfolio Optimization Cookbook 1.6.0, transaction-cost chapter.
- URL or DOI: MOSEK transaction costs and cardinality
- Accessed: 2026-09-14.
- Jurisdiction: None; this is mathematical/solver documentation.
- Supports: the source's §6.4 binary linking model for limiting nonzero coordinates of the trade vector, explicit upper bounds on trades, and §6.5 semi-continuous buy-in thresholds.
- Limitations: MOSEK's §6.4 variable is tilde-x=x-x0, so its cardinality is trade cardinality, not final-holdings cardinality. This package adapts the binary-link idea to final x_i with l_i y_i ≤ x_i ≤ u_i y_i; that adaptation is separately stated and checked by exhaustive support enumeration. The source does not certify this package's exact numerical fixture or a universal production-scale MIP solver.
Moehle, Kochenderfer, Boyd, and Ang, “Tax-aware portfolio optimization”
- Organization or authors: Nicholas Moehle, Mykel J. Kochenderfer, Stephen Boyd, and Andrew Ang.
- Source type: Original research paper/preprint.
- Publication or effective date: Journal version 2021; arXiv manuscript revised February 2021.
- Version: DOI 10.1007/s10957-021-01823-0; arXiv 2008.04985.
- URL or DOI: Stanford paper page; arXiv PDF; DOI.
- Accessed: 2026-09-14.
- Jurisdiction: None for the cardinality contract; tax-specific policy is out of scope here.
- Supports: integer/binary practical constraints can make otherwise convex portfolio optimization nonconvex and require explicit mixed-integer treatment.
- Limitations: this paper is contextual evidence about nonconvex portfolio constraints, not a direct definition of final-holdings cardinality. It does not supply the package fixture or validate the trade/holding adaptation.
Google discovery status
The A03 Google top-ten ledger remains queued in the bounded wave. When run, it must record the exact query, locale, date/time, Google surface, observed order and page-review extent; snippets, AI answers, and other engines cannot be relabeled as Google organic results.
Full dependency-light reference implementations in both supported languages.
/**
* Exact small-cardinality long-only portfolio reference.
*
* The binary support variable marks final holdings. It is deliberately not
* a trade-count variable. For the small declared scope, every support is
* enumerated and its bounded concave quadratic subproblem is solved by
* bound-active KKT enumeration.
*/
export const ABSOLUTE_TOLERANCE = 1e-9;
export const RELATIVE_TOLERANCE = 2e-8;
export const FEASIBILITY_TOLERANCE = 2e-9;
export const KKT_TOLERANCE = 2e-7;
export const BOUND_TOLERANCE = 2e-8;
export const SUPPORT_TOLERANCE = 1e-8;
export const PSD_RELATIVE_TOLERANCE = 1e-10;
export const MAX_ASSETS = 12;
export const MAX_CARDINALITY = 6;
export class CardinalityValidationError extends Error {
readonly code: string;
constructor(code: string, message: string) {
super(message);
this.name = "CardinalityValidationError";
this.code = code;
}
}
export class CardinalityNumericalError extends Error {
readonly code: string;
constructor(code: string, message: string) {
super(message);
this.name = "CardinalityNumericalError";
this.code = code;
}
}
export interface CardinalityCertificate {
kind: "exhaustive-support-enumeration";
complete: boolean;
continuousKkt?: boolean;
globalOptimality: string;
}
export interface CardinalityResult {
status: "optimal" | "feasible" | "incomplete" | "infeasible";
weights: number[] | null;
selected: number[] | null;
trade: number[] | null;
objective: number | null;
bestBound: number | null;
relativeGap: number | null;
enumeratedSupports: number;
totalSupports: number;
realizedSupport?: number;
tradeSupport?: number;
certificate: CardinalityCertificate;
}
interface Validated {
mean: number[]; sigma: number[][]; risk: number; current: number[];
cardinality: number; lower: number[]; upper: number[];
}
interface BestPoint {
value: number;
point: number[];
}
function invalid(code: string, message: string): never {
throw new CardinalityValidationError(code, message);
}
function asReal(value: unknown, name: string): number {
if (typeof value !== "number" || !Number.isFinite(value)) {
invalid("invalid-input", name + " must be a finite real number.");
}
return value;
}
function vector(value: unknown, name: string): number[] {
if (!Array.isArray(value) || value.length === 0) {
invalid("invalid-input", name + " must be a non-empty numeric array.");
}
return value.map((entry, index) => asReal(entry, name + "[" + index + "]"));
}
function matrix(value: unknown, name: string, size: number): number[][] {
if (!Array.isArray(value) || value.length !== size) {
invalid("invalid-input", name + " must have " + size + " rows.");
}
return value.map((row, rowIndex) => {
if (!Array.isArray(row) || row.length !== size) {
invalid("invalid-input", name + " must be " + size + " by " + size + ".");
}
return row.map((entry, columnIndex) => asReal(entry, name + "[" + rowIndex + "][" + columnIndex + "]"));
});
}
function stableSum(values: readonly number[]): number {
let total = 0;
let compensation = 0;
for (const value of values) {
const candidate = total + value;
if (Math.abs(total) >= Math.abs(value)) compensation += (total - candidate) + value;
else compensation += (value - candidate) + total;
total = candidate;
}
return total + compensation;
}
function dot(left: readonly number[], right: readonly number[]): number {
return stableSum(left.map((value, index) => value * right[index]));
}
function objective(point: readonly number[], mean: readonly number[], covariance: readonly (readonly number[])[], risk: number): number {
return dot(mean, point) - risk * dot(point, covariance.map((row) => dot(row, point)));
}
function minimumEigenvalue(matrixValue: readonly (readonly number[])[]): number {
const size = matrixValue.length;
const work = matrixValue.map((row) => [...row]);
if (size === 1) return work[0][0];
const scale = Math.max(...work.flat().map((value) => Math.abs(value)));
if (scale === 0) return 0;
for (let sweep = 0; sweep < Math.max(60, 20 * size * size); sweep += 1) {
let p = 0; let q = 1; let largest = 0;
for (let row = 0; row < size; row += 1) {
for (let column = row + 1; column < size; column += 1) {
const magnitude = Math.abs(work[row][column]);
if (magnitude > largest) { largest = magnitude; p = row; q = column; }
}
}
if (largest <= 1e-14 * scale) 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];
work[index][p] = cosine * aip - sine * aiq; work[p][index] = work[index][p];
work[index][q] = sine * aip + cosine * aiq; work[q][index] = work[index][q];
}
work[p][p] = cosine * cosine * app - 2 * cosine * sine * apq + sine * sine * aqq;
work[q][q] = sine * sine * app + 2 * cosine * sine * apq + cosine * cosine * aqq;
work[p][q] = 0; work[q][p] = 0;
}
return Math.min(...work.map((row, index) => row[index]));
}
function solveLinearSystem(matrixValue: readonly (readonly number[])[], right: readonly number[]): number[] | undefined {
const size = right.length;
const work = matrixValue.map((row, index) => [...row, right[index]]);
const scale = Math.max(...matrixValue.flat().map((value) => Math.abs(value)));
if (scale === 0) return undefined;
const pivotTolerance = 1e-12 * scale;
for (let column = 0; column < size; column += 1) {
let pivot = column;
for (let row = column + 1; row < size; row += 1) {
if (Math.abs(work[row][column]) > Math.abs(work[pivot][column])) pivot = row;
}
if (Math.abs(work[pivot][column]) <= pivotTolerance) return undefined;
if (pivot !== column) [work[column], work[pivot]] = [work[pivot], work[column]];
const pivotValue = work[column][column];
for (let row = column + 1; row < size; row += 1) {
const factor = work[row][column] / pivotValue;
if (factor === 0) continue;
for (let entry = column; entry <= size; entry += 1) work[row][entry] -= factor * work[column][entry];
}
}
const solution = Array.from({ length: size }, () => 0);
for (let row = size - 1; row >= 0; row -= 1) {
let remainder = work[row][size];
for (let entry = row + 1; entry < size; entry += 1) remainder -= work[row][entry] * solution[entry];
solution[row] = remainder / work[row][row];
}
return solution.every((value) => Number.isFinite(value)) ? solution : undefined;
}
function solveKktSystem(
hessian: readonly (readonly number[])[], linear: readonly number[], budget: number,
): number[] | undefined {
// The budget row has unit coefficients while the stationarity rows may be
// 1e-100 or 1e100. Scale stationarity by the actual objective magnitude,
// then restore the multiplier to original units after elimination.
const size = linear.length;
let stationarityScale = 0;
for (const row of hessian) for (const value of row) stationarityScale = Math.max(stationarityScale, Math.abs(value));
for (const value of linear) stationarityScale = Math.max(stationarityScale, Math.abs(value));
if (stationarityScale === 0 || !Number.isFinite(stationarityScale)) return undefined;
const system = Array.from({ length: size + 1 }, () => Array.from({ length: size + 1 }, () => 0));
const right = Array.from({ length: size + 1 }, () => 0);
for (let row = 0; row < size; row += 1) {
for (let column = 0; column < size; column += 1) system[row][column] = hessian[row][column] / stationarityScale;
system[row][size] = 1;
right[row] = linear[row] / stationarityScale;
}
for (let column = 0; column < size; column += 1) system[size][column] = 1;
right[size] = budget;
const solution = solveLinearSystem(system, right);
if (!solution) return undefined;
solution[size] *= stationarityScale;
return solution;
}
function validateInputs(
meanValue: unknown, covarianceValue: unknown, lambdaValue: unknown,
currentValue: unknown, cardinalityValue: unknown, lowerValue: unknown, upperValue: unknown,
): Validated {
const mean = vector(meanValue, "mu");
const size = mean.length;
const sigma = matrix(covarianceValue, "covariance", size);
const risk = asReal(lambdaValue, "lambda_risk");
const current = vector(currentValue, "current");
if (current.length !== size) invalid("invalid-input", "mu, covariance, and current must align.");
if (typeof cardinalityValue !== "number" || !Number.isInteger(cardinalityValue)) invalid("invalid-input", "K must be an integer.");
const cardinality = cardinalityValue;
if (size > MAX_ASSETS || cardinality > MAX_CARDINALITY) {
invalid("scope-limit", "the exact reference is limited to n≤" + MAX_ASSETS + " and K≤" + MAX_CARDINALITY + "; use a certified MIP solver for larger inputs.");
}
if (cardinality < 1 || cardinality > size) invalid("invalid-cardinality", "K must satisfy 1 ≤ K ≤ n.");
if (risk < 0) invalid("invalid-input", "lambda_risk must be nonnegative.");
if (current.some((value) => value < 0)) invalid("invalid-input", "current holdings must be nonnegative.");
const currentTotal = stableSum(current);
if (!Number.isFinite(currentTotal)) invalid("wealth-overflow", "current holdings sum must be finite.");
if (Math.abs(currentTotal - 1) > 10 * FEASIBILITY_TOLERANCE) invalid("current-not-normalized", "current holdings must sum to one.");
const lower = lowerValue === undefined || lowerValue === null ? Array.from({ length: size }, () => 0) : vector(lowerValue, "lower");
const upper = upperValue === undefined || upperValue === null ? Array.from({ length: size }, () => 1) : vector(upperValue, "upper");
if (lower.length !== size || upper.length !== size) invalid("invalid-input", "lower and upper bounds must align with mu.");
lower.forEach((value, index) => {
const high = upper[index];
if (value < 0 || high < 0 || value > high || high > 1 || !Number.isFinite(value + high)) {
invalid("invalid-bound", "lower[" + index + "] and upper[" + index + "] must satisfy 0 ≤ lower ≤ upper ≤ 1.");
}
});
const scale = Math.max(...sigma.flat().map((value) => Math.abs(value)));
if (scale === 0) {
if (risk > 0) invalid("singular-covariance", "positive risk aversion requires a positive-definite covariance in the exact reference scope.");
} else {
const symmetryTolerance = 1e-12 * scale;
for (let row = 0; row < size; row += 1) {
for (let column = 0; column < row; column += 1) {
if (Math.abs(sigma[row][column] - sigma[column][row]) > symmetryTolerance) invalid("indefinite-covariance", "covariance must be symmetric at the matrix scale.");
}
}
const minimum = minimumEigenvalue(sigma);
if (minimum < -PSD_RELATIVE_TOLERANCE * scale) invalid("indefinite-covariance", "covariance must be positive semidefinite.");
if (risk > 0 && minimum <= PSD_RELATIVE_TOLERANCE * scale) invalid("singular-covariance", "positive risk aversion requires a positive-definite covariance in the exact active-set reference.");
}
return { mean, sigma, risk, current, cardinality, lower, upper };
}
function pointIsBetter(
point: readonly number[], best: BestPoint | undefined,
mean: readonly number[], sigma: readonly (readonly number[])[], risk: number,
): BestPoint | undefined {
const value = objective(point, mean, sigma, risk);
if (!Number.isFinite(value)) throw new CardinalityNumericalError("nonfinite-objective", "continuous subproblem produced a non-finite objective.");
const comparisonTolerance = best === undefined ? 0 : RELATIVE_TOLERANCE * Math.max(Math.abs(value), Math.abs(best.value));
return best === undefined || value > best.value + comparisonTolerance ? { value, point: [...point] } : best;
}
function continuousLinearSupport(
mean: readonly number[], sigma: readonly (readonly number[])[], risk: number,
lower: readonly number[], upper: readonly number[], support: readonly number[],
): BestPoint | undefined {
let best: BestPoint | undefined;
const dimension = support.length;
for (let freePosition = 0; freePosition < dimension; freePosition += 1) {
const fixedPositions = Array.from({ length: dimension }, (_, position) => position).filter((position) => position !== freePosition);
for (let mask = 0; mask < (1 << fixedPositions.length); mask += 1) {
const point = Array.from({ length: mean.length }, () => 0);
fixedPositions.forEach((position, bit) => {
const index = support[position];
point[index] = mask & (1 << bit) ? upper[index] : lower[index];
});
const freeIndex = support[freePosition];
point[freeIndex] = 1 - stableSum(point);
if (point[freeIndex] < lower[freeIndex] - FEASIBILITY_TOLERANCE || point[freeIndex] > upper[freeIndex] + FEASIBILITY_TOLERANCE) continue;
point[freeIndex] = Math.min(upper[freeIndex], Math.max(lower[freeIndex], point[freeIndex]));
if (Math.abs(stableSum(point) - 1) > 10 * FEASIBILITY_TOLERANCE) continue;
best = pointIsBetter(point, best, mean, sigma, risk);
}
}
return best;
}
function continuousSupport(
mean: readonly number[], sigma: readonly (readonly number[])[], risk: number,
lower: readonly number[], upper: readonly number[], support: readonly number[],
): BestPoint | undefined {
const supportLower = stableSum(support.map((index) => lower[index]));
const supportUpper = stableSum(support.map((index) => upper[index]));
if (supportLower > 1 + FEASIBILITY_TOLERANCE || supportUpper < 1 - FEASIBILITY_TOLERANCE) return undefined;
if (risk === 0) return continuousLinearSupport(mean, sigma, risk, lower, upper, support);
let best: BestPoint | undefined;
const dimension = support.length;
const stateCount = 3 ** dimension;
for (let code = 0; code < stateCount; code += 1) {
let stateCode = code;
const fixed = Array.from({ length: mean.length }, () => 0);
const fixedState = new Map<number, number>();
const free: number[] = [];
support.forEach((index) => {
const state = stateCode % 3; stateCode = Math.floor(stateCode / 3);
if (state === 1) free.push(index);
else { fixed[index] = state === 0 ? lower[index] : upper[index]; fixedState.set(index, state); }
});
const fixedTotal = stableSum(fixed);
if (free.length === 0) {
if (Math.abs(fixedTotal - 1) <= 10 * FEASIBILITY_TOLERANCE) best = pointIsBetter(fixed, best, mean, sigma, risk);
continue;
}
const dimensionFree = free.length;
const hessian = Array.from({ length: dimensionFree }, () => Array.from({ length: dimensionFree }, () => 0));
const linear = Array.from({ length: dimensionFree }, () => 0);
free.forEach((index, row) => {
free.forEach((other, column) => { hessian[row][column] = 2 * risk * sigma[index][other]; });
linear[row] = mean[index] - 2 * risk * dot(sigma[index], fixed);
});
const solution = solveKktSystem(hessian, linear, 1 - fixedTotal);
if (!solution) continue;
const point = [...fixed];
free.forEach((index, row) => { point[index] = solution[row]; });
if (Math.abs(stableSum(point) - 1) > 10 * FEASIBILITY_TOLERANCE) continue;
if (support.some((index) => point[index] < lower[index] - BOUND_TOLERANCE || point[index] > upper[index] + BOUND_TOLERANCE)) continue;
support.forEach((index) => { point[index] = Math.min(upper[index], Math.max(lower[index], point[index])); });
const multiplier = solution[dimensionFree];
const gradients = support.map((index) => mean[index] - 2 * risk * dot(sigma[index], point));
const stationarityScale = Math.max(
...gradients.map((value) => Math.abs(value)),
Math.abs(multiplier),
...support.map((index) => Math.abs(mean[index]) + 2 * risk * stableSum(sigma[index].map((value, column) => Math.abs(value) * Math.abs(point[column])))),
);
const kktTolerance = KKT_TOLERANCE * stationarityScale;
const valid = support.every((index, supportPosition) => {
const gradient = gradients[supportPosition];
const state = fixedState.get(index);
const width = upper[index] - lower[index];
if (state === 0 && width > BOUND_TOLERANCE) return gradient <= multiplier + kktTolerance;
if (state === 2 && width > BOUND_TOLERANCE) return gradient >= multiplier - kktTolerance;
return state === undefined ? Math.abs(gradient - multiplier) <= kktTolerance : true;
});
if (valid) best = pointIsBetter(point, best, mean, sigma, risk);
}
return best;
}
function supportTotal(size: number, cardinality: number): number {
let total = 0;
for (let count = 1; count <= cardinality; count += 1) {
let choose = 1;
for (let factor = 1; factor <= count; factor += 1) choose = choose * (size - count + factor) / factor;
total += Math.round(choose);
}
return total;
}
function supports(size: number, cardinality: number): number[][] {
const result: number[][] = [];
const visit = (start: number, target: number, prefix: number[]) => {
if (prefix.length === target) { result.push([...prefix]); return; }
for (let index = start; index < size; index += 1) visit(index + 1, target, [...prefix, index]);
};
for (let count = 1; count <= cardinality; count += 1) visit(0, count, []);
return result;
}
export function optimizeCardinalityConstrained(
mu: unknown, covariance: unknown, lambdaRisk: unknown, current: unknown, K: unknown,
options: { lower?: unknown; upper?: unknown; maxSupports?: unknown } = {},
): CardinalityResult {
const checked = validateInputs(mu, covariance, lambdaRisk, current, K, options.lower, options.upper);
const allSupports = supports(checked.mean.length, checked.cardinality);
const totalSupports = allSupports.length;
let supportLimit = totalSupports;
if (options.maxSupports !== undefined) {
if (typeof options.maxSupports !== "number" || !Number.isInteger(options.maxSupports) || options.maxSupports <= 0) invalid("invalid-input", "max_supports must be a positive integer.");
supportLimit = Math.min(options.maxSupports, totalSupports);
}
let best: BestPoint | undefined;
let enumerated = 0;
for (const support of allSupports) {
if (enumerated >= supportLimit) break;
enumerated += 1;
const candidate = continuousSupport(checked.mean, checked.sigma, checked.risk, checked.lower, checked.upper, support);
const comparisonTolerance = best === undefined ? 0 : RELATIVE_TOLERANCE * Math.max(Math.abs(candidate?.value ?? 0), Math.abs(best.value));
if (candidate && (best === undefined || candidate.value > best.value + comparisonTolerance)) best = candidate;
}
const complete = enumerated === totalSupports;
if (best === undefined) {
return {
status: complete ? "infeasible" : "incomplete", weights: null, selected: null, trade: null,
objective: null, bestBound: null, relativeGap: null, enumeratedSupports: enumerated,
totalSupports, certificate: {
kind: "exhaustive-support-enumeration", complete,
globalOptimality: complete ? "all supports were checked; no continuous support was feasible" : "no incumbent was found before the support limit",
},
};
}
const weights = best.point;
const selected = weights.map((value) => value > 0 ? 1 : 0);
const trade = weights.map((value, index) => value - checked.current[index]);
const realizedSupport = selected.reduce((total: number, value: number) => total + value, 0);
const tradeSupport = trade.filter((value) => Math.abs(value) > SUPPORT_TOLERANCE).length;
return {
status: complete ? "optimal" : "feasible", weights, selected, trade,
objective: best.value, bestBound: complete ? best.value : null, relativeGap: complete ? 0 : null,
enumeratedSupports: enumerated, totalSupports, realizedSupport, tradeSupport,
certificate: {
kind: "exhaustive-support-enumeration", complete, continuousKkt: true,
globalOptimality: complete ? "all supports and continuous KKT active sets were checked" : "incumbent only; support enumeration stopped before a proof",
},
};
}
Choose a scenario, step through the decisions and compare the result with the worked example.