Tensor trains vs. Monte Carlo

The price of a European option on a basket of $d$ assets is an expectation over $d$ correlated random variables. Monte Carlo is indifferent to $d$ but converges as $N^{-1/2}$, so each additional digit costs a hundred times the samples. This note instead approximates the integrand by a tensor train, a product of small matrices with one factor per variable, and integrates the approximation exactly¹. We compare it with Monte Carlo and randomized quasi-Monte Carlo on baskets of 5, 10 and 20 assets, charging the train for every evaluation spent building it.

At five assets the train reaches a relative error of $10^{-6}$ with about a hundredth of the evaluations quasi-Monte Carlo needs; at ten and twenty, quasi-Monte Carlo is cheaper at every accuracy both reached. The rank of the train, which grows with the number of assets under every correlation structure we tried, explains both.

The integral

Let $S_1,\dots,S_d$ follow correlated geometric Brownian motions under the risk-neutral measure $\mathbb Q$,

with risk-free rate $r$, maturity $T$, volatilities $\sigma_i$ and correlation matrix $R$. With weights $w_i$ and $x^+=\max(x,0)$, a call on the basket $\sum_i w_iS_i(T)$ with strike $K$ is worth

Writing $Y = Az$ for any $A$ with $AA^{\top}=R$, where $z$ is standard normal in $\mathbb R^d$, gives

where $\varphi_d(z)=\prod_k\varphi(z_k)$ is the standard normal density and $f$ the discounted payoff as a function of $z$. The choice of $A$ leaves the price unchanged but orders the variables, which matters to quasi-Monte Carlo and to the train. We take $A = V\Lambda^{1/2}$ from $R = V\Lambda V^{\top}$, eigenvalues decreasing, so that $z_1$ is the first principal component, for an equity basket the market factor.

The test problem is an equally weighted at-the-money call with $S_i(0)=K=100$, $T=1$ and $r=3\%$, volatilities evenly spaced from 15% to 40%, and correlation 0.6 within each of three sectors and 0.3 across, slightly perturbed. The reference prices are 9.04975, 8.72947 and 8.55218 for $d=5$, 10 and 20, each with a relative standard error below $3\times10^{-7}$.

A quadrature rule replaces $\int g(x)\,\varphi(x)\,dx$ by $\sum_i q_i\,g(x_i)$, with $n$ nodes $x_i$ and weights $q_i$. Applied along every axis it approximates $V_0$ by a sum over $n^d$ grid points,

Monte Carlo

Monte Carlo replaces the grid with independent draws $Z_n\sim\mathcal N(0,I_d)$,

The dimension enters only through $\sigma_f$, the standard deviation of the payoff, 1.49 times the price for ten assets. Figure i runs the estimator seven times on independent seeds, and their spread at each $N$ shows that standard deviation.

At $\htmlClass{vbind vbind-mc-n}{N}=$ 262,144 the estimate is 8.701 with a 95% half-width of 0.050, and 6 of the seven runs lie inside the envelope. A relative standard error of $10^{-6}$ would take $2.2\times10^{12}$ samples, two days on an Apple M4.

Randomized quasi-Monte Carlo draws from a scrambled Sobol sequence, which fills the unit cube more evenly, and estimates its error from independent scramblings, each unbiased. It does best on smooth integrands of low effective dimension; on the smoothed integrands used later its error falls roughly as $1/N$, which makes it the stronger baseline.

Low rank in two dimensions

Take first the two-asset integrand, payoff times correlated normal density, on a grid as a $64\times64$ matrix $F$. Its singular value decomposition

with singular values $\sigma_1\ge\sigma_2\ge\dots\ge0$ and orthonormal vectors $u_\alpha$, $v_\alpha$, gives the best rank-$r$ approximation in the Frobenius norm when cut off after $r$ terms. The singular values fall quickly and then slowly (figure ii), because of the kink where the basket crosses the strike.

With $\htmlClass{vbind vbind-svd-r}{r}=$ 3 terms the approximation stores 387 numbers instead of 4,096, with relative Frobenius error 4.1×10⁻². Low rank also makes weighted sums cheap:

so a two-dimensional quadrature becomes $r$ pairs of one-dimensional ones. With the density inside $F$ the weights are equal; the two sides differ by a relative 5.0×10⁻⁴, an error that need not fall monotonically with $r$.

The tensor-train format

A tensor train approximates a function on a product grid by a product of matrices, one per variable²,

Each core $G_k$ is an $r_{k-1}\times n\times r_k$ array holding a matrix for each of the $n$ grid values of $x_k$. Evaluating $f$ at a grid point multiplies one matrix from each core, passing a row vector of length $r_k$ from core $k$ to core $k+1$ until a number remains, at a cost of $O(dr^2)$ with $r=\max_k r_k$. The train stores $O(dnr^2)$ numbers instead of $n^d$, and the $r_k$ are its TT ranks.

The ten-asset train below has nine cores, for $z_2,\dots,z_{10}$, because $z_1$ is integrated out analytically (see the section on the kink). With ten grid points per core and ranks at most eight, it stores 4,640 numbers in place of $10^9$.

At the point of figure iii the train returns 8.481 after 464 multiply–adds, against 8.481 for the smoothed integrand itself, a relative error of 5.3×10⁻⁵. Drag any one lane and the error stays below 0.3%; drag several to the outermost nodes and the train breaks down, even to negative values.

What the rank measures

The ranks a train needs depend on $f$ and on the accuracy required. Arrange the grid values of $f$ into the matrix $F^{\langle k\rangle}$, the $k$-th unfolding, with rows indexed by $(x_1,\dots,x_k)$ and columns by $(x_{k+1},\dots,x_d)$. Then $f$ has an exact tensor train with $r_k=\operatorname{rank}F^{\langle k\rangle}$ and none with smaller ranks²: $r_k$ is the number of values the first $k$ variables must pass on for the rest to finish evaluating $f$.

A sum $g_1(x_1)+\dots+g_d(x_d)$ has rank 2 at every cut, since the left half need pass on only its running total and the constant 1:

The program below checks this for a basket of independent assets, $g_k=w_kS_k$, then computes the rank of one unfolding for the basket value $s$ and for the call payoff $(s-100)^+$.

// A basket of d independent assets S_k = 100 e^{0.2 x_k − 0.02}, held with equal
// weights w = 1/d, is worth s = g_1(x_1) + … + g_d(x_d) with g_k = w S_k. At the
// point x_k = cos k we multiply out the cores of the display, a 1 × 2 matrix,
// d − 2 matrices of size 2 × 2 and a 2 × 1, and compare the product with the sum.
const strike = 100;
const asset = (x: number): number => 100 * Math.exp(0.2 * x - 0.02);

type Matrix = number[][];
const times = (a: Matrix, b: Matrix): Matrix =>
	a.map((row) => (b[0] ?? []).map((_, j) => row.reduce((s, v, i) => s + v * (b[i]?.[j] ?? 0), 0)));

const cores = (g: readonly number[]): Matrix[] =>
	g.map((gk, k) => {
		if (k === 0) return [[gk, 1]];
		if (k === g.length - 1) return [[1], [gk]];
		return [
			[1, 0],
			[gk, 1],
		];
	});

for (const d of [5, 10, 20]) {
	const g = Array.from({ length: d }, (_, k) => asset(Math.cos(k + 1)) / d);
	const train = cores(g);
	const rank = Math.max(...train.map((core) => core[0]?.length ?? 0));
	const product = train.reduce(times)[0]?.[0] ?? 0;
	const sum = g.reduce((total, gk) => total + gk, 0);
	const label = `d = ${d}`.padEnd(9);
	console.log(`${label}rank ${rank}   product ${product.toFixed(8)}   sum ${sum.toFixed(8)}`);
}

// The rank of a matrix by elimination with full pivoting: subtract the rank-one
// cross through the largest remaining entry until no entry is left above the
// rounding level.
function rankOf(matrix: Matrix): number {
	const a = matrix.map((row) => [...row]);
	const tolerance = 1e-10 * Math.max(...a.flat().map(Math.abs));
	for (let r = 0; ; r++) {
		let [pi, pj, pivot] = [0, 0, 0];
		for (const [i, row] of a.entries())
			for (const [j, v] of row.entries())
				if (Math.abs(v) > Math.abs(pivot)) [pi, pj, pivot] = [i, j, v];
		if (Math.abs(pivot) <= tolerance) return r;
		const pivotRow = [...(a[pi] ?? [])];
		for (const row of a) {
			const f = (row[pj] ?? 0) / pivot;
			for (const [j, v] of row.entries()) row[j] = v - f * (pivotRow[j] ?? 0);
		}
	}
}

// The unfolding at the middle cut for d = 4 on 8 points per variable, rows
// indexed by (x_1, x_2) and columns by (x_3, x_4). Each half enters only through
// its partial sum. The call payoff (s − 100)⁺ is not affine in the left partial
// sum, so the pair (s, 1) can no longer carry it across the cut.
const terms = Array.from({ length: 8 }, (_, i) => asset(-2 + (4 * i) / 7) / 4);
const halves = terms.flatMap((a) => terms.map((b) => a + b));
const unfolding = (f: (s: number) => number): Matrix =>
	halves.map((left) => halves.map((right) => f(left + right)));
const linear = rankOf(unfolding((s) => s));
const call = rankOf(unfolding((s) => Math.max(s - strike, 0)));
console.log("middle cut of d = 4 on 8 points per variable");
console.log(`rank ${linear} for s, rank ${call} for (s − 100)⁺`);

A basket payoff depends on the first $k$ assets only through the partial sum $s_k=\sum_{i\le k}w_iS_i$. For independent assets, rounding $s_k$ to $m$ levels makes evaluating $f$ a dynamic program with rank at most $m$ at every cut. In the program each half of the grid has 36 distinct partial sums, and the call needs 30 of the 36 ranks they allow. Without the kink a few levels suffice; with correlation every asset depends on several coordinates, and more must cross each cut.

Figure iv draws, at each cut of the smoothed ten-asset integrand, the relative error of the best rank-$r$ approximation of $F^{\langle \htmlClass{vbind vbind-cut-k}{k}\rangle}$, $\htmlClass{vbind vbind-cut-rank}{\big(\sum_{\alpha>r}\sigma_\alpha^2\big)^{1/2}\big/\big(\sum_{\alpha}\sigma_\alpha^2\big)^{1/2}}$.

At cut 4 the unfolding needs rank 14 for a relative Frobenius error of $10^{-4}$; the ranks peak in the middle, where the cut divides the variables most evenly. This criterion is stricter than pricing requires.

Integration by contraction

On each axis we use the $n$-point Gauss–Hermite rule, for which $\sum_i q_i\,p(x_i)=\int p(x)\,\varphi(x)\,dx$ for every polynomial $p$ of degree below $2n$. For a tensor train the product-rule sum factorizes, and computing the right-hand side is called contracting the train:

Each factor costs $O(nr^2)$ and each product from the left $O(r^2)$, so the integral costs $O(dnr^2)$ instead of $O(n^d)$. The program below sums the $4^6$ weighted entries of a random train with $d=6$, $n=4$ and ranks 3 one at a time, then contracts it, counting multiply–adds.

// A random tensor train with d variables, n values each and every rank r.
// Core k stores one matrix per value of x_k, row-major. Raise d and run again:
// the brute-force sum grows as n^d, the contraction as d·n·r².
const d = 6;
const n = 4;
const r = 3;
const weights = [0.1, 0.4, 0.4, 0.1]; // quadrature weights q_i, one per value of x_k
const ranks = Array.from({ length: d + 1 }, (_, k) => (k === 0 || k === d ? 1 : r));

let seed = 1;
const random = (): number => {
	seed = (seed * 48271) % 2147483647;
	return seed / 2147483647;
};

const cores = Array.from({ length: d }, (_, k) =>
	Array.from({ length: n }, () =>
		Float64Array.from({ length: (ranks[k] ?? 1) * (ranks[k + 1] ?? 1) }, random),
	),
);

// v ← v · M for a row vector v and a matrix M with `cols` columns.
let multiplyAdds = 0;
function times(v: Float64Array, m: Float64Array, cols: number): Float64Array {
	const out = new Float64Array(cols);
	v.forEach((vi, i) => {
		for (let j = 0; j < cols; j++) out[j] = (out[j] ?? 0) + vi * (m[i * cols + j] ?? 0);
	});
	multiplyAdds += v.length * cols;
	return out;
}

// One entry: each core contributes the matrix its index selects.
function entry(index: readonly number[]): number {
	let v: Float64Array = Float64Array.of(1);
	cores.forEach((core, k) => {
		v = times(v, core[index[k] ?? 0] ?? new Float64Array(), ranks[k + 1] ?? 1);
	});
	return v[0] ?? 0;
}

// Brute force: every one of the n^d entries, weighted. Each term costs the d
// vector–matrix products inside entry(), which times() counts, and d more:
// d − 1 multiplications for its weight and one multiply–add into the sum.
let brute = 0;
for (let flat = 0; flat < n ** d; flat++) {
	const index = Array.from({ length: d }, (_, k) => Math.floor(flat / n ** k) % n);
	brute += index.reduce((q, i) => q * (weights[i] ?? 0), 1) * entry(index);
	multiplyAdds += d;
}
const bruteCost = multiplyAdds;

// Contraction: collapse each core to Σ_i q_i G_k(i), then run one row vector
// down the train.
multiplyAdds = 0;
let carried: Float64Array = Float64Array.of(1);
cores.forEach((core, k) => {
	const summed = new Float64Array(core[0]?.length ?? 0);
	core.forEach((m, i) => {
		m.forEach((x, e) => {
			summed[e] = (summed[e] ?? 0) + (weights[i] ?? 0) * x;
		});
	});
	multiplyAdds += n * summed.length;
	carried = times(carried, summed, ranks[k + 1] ?? 1);
});
const integral = carried[0] ?? 0;

const count = (x: number): string => x.toLocaleString("en-US");
const report = (label: string, value: number, cost: number): void =>
	console.log(`${label}   ${value.toPrecision(15)}   ${count(cost)} multiply–adds`);
console.log(`entries       ${count(n ** d)}`);
report("brute force", brute, bruteCost);
report("contraction", integral, multiplyAdds);

Contracting the ten-asset train takes 5,104 multiply–adds. Its relative error of $2.2\times10^{-5}$ comes from the ten-point grid and the rank-8 approximation; the contraction adds only rounding error.

Construction by cross approximation

The textbook construction, TT-SVD, truncates SVDs of the unfoldings and so needs all $n^d$ values; cross approximation does not. For a matrix $A$ of rank $r$, any $r$ rows $I$ and $r$ columns $J$ whose intersection $A[I,J]$ is nonsingular determine $A$:

If $A$ is only approximately of rank $r$, an intersection of maximal volume³, the largest $|\det|$, bounds the entrywise error by $(r+1)\,\sigma_{r+1}(A)$, where $\sigma_{r+1}(A)$ is the spectral-norm error of the best rank-$r$ approximation. Figure v builds a cross greedily: each step requests a row, pivots on its largest residual and requests the column through the pivot, whose largest residual among unused rows picks the next row. No maximal-volume bound covers this partial pivoting.

After 3 steps it has evaluated 375 of the 4,096 entries for a relative error of 9.5×10⁻², against 4.1×10⁻² for the truncated SVD of the same rank, which needs every entry. Beyond five steps the greedy approximation stalls.

TT-cross applies the idea at every cut⁴, keeping $r_k$ settings of $(x_1,\dots,x_k)$ and $r_k$ of $(x_{k+1},\dots,x_d)$. At core $k$ it evaluates $f$ at every combination of a kept setting on the left, a grid value of $x_k$ and a kept setting on the right, and picks the settings for the next cut with maxvol⁵, which finds a submatrix of locally maximal volume by row exchanges. A sweep along the train and back visits every core twice, so it costs twice the numbers the train stores. Exactly maximal volume at every step would keep the entrywise error within a factor of the best at the same ranks, a factor depending on $r$ and $d$ but not exponentially on $d$⁶; maxvol finds only a local maximum, and on a black box cross approximation can miss features it never samples.

Each entry is first multiplied by $\prod_k\sqrt{q_{i_k}}$, so pivoting on large entries means pivoting on those that carry the integral, and contracting against $\sqrt q$ recovers it. As the weights on each axis sum to one, the Cauchy–Schwarz inequality bounds the error the train adds to the quadrature sum by the Frobenius error of the scaled train. The TT-cross listed in the appendix reaches a relative error of $6\times10^{-9}$ on the five-asset basket from 19,968 evaluations, 30% of its $16^4$-point grid.

The kink

The raw payoff $(B(z)-K)^+$, with $B(z)=\sum_i w_iS_i(z)$, is continuous but not differentiable on the exercise boundary $B(z)=K$. Across this kink polynomial quadrature converges only as a power of the number of nodes instead of exponentially, and the singular values of the unfoldings decay slowly. Switched to raw, figure iv adds $z_1$ and a ninth cut, and every cut shared with the smoothed integrand needs a higher rank, the largest rising from 14 to 23; raw ranks of 22 and 23, just under the rank-24 cap of the plate’s trains, are lower bounds.

The remedy is to integrate the kink out analytically, a device known as pre-integration⁷ or conditional smoothing⁸. Write $w_iS_i(T)=c_i\,e^{b_iz_1}$, where $b_i=\sigma_i\sqrt T\,A_{i1}$ is the loading of asset $i$ on the market factor and $c_i$ depends on $z_{2:d}=(z_2,\dots,z_d)$. Every loading is positive, so for fixed $z_{2:d}$ the basket increases in $z_1$ and crosses $K$ at a single point $z_1^\star(z_{2:d})$, found by Newton’s method. Each term is lognormal in $z_1$, and with $\Phi$ the standard normal distribution function, $\int_t^\infty e^{bz}\varphi(z)\,dz=e^{b^2/2}\,\Phi(b-t)$ gives a formula of Black–Scholes type, analytic in the remaining $d-1$ variables and with the same integral:

In figure vi, at $\htmlClass{vbind vbind-smooth-z2}{z_2}=$ 0.60, the boundary sits at $\htmlClass{vbind vbind-smooth-boundary}{z_1^\star}=$ 0.07 and the conditional expectation is 7.29. On five assets the best raw train reached a relative error of $1.4\times10^{-3}$ from $1.9\times10^{5}$ evaluations, and TT-cross on the smoothed integrand $6\times10^{-9}$ from a tenth as many. Monte Carlo and quasi-Monte Carlo get the same smoothed integrand, which for ten assets cuts the Monte Carlo standard deviation from 1.49 to 0.17 times the price, saving a factor of about 80 in samples and 16 in time.

Accuracy versus cost

Figure vii compares the errors of Monte Carlo (MC), randomized quasi-Monte Carlo (RQMC) and tensor trains (TT) on the three baskets, raw and smoothed, counting as cost every evaluation of the integrand, repeats included. Train errors are not monotone in rank, so the cheapest train to reach each error is picked with the reference in hand; in practice a second train or the control variate below would have to confirm it, at a cost left out here.

For $\htmlClass{vbind vbind-cost-d}{d}=$ 5, a relative error of $10^{-4}$ costs the tensor train 1,824 evaluations against 3,700 for smoothed quasi-Monte Carlo, and $10^{-6}$ costs it 4,840 against 490,000. Monte Carlo on the smoothed integrand would need 6.9×10¹⁰ for $10^{-6}$.

At five assets the train loses only at low accuracy, with seven times the evaluations of quasi-Monte Carlo at $10^{-3}$. At ten it needs 55 times fewer than smoothed Monte Carlo at $10^{-4}$, but quasi-Monte Carlo is cheaper than our trains at every accuracy both reached, thirteen times at $10^{-4}$ and three at $10^{-6}$. At twenty the errors no longer fall with rank: from rank 2 to 40 they scatter between $6.5\times10^{-6}$ and $1.4\times10^{-3}$, and the rank-40 train, built from $7\times10^{6}$ evaluations, is less accurate than the rank-2 train built from 3,456. Quasi-Monte Carlo is cheaper than every train we built.

An independent implementation, teneva⁹, slows the same way: at rank 16 and $d=20$ it needs 2.8 million evaluations for $9.4\times10^{-5}$, where quasi-Monte Carlo needs 4,096 for $7.8\times10^{-5}$. Wall-clock times, measured for every method in one Python process, give the same ordering at errors of $10^{-5}$ and below.

Why the rank grows

The eigenvalues of the correlation matrix explain the twenty-asset result (figure viii). For $\htmlClass{vbind vbind-growth-d}{d}=$ 20, the market and two sector directions carry 65% of the variance, and the remaining 17 directions are of similar size. Each moves every asset, so each cut must carry information about many assets. The largest rank any cut needs for $10^{-4}$ grows with $d$ under all three correlation structures we tried: the sector blocks, a single factor with every pairwise correlation 0.5, and random positive loadings on $d/2$ factors.

The same spectrum suits quasi-Monte Carlo, whose effective dimension the three dominant directions keep small: its cost to reach $10^{-6}$ grew by less than a factor of two from ten to twenty assets. Known rank bounds cover Gaussian densities in correlated coordinates¹⁰; in principal-component coordinates the density has rank 1, so the ranks here come from the payoff, for which we know of no bound.

A tensor train therefore needs an integrand that is smooth, or can be made so, and dominated by a few interactions, with strongly coupled variables adjacent, and its construction must be repaid by high accuracy or by many integrals sharing one train. A single ten-asset price to $10^{-3}$ took smoothed quasi-Monte Carlo 256 evaluations, our smallest sample, and the cheapest train 4,608.

Tensor trains as control variates

A contraction returns a number with no error estimate, and TT-cross supplies none. The usual rule stops the iteration once the integral changes little from one sweep to the next; ours ran a fixed number of sweeps, and figure ix replays the twelve that built the ten-asset train.

After sweep 12 and 110,800 evaluations the train integrates to 8.72928, a relative error of 2.2×10⁻⁵, and its relative change from the previous sweep is 1.3×10⁻⁴. The ranks settle at sweep 3; each later sweep costs 9,280 evaluations, and the error moves between $1.3\times10^{-6}$ and $2.5\times10^{-4}$. Drag the sweep and compare the dashed step from the previous sweep, all a stopping rule sees, with the solid drop to the reference: at sweep 6 the step is $3.7\times10^{-6}$ and the error $1.0\times10^{-4}$.

As a control variate the train gets an error estimate from Monte Carlo, which then estimates only what the train got wrong. Let $f$ now be the smoothed integrand, $f_{TT}$ the train extended off the grid by polynomial interpolation through the Gauss–Hermite nodes, and $I_{TT}$ its integral, which the contraction computes exactly because the rule integrates that polynomial without error. Then

is unbiased for any train, however crude. Its variance is $\htmlClass{vbind vbind-cv-ratio}{\operatorname{Var}(f-f_{TT})}/N$, small when the train is accurate, and the central limit theorem gives a confidence interval as for plain Monte Carlo.

After 32,768 samples the 95% intervals of figure x are 8.72948 ± 7.8×10⁻⁵ with the train and 8.71382 ± 0.01548 without it, against the reference 8.72947. At the estimated variance ratio of 39,408 the band with the train is too narrow to see. Over five trains of ranks 2 to 16, each the control variate for 4,096 samples, the variance fell by factors from 24 to $1.2\times10^6$, and all five 95% intervals contained the reference.

The program below shows the mechanism in one dimension, with Taylor polynomials of rising degree in place of the train. Every interval covers the exact value. The degree-8 control divides the variance by 4,190, and the degree-1 control increases it.

// E[f(Z)] for f(z) = e^{0.3z} cos z and Z ~ N(0, 1); the exact value is
// e^{(0.09 − 1)/2} cos 0.3. Each control p is a Taylor polynomial of f at 0.
// Its mean under N(0, 1) follows exactly from the moments of Z, and it enters
// with coefficient one, as the train does. Change the degrees or the number of
// draws and run it again.
const degrees = [1, 2, 4, 8];
const draws = 20000;
const a = 0.3;
const f = (z: number): number => Math.exp(a * z) * Math.cos(z);
const exact = Math.exp((a * a - 1) / 2) * Math.cos(a);

// Taylor coefficients of e^{az} cos z = Re e^{(a+i)z}: c_k = Re (a+i)^k / k!.
const coefficients: number[] = [];
let re = 1;
let im = 0;
let factorial = 1;
for (let k = 0; k <= Math.max(...degrees); k++) {
	coefficients.push(re / factorial);
	[re, im] = [re * a - im, re + im * a];
	factorial *= k + 1;
}

// E[Z^k] = (k − 1)!! for even k, 0 for odd k.
function moment(k: number): number {
	if (k % 2 === 1) return 0;
	return k === 0 ? 1 : (k - 1) * moment(k - 2);
}

let seed = 7;
const uniform = (): number => {
	seed = (seed * 48271) % 2147483647;
	return seed / 2147483647;
};
const normal = (): number =>
	Math.sqrt(-2 * Math.log(uniform())) * Math.cos(2 * Math.PI * uniform());
const sample = Array.from({ length: draws }, normal);

// Sample mean, 95% half-width and variance, the variance computed in two passes.
function stats(values: readonly number[]): { mean: number; half: number; variance: number } {
	const mean = values.reduce((s, v) => s + v, 0) / values.length;
	const variance = values.reduce((s, v) => s + (v - mean) ** 2, 0) / (values.length - 1);
	return { mean, half: 1.96 * Math.sqrt(variance / values.length), variance };
}

// Every estimator uses the same draws. With a control, the estimate is the exact
// mean of p plus the sample mean of f − p.
const row = (label: string, mean: number, half: number): string =>
	`${label.padEnd(17)}${mean.toFixed(6)} ± ${half.toFixed(6)}`;
const plain = stats(sample.map(f));
console.log(`exact            ${exact.toFixed(6)}`);
console.log(row("plain MC", plain.mean, plain.half));
for (const degree of degrees) {
	const c = coefficients.slice(0, degree + 1);
	const p = (z: number): number => c.reduceRight((acc, ck) => acc * z + ck, 0);
	const meanOfP = c.reduce((s, ck, k) => s + ck * moment(k), 0);
	const { mean, half, variance } = stats(sample.map((z) => f(z) - p(z)));
	const ratio = plain.variance / variance;
	const shown = ratio < 10 ? ratio.toPrecision(2) : Math.round(ratio).toLocaleString("en-US");
	const label = `degree ${degree}`;
	console.log(`${row(label, meanOfP + mean, half)}   variance ÷ ${shown}`);
}
console.log(`± is a 95% half-width from ${draws.toLocaleString("en-US")} draws`);

A train used as the proposal in importance sampling or Metropolis–Hastings likewise costs efficiency when inaccurate and leaves the estimator correct. Such a train can be sampled exactly¹¹ and kept nonnegative¹², and a defensive mixture with the original density bounds the importance weights¹³ and makes an independence chain converge geometrically from any start¹⁴.

Related work

Two financial problems suited to these uses remain untested, credit-portfolio tail probabilities and default contagion²¹.

Reproducing

The experiments and a Python tensor-train toolkit are at github.com/conorscheidt/tensor-trains. In about ten minutes on an Apple M4, run_pricing.py regenerates every measured number here except the teneva comparison, which comes from crosscheck_teneva.py.

Appendix: a complete TT-cross

The program below is a one-site TT-cross, updating one core at a time, for the five-asset basket. It approximates the smoothed integrand, scaled by $\sqrt q$, on a 16-point Gauss–Hermite grid in each of $z_2,\dots,z_5$ and checks the result against the dense sum over all 65,536 points.

// One-site TT-cross for a call on an equally weighted basket of five assets,
// checked and timed against the dense quadrature sum on the same grid. The
// basket is the five-asset test problem of the post: S_i(0) = K = 100, T = 1,
// r = 3%, volatilities evenly spaced from 15% to 40%, and correlation 0.6 within
// each of three sectors and 0.3 across, slightly perturbed. In principal-component
// coordinates z_1, …, z_5, with z_1 the market factor, the payoff is integrated
// over z_1 in closed form, which leaves a smooth integrand in z_2, …, z_5. Both
// methods use the n-point Gauss–Hermite rule on each of these four axes. The
// dense sum evaluates the integrand at all n^4 grid points. TT-cross evaluates
// it only on the fibers its sweeps select and integrates the resulting train by
// contraction; once n is at least the rank, its evaluation count grows linearly
// in n. Try n = 8 and n = 32.
const rank = 12; // the rank at each of the three interior cuts, capped by the grid
const n = 16; // Gauss–Hermite points per variable
const sweeps = 2; // each sweep runs left to right, then right to left

// log S_T = mu + B z, taken from Basket(5) in basket.py of the experiment
// repository; each asset carries weight w in the basket. The reference price is
// the one the post reports for five assets, to ten significant digits.
// biome-ignore format: a table
const mu = [4.623920185988092, 4.612592060988092, 4.597357685988092, 4.578217060988091, 4.555170185988092];
// biome-ignore format: a table
const B = [
	[0.10682750963810546, 0.07906726507684363, 0.01325677104506904, 0.045081827616054324, 0.051265364702141866],
	[0.15339005657210322, -0.10539738439838546, 0.023223121065429723, 0.08096189595219622, -0.05852341376399753],
	[0.1678260686289569, -0.017291345820175212, -0.21520247925860803, -0.028885788159002907, 0.0037315159998804665],
	[0.2364057394315545, 0.17819491404058216, 0.06611260640778252, -0.09299427752734549, -0.11509273041780789],
	[0.27643173655071274, -0.20191106730777883, 0.11491252243985428, -0.13425588510204967, 0.10764704096335487],
];
const w = 0.2;
const strike = 100;
const discount = Math.exp(-0.03);
const reference = 9.049751584;
const d = mu.length - 1; // the train's variables are z_2, …, z_5

const at = (a: ArrayLike<number>, i: number): number => a[i] ?? 0;

// The n-point Gauss–Hermite rule for N(0, 1). Its nodes are the roots of the
// orthonormal Hermite polynomial p_n, all inside ±√(4n + 2), bracketed on a fine
// grid and bisected to machine precision. Its weights q_i = 1/(n p_{n−1}(x_i)²)
// sum to 1.
function hermite(x: number): [number, number] {
	let [previous, current] = [0, 1];
	for (let k = 0; k < n; k++)
		[previous, current] = [current, (x * current - Math.sqrt(k) * previous) / Math.sqrt(k + 1)];
	return [current, previous];
}
const negative = (x: number): boolean => hermite(x)[0] < 0;
const edge = Math.sqrt(4 * n + 2);
const nodes: number[] = [];
for (let i = 0; i * 1e-3 < 2 * edge; i++) {
	let lo = -edge + i * 1e-3;
	let hi = lo + 1e-3;
	if (negative(lo) === negative(hi)) continue;
	for (let mid = (lo + hi) / 2; mid !== lo && mid !== hi; mid = (lo + hi) / 2)
		if (negative(mid) === negative(lo)) lo = mid;
		else hi = mid;
	nodes.push(lo);
}
const weights = nodes.map((x) => 1 / (n * hermite(x)[1] ** 2));

// Φ(x) = ½ + φ(x)(x + x³/3 + x⁵/15 + …), Marsaglia's series, summed until the
// terms stop changing the total. Its absolute error is below 1e-15.
function cdf(x: number): number {
	if (Math.abs(x) > 8) return x > 0 ? 1 : 0;
	let sum = x;
	let term = x;
	for (let k = 3; sum + (term * x * x) / k !== sum; k += 2) {
		term *= (x * x) / k;
		sum += term;
	}
	return 0.5 + (sum * Math.exp(-0.5 * x * x)) / Math.sqrt(2 * Math.PI);
}

// E[payoff | z_2, …, z_5]: Newton finds the market level z* at which the basket
// hits the strike (it is increasing in z_1), then z_1 is integrated in closed form.
function smoothedPayoff(rest: readonly number[]): number {
	const b = B.map((loadings) => at(loadings, 0));
	const logc = mu.map(
		(m, i) => Math.log(w) + m + rest.reduce((s, z, j) => s + at(B[i] ?? [], j + 1) * z, 0),
	);
	let t = 0;
	for (let it = 0; it < 50; it++) {
		const terms = logc.map((lc, i) => Math.exp(lc + at(b, i) * t));
		const sum = terms.reduce((s, x) => s + x, 0);
		const step =
			(Math.log(sum) - Math.log(strike)) / (terms.reduce((s, x, i) => s + x * at(b, i), 0) / sum);
		t -= step;
		if (Math.abs(step) < 1e-13) break;
	}
	const value = logc.reduce(
		(s, lc, i) => s + Math.exp(lc + at(b, i) ** 2 / 2) * cdf(at(b, i) - t),
		0,
	);
	return discount * Math.max(value - strike * cdf(-t), 0);
}

type Matrix = { rows: number; cols: number; data: Float64Array };
const matrix = (rows: number, cols: number): Matrix => ({
	rows,
	cols,
	data: new Float64Array(rows * cols),
});
const get = (m: Matrix, i: number, j: number): number => at(m.data, i * m.cols + j);
const column = (m: Matrix, j: number): number[] =>
	Array.from({ length: m.rows }, (_, i) => get(m, i, j));

// An orthonormal basis U for the columns (modified Gram–Schmidt).
function orthonormal(m: Matrix): Matrix {
	const u = matrix(m.rows, m.cols);
	u.data.set(m.data);
	for (let j = 0; j < u.cols; j++) {
		for (let p = 0; p < j; p++) {
			const dot = column(u, p).reduce((s, x, i) => s + x * get(u, i, j), 0);
			for (let i = 0; i < u.rows; i++) u.data[i * u.cols + j] = get(u, i, j) - dot * get(u, i, p);
		}
		const norm = Math.hypot(...column(u, j)) || 1;
		for (let i = 0; i < u.rows; i++) u.data[i * u.cols + j] = get(u, i, j) / norm;
	}
	return u;
}

// X = U · U[rows]⁻¹, by Gauss–Jordan on the transposed system U[rows]ᵀ Xᵀ = Uᵀ.
function interpolate(u: Matrix, rows: readonly number[]): Matrix {
	const r = u.cols;
	const a = Array.from({ length: r }, (_, i) => rows.map((row) => get(u, row, i)));
	const x = Array.from({ length: r }, (_, i) => column(u, i));
	for (let c = 0; c < r; c++) {
		let p = c;
		for (let i = c + 1; i < r; i++) if (Math.abs(a[i]?.[c] ?? 0) > Math.abs(a[p]?.[c] ?? 0)) p = i;
		[a[c], a[p]] = [a[p] ?? [], a[c] ?? []];
		[x[c], x[p]] = [x[p] ?? [], x[c] ?? []];
		for (let i = 0; i < r; i++) {
			if (i === c) continue;
			const f = (a[i]?.[c] ?? 0) / (a[c]?.[c] ?? 1);
			a[i] = (a[i] ?? []).map((v, j) => v - f * (a[c]?.[j] ?? 0));
			x[i] = (x[i] ?? []).map((v, k) => v - f * (x[c]?.[k] ?? 0));
		}
	}
	const out = matrix(u.rows, r);
	for (let i = 0; i < r; i++)
		for (let k = 0; k < u.rows; k++) out.data[k * r + i] = (x[i]?.[k] ?? 0) / (a[i]?.[i] ?? 1);
	return out;
}

// maxvol: r rows of a tall U whose r × r block has locally maximal |det|. Start
// from the rows partial pivoting picks, then swap rows while some entry of
// U · U[rows]⁻¹ exceeds 1.05 in modulus.
function maxvol(u: Matrix): number[] {
	const work = Array.from({ length: u.rows }, (_, i) =>
		Array.from(u.data.subarray(i * u.cols, (i + 1) * u.cols)),
	);
	const rows: number[] = [];
	for (let j = 0; j < u.cols; j++) {
		let p = -1;
		for (const [i, row] of work.entries()) {
			if (!rows.includes(i) && (p < 0 || Math.abs(row[j] ?? 0) > Math.abs(work[p]?.[j] ?? 0)))
				p = i;
		}
		rows.push(p);
		const pivot = work[p] ?? [];
		for (const [i, row] of work.entries()) {
			const f = (row[j] ?? 0) / (pivot[j] ?? 1);
			if (i !== p) work[i] = row.map((v, c) => v - f * (pivot[c] ?? 0));
		}
	}
	for (let swap = 0; swap < 100; swap++) {
		const x = interpolate(u, rows);
		let best = { i: 0, j: 0, v: 0 };
		for (let i = 0; i < x.rows; i++) {
			for (let j = 0; j < x.cols; j++)
				if (Math.abs(get(x, i, j)) > best.v) best = { i, j, v: Math.abs(get(x, i, j)) };
		}
		if (best.v <= 1.05) break;
		rows[best.j] = best.i;
	}
	return rows;
}

const sqrtWeights = weights.map(Math.sqrt);

// Entries are scaled by Π √q, so pivoting on large entries pivots on the ones
// that carry the integral; the contraction below supplies the other √q.
// Sweeps request some grid points more than once, so both counts are kept.
let evaluations = 0;
const points = new Set<number>();
function entry(index: readonly number[]): number {
	evaluations++;
	points.add(index.reduce((flat, i) => flat * n + i, 0));
	return (
		smoothedPayoff(index.map((i) => at(nodes, i))) *
		index.reduce((p, i) => p * at(sqrtWeights, i), 1)
	);
}

// Left index sets I[k] (prefixes) start empty; right sets J[k] (suffixes) start
// random, as many as the rank allows: no more than n^k prefixes can sit on the left.
let seed = 12345;
const draw = (): number => {
	seed = (seed * 48271) % 2147483647;
	return seed % n;
};
function suffixes(length: number): number[][] {
	const seen = new Map<string, number[]>();
	while (seen.size < Math.min(rank, n ** length, n ** (d - length))) {
		const s = Array.from({ length }, draw);
		seen.set(s.join(), s);
	}
	return [...seen.values()];
}
const left: number[][][] = Array.from({ length: d + 1 }, () => [[]]);
const right: number[][][] = Array.from({ length: d + 1 }, (_, k) =>
	k === 0 || k === d ? [[]] : suffixes(d - k),
);
const cores: { r0: number; r1: number; data: Float64Array }[] = [];

// The block of entries with a prefix from I[k], any x_k, and a suffix from J[k+1],
// as an (r0 · n) × r1 matrix.
function fiber(k: number): Matrix {
	const pre = left[k] ?? [[]];
	const suf = right[k + 1] ?? [[]];
	const m = matrix(pre.length * n, suf.length);
	for (const [a, p] of pre.entries()) {
		for (let i = 0; i < n; i++)
			for (const [b, s] of suf.entries())
				m.data[(a * n + i) * suf.length + b] = entry([...p, i, ...s]);
	}
	return m;
}

const crossStart = performance.now();
for (let sweep = 0; sweep < sweeps; sweep++) {
	// left to right: maxvol on each fiber picks the prefixes I[k+1]
	for (let k = 0; k < d - 1; k++) {
		const u = orthonormal(fiber(k));
		const rows = maxvol(u);
		cores[k] = { r0: (left[k] ?? []).length, r1: u.cols, data: interpolate(u, rows).data };
		left[k + 1] = rows.map((s) => [...((left[k] ?? [])[Math.floor(s / n)] ?? []), s % n]);
	}
	cores[d - 1] = { r0: (left[d - 1] ?? []).length, r1: 1, data: fiber(d - 1).data };
	// right to left: the same step on the transposed fiber picks the suffixes J[k]
	for (let k = d - 1; k > 0; k--) {
		const f = fiber(k);
		const r0 = f.rows / n;
		const r1 = f.cols;
		const t = matrix(n * r1, r0);
		for (let a = 0; a < r0; a++)
			for (let i = 0; i < n; i++)
				for (let b = 0; b < r1; b++) t.data[(i * r1 + b) * r0 + a] = get(f, a * n + i, b);
		const u = orthonormal(t);
		const rows = maxvol(u);
		const x = interpolate(u, rows);
		const data = new Float64Array(u.cols * n * r1);
		for (let a = 0; a < u.cols; a++)
			for (let i = 0; i < n; i++)
				for (let b = 0; b < r1; b++) data[(a * n + i) * r1 + b] = get(x, i * r1 + b, a);
		cores[k] = { r0: u.cols, r1, data };
		right[k] = rows.map((s) => [Math.floor(s / r1), ...((right[k + 1] ?? [])[s % r1] ?? [])]);
	}
	cores[0] = { r0: 1, r1: (right[1] ?? []).length, data: fiber(0).data };
}

// Contract against √q, one core at a time.
let carried = Float64Array.of(1);
for (const { r0, r1, data } of cores) {
	const next = new Float64Array(r1);
	for (let a = 0; a < r0; a++)
		for (let i = 0; i < n; i++)
			for (let b = 0; b < r1; b++)
				next[b] =
					at(next, b) + at(carried, a) * at(sqrtWeights, i) * at(data, (a * n + i) * r1 + b);
	carried = next;
}
const price = at(carried, 0);
const crossMs = performance.now() - crossStart;

// The dense sum: every grid point, weighted by Π q.
const denseStart = performance.now();
let dense = 0;
for (let flat = 0; flat < n ** d; flat++) {
	const index = Array.from({ length: d }, (_, k) => Math.floor(flat / n ** k) % n);
	dense +=
		index.reduce((p, i) => p * at(weights, i), 1) * smoothedPayoff(index.map((i) => at(nodes, i)));
}
const denseMs = performance.now() - denseStart;

const count = (x: number): string => x.toLocaleString("en-US");
const error = (x: number): string => (Math.abs(x - reference) / reference).toExponential(1);
const line = (label: string, value: string, err = "", cost = "", ms = ""): void => {
	const numbers = err.padStart(16) + cost.padStart(14) + ms.padStart(7);
	console.log((label.padEnd(11) + value.padEnd(10) + numbers).trimEnd());
};
line("", "price", "relative error", "evaluations", "ms");
line("reference", reference.toFixed(8));
line("dense sum", dense.toFixed(8), error(dense), count(n ** d), denseMs.toFixed(0));
line("TT-cross", price.toFixed(8), error(price), count(evaluations), crossMs.toFixed(0));
console.log(`TT-cross evaluated ${count(points.size)} distinct grid points of ${count(n ** d)}`);

The dense sum agrees with the reference within the reference’s standard error, so the $5.6\times10^{-9}$ of TT-cross is compression error. TT-cross requested 19,968 entries, 11,275 of them distinct, in about the dense sum’s time. Its count grows linearly in $n$ once $n$ reaches the rank: at $n=32$ it needs 39,936 evaluations against 1,048,576 and runs over ten times faster, while at $n=8$ its 6,656 exceed the 4,096 grid points.