Quant Finance: Solving the High-Dimensional Black-Scholes PDE via Quantized Tensor Trains (QTT)
A complete Python guide to overcoming the curse of dimensionality in multi-asset option pricing using tensor decompositions.
Use the button at the end of this article to download the source code
This is not an official implementation. Since no official public code is available, I wrote this version based on my understanding of the approach. It has not been thoroughly tested, so please treat it as a reference for understanding the implementation rather than production-ready code.
Paper we are implementing today: https://arxiv.org/abs/2601.00009
This paper presents quantized tensor train (QTT) solvers for the multi-asset Black-Scholes PDE, overcoming the curse of dimensionality. Two solvers are developed: a time-stepping algorithm for European and American options, and a space-time algorithm for European options. Both provide full-grid solutions and Greeks with polynomial scaling in the number of assets and polylogarithmic scaling in grid size. Numerical experiments for basket and max-min options in 3-5 dimensions demonstrate high accuracy on a personal computer.
Implementation Assumptions
Fully implicit scheme (θ=1) for time-stepping.
TT-Cross sweeps default to 3-4.
Time-stepping solver is primary; space-time solver for d≤3.
Boundary conditions for basket put: lower faces as Eqs. (A20)-(A22), upper faces homogeneous Dirichlet (Eq. A23).
For max-min put: lower boundaries Ke^{-rt}, upper boundaries zero.
Grid dimensions: c=7-9 spatial cores per dimension; domain determined adaptively from pilot coarse run.
Market parameters from Appendix B.2.
ALS/MALS solver uses 2 sweeps by default.
Rank reduction after early-exercise uses TT-SVD with tolerance or target rank.
Introduction to Multi-Asset Option Pricing and the Curse of Dimensionality
Pricing a financial option on a single underlying asset is a textbook problem. The Black–Scholes model gives a closed-form solution for European calls and puts, and a simple one-dimensional finite-difference grid can handle American options or more complex payoffs with ease. But the moment you have two, three, or five underlying assets, the situation changes dramatically. The number of grid points needed to represent the solution grows exponentially with the number of assets, quickly exceeding the memory of any computer. This is the curse of dimensionality, and it is the central challenge that the paper by Kazeev, Khoromskij, and Tyrtyshnikov (2601.00009v2) overcomes using quantized tensor trains (QTT).
The Multi-Asset Black–Scholes PDE
The starting point is the d-dimensional Black–Scholes partial differential equation. For an option whose value depends on d asset prices S_1, S_2, ..., S_d and time t, the PDE is:
Here V(S,t) is the option price, a scalar function of the d asset prices and time. Each σ_i is the volatility of asset i, ρ_{ij} is the correlation between assets i and j, and r is the risk-free interest rate. The double sum over i and j couples every pair of assets through the covariance matrix σ_i σ_j ρ_{ij}. This coupling is what makes the problem genuinely high-dimensional: the price of a basket option depends on the joint distribution of all underlyings, not just their marginal distributions.
Why Classical Methods Fail
A standard approach to solving this PDE is to discretise the spatial domain with a uniform grid. If we use N grid points per asset, the total number of unknowns is N^d. For d=5 assets and a modest N=256 points per dimension, that is 256^5 ≈ 1.1 × 10^12 unknowns — over a trillion grid points. Storing a single vector of that size would require terabytes of memory, and solving the resulting linear system is completely infeasible on any personal computer.
Monte Carlo simulation avoids the grid entirely by sampling random paths, but it has its own limitations. It provides only point estimates of the option price at a given initial spot vector, not the full solution surface. Computing Greeks (Delta, Gamma) requires additional simulation or finite-difference approximations, which introduce noise and extra computational cost. Moreover, Monte Carlo converges slowly, at a rate of 1/√(number of paths), making high accuracy expensive.
The Promise of Quantized Tensor Trains
Quantized tensor trains (QTT) offer a way to represent the full solution vector with a number of parameters that grows only polynomially with the number of assets d and polylogarithmically with the grid size N. The key idea is to break the huge vector of length N^d into a product of many small tensors, each of which lives on a binary (quantized) representation of its index.
Concretely, suppose we use c QTT cores per spatial dimension, so that N = 2^c. The total number of cores is d × c. Each core has a small bond dimension r (typically between 5 and 20). The total number of parameters is then O(d × c × r^2). For d=5, c=8, and r=15, that is about 5 × 8 × 15^2 = 9000 parameters — a far cry from the trillion points of the full grid.
This compression is not just a storage trick. The paper shows that the Black–Scholes operator itself has an exact QTT representation with rank bounded by O(d), and that the solution can be computed by solving linear systems entirely within the QTT format. The result is a solver that runs on a personal computer for d=3,4,5 assets, producing full-grid prices and Greeks with high accuracy.
What This Tutorial Covers
In the sections that follow, we will build this QTT solver from the ground up. We will start with the core data structure (the quantized tensor train) and its basic arithmetic. Then we will construct the analytic building blocks — the exponential function, boundary selectors, and tridiagonal matrices — that form the operator and payoff. We will implement the TT-Cross algorithm for approximating the max and min functions needed for payoffs and early exercise. We will develop the ALS/MALS solver for linear systems in QTT format. Finally, we will assemble everything into the time-stepping and space-time solvers, and compute Greeks from the solution.
By the end, you will understand not only how the paper achieves its remarkable results, but also how to implement a QTT-based PDE solver yourself.
QTT Solver Methodology
This section details the core algorithmic components of the quantized tensor train (QTT) solvers for multi-asset option pricing. The key idea is to represent the discretized solution and operators as QTT tensors, enabling efficient storage and computation that scales polynomially with the number of assets.
Quantized Tensor Trains (QTT)
A tensor train (TT) decomposes a d-dimensional tensor into a product of 3-dimensional core tensors. The QTT format further applies this decomposition to a quantized (binary) representation of each dimension, effectively converting an N-point grid into a tensor of order log2(N). For a grid with 2^c points per dimension, the QTT representation uses c cores per dimension, each of small rank.
TT-Cross Algorithm
The TT-Cross algorithm approximates a tensor by sampling its entries along selected fibers (rows/columns). It iteratively refines the approximation by finding the most informative cross-sections, requiring only O(d r^2 N) evaluations of the original tensor function, where r is the TT rank.
Alternating Linear Scheme (ALS)
ALS solves linear systems in the TT format by optimizing one core at a time while keeping others fixed. It cycles through all cores, solving local least-squares problems, and converges to a solution of the global system. The Multigrid ALS (MALS) variant improves convergence by adaptively adjusting ranks.
Time-Stepping Solver Pipeline
Discretization: Apply backward Euler in time and central differences in space to the log-transformed Black-Scholes PDE.
QTT Construction: Build the finite-difference operator and initial condition (payoff) as QTT tensors with known low-rank structures.
Time Stepping: For each time step, solve the linear system (I - Δt A) v{n+1} = vn using ALS, where A is the spatial discretization operator.
Boundary Conditions: Impose Dirichlet conditions using an eraser MPO that zeros out boundary points.
Output: The solution at each time step is a QTT tensor representing the option price on the full grid.
Space-Time Formulation
An alternative approach solves for all time steps simultaneously by forming a block-bidiagonal system that couples spatial points across time. This yields a single large QTT system that can be solved with ALS, avoiding sequential time stepping. The space-time operator has a known QTT rank bound that grows only polynomially with the number of assets.
Comparison
Time-stepping: More flexible (handles American options via early exercise checks), but requires solving a linear system at each step.
Space-time: Faster for European options (single solve), but cannot easily incorporate early exercise.
Both methods achieve full-grid solutions with polynomial scaling in the number of assets and polylogarithmic scaling in grid size.
Building Blocks: Analytic QTT Representations
Now that we have a feel for the QTT data structure, we need the concrete building blocks that the multi-asset solver uses. The Black–Scholes operator, its boundary conditions, and the payoff functions all rely on a few special functions and matrices that can be represented exactly and with tiny rank in the QTT format. This section implements three of those building blocks: the exponential function (rank 1), the left/right boundary selection vectors (rank 2), and the tridiagonal Toeplitz matrix (rank 3). These constructions follow Appendices C.1, C.2, and C.3 of the paper.
Exponential Function (rank 1)
Why is the exponential function important? In the log-price coordinates x = ln S, the payoff of a basket option involves terms like e^x. The QTT representation of e^{αx} on a dyadic grid is exact and has rank 1. The key insight is the binary expansion of the grid index: on the interval (0,1), a grid point x is represented by a binary fraction 0.b_1 b_2 ... b_c, where b_i is the i-th bit of the index. Then
e^{αx} = e^{α * Σ b_i * 2^{-i}} = ∏_i e^{α * b_i * 2^{-i}}.
This product form translates directly into a QTT with cores F_i where F_i[0] = 1 and F_i[1] = exp(α * 2^{-i}). The function analytic_qtt_exponential in qtt_analytic.py builds exactly this representation. The core loop is:
cores = []
for i in range(c):
factor = np.exp(scaled_alpha * 2.0 ** (-(i + 1)))
core = np.array([[[1.0]], [[factor]]], dtype=np.float64)
core = core.reshape(1, 2, 1)
cores.append(core)Each core has shape (1, 2, 1) – a rank-1 MPS. To map to an arbitrary interval (a, b), we first scale alpha by (b-a), then multiply every entry of the first core by e^a. The result is a QTT vector of length 2^c that matches the exponential function at every grid point.
Boundary Selection Vectors vL and vR (rank 2)
When imposing Dirichlet boundary conditions, we need to modify the solution at the leftmost and rightmost grid points. The paper defines two vectors: vL has 1 at every position except the first (index 0), where it is 0; vR has 1 at every position except the last (index 2^c-1), where it is 0. These are building blocks for the “eraser” MPO that zeroes out boundary rows and columns.
The QTT representation of vL and vR is exact and has rank 2. The construction of the first core is instructive. For vL:
first = np.zeros((1, 2, 2), dtype=np.float64)
first[0, 0, 0] = 1.0
first[0, 1, 1] = 1.0For vL, the first core distinguishes between the grid index 0 (bit 0) and all others (bit 1). The middle cores propagate the two “states” (leftmost or not) through the bit string, and the last core enforces the final zero or one. The result is a rank-2 MPS with cores of shape (1,2,2), (2,2,2), ..., (2,2,1). The function qtt_vL_vR(c, side='L') returns this representation.
Tridiagonal Toeplitz Matrix (rank 3) – Lemma 1
The finite-difference discretisation of the Black–Scholes operator requires tridiagonal matrices for the second and first derivative terms. These matrices are Toeplitz (constant diagonals) and have size 2^c × 2^c. Lemma 1 in the paper (Appendix C.3) gives an exact QTT representation with bond dimension 3.
Let the matrix have diagonal entries α, superdiagonal entries β, and subdiagonal entries γ. The QTT uses three 2×2 building blocks:
I = np.array([[1.0, 0.0], [0.0, 1.0]], dtype=np.float64)
J = np.array([[0.0, 1.0], [0.0, 0.0]], dtype=np.float64)
Jp = np.array([[0.0, 0.0], [1.0, 0.0]], dtype=np.float64) # J'The first core has shape (1, 2, 2, 3) and its three slices are αI, βJ, and γJ':
first = np.zeros((1, 2, 2, 3), dtype=np.float64)
first[0, :, :, 0] = alpha * I
first[0, :, :, 1] = beta * J
first[0, :, :, 2] = gamma * JpMiddle cores (shape (3, 2, 2, 3)) repeat the same pattern: each of the three left bond indices corresponds to a “state” that routes the three diagonal contributions. The last core (shape (3, 2, 2, 1)) again contains αI, βJ, and γJ' but now projects to a single right bond. The function qtt_toeplitz_tridiagonal(c, alpha, beta, gamma) returns this rank-3 MPO.
For completeness, qtt_identity_mpo(c) provides the identity matrix as a rank-1 MPO, simply formed by repeating the 2×2 identity reshaped to (1,2,2,1) for each core.
Worked Example: Exponential on a 3-Core Grid
Take c=3 cores (grid of 8 points in (0,1)), α=0.5. The cores are:
Core 0:
[1, exp(0.5 * 0.5)] ≈ [1, 1.284]Core 1:
[1, exp(0.5 * 0.25)] ≈ [1, 1.133]Core 2:
[1, exp(0.5 * 0.125)] ≈ [1, 1.064]
Contracting the cores (each of shape (1,2,1)) gives the vector of length 8: [1, 1.064, 1.133, 1.203, 1.284, 1.365, 1.455, 1.549].
This matches the exact values of e^{0.5 x} at the points x = 0, 1/8, ..., 7/8.
For the tridiagonal Toeplitz matrix with α=2, β=-1, γ=-1 (the standard second-derivative stencil), the QTT representation has rank 3 and is exact. When applied to a vector, it produces the same result as the full matrix-vector product, up to machine precision.
Advanced Notes
The exponential QTT cores F_i[0] = 1, F_i[1] = exp(α * 2^{-i}) follow directly from the binary expansion of the grid index. Because the product of the cores reproduces the product of the exponentials, the representation is exact and rank 1. The same idea extends to higher dimensions: for a d-dimensional tensor product grid, the Kronecker product of the 1D QTTs yields a rank-1 d-dimensional MPS.
The tridiagonal Toeplitz representation (Lemma 1) uses the three 2×2 matrices I, J, J' to encode the three non-zero diagonals. The bond dimension of 3 arises because the matrix has three independent “streams” – the diagonal, the superdiagonal, and the subdiagonal – that must be carried through the tensor network. The middle cores are all identical, which makes the construction highly efficient for large c.
TT-Cross: Approximating Functions with Controlled Rank
Use the button or URL below to download the source code.
In the previous sections we built exact, low-rank QTT representations for the exponential function, boundary vectors, and tridiagonal matrices. But the payoff of an option — max(K - S, 0) for a put, or max(S_1, S_2, ..., S_d) for a max-min contract — is not an exponential or a linear function. It is a piecewise-linear function with a kink, and its QTT rank can be large if we try to represent it exactly. The same problem appears in the early-exercise step for American options, where we need to compute the element-wise maximum of the solution and the payoff at every time step.
To handle these operations, the paper uses the TT-Cross algorithm. Instead of building an exact representation, TT-Cross approximates a function on a tensor grid with a controlled rank. It does this by evaluating the function on a carefully chosen subset of grid points — the cross — and then constructing a low-rank tensor that matches the function at those points. The result is a QTT tensor whose bond dimensions never exceed a user-specified cap χ, and whose accuracy is sufficient for pricing errors below 1%.
The Cross Approximation Idea
Suppose you have a matrix A of size m × n and you want a low-rank approximation. The skeleton decomposition says that if A has rank r, you can pick r rows and r columns such that the submatrix at their intersection (the cross) determines the whole matrix. Concretely, if you select row indices I and column indices J, then
A ≈ A[:, J] * (A[I, J])^{-1} * A[I, :].
For a tensor, the same idea extends: you select fibers (the tensor analogue of rows and columns) along each dimension, and build a low-rank tensor train that interpolates the original tensor on those fibers. The TT-Cross algorithm automates the selection of these fibers by sweeping through the dimensions, using a routine called maxvol to find the most informative indices.
The tt_cross Function
The file tt_cross.py implements the TT-Cross algorithm for general function approximation on a tensor grid. The main entry point is:
def tt_cross(function_handle, shape, max_rank, max_sweeps=4, tol=1e-6):
"""TT-Cross algorithm for approximating a function on a tensor grid.
Parameters
----------
function_handle : callable
A function that takes a tuple of integer indices (one per dimension)
and returns a scalar float.
shape : list of int
Shape of the tensor grid. Each entry must be a power of two.
max_rank : int
Maximum allowed bond dimension (rank) of the approximation.
max_sweeps : int
Maximum number of alternating sweeps.
tol : float
Tolerance for early stopping.
Returns
-------
QTTTensor
A QTT tensor (MPS) approximating the function on the given grid.
"""The algorithm works as follows:
Initialise a random tensor train with bond dimensions capped at
max_rank.Left-to-right sweep: for each dimension
k, fix all cores except thek-th. Evaluate the function on a set of indices selected by the maxvol procedure on the unfolding matrix. Update thek-th core to match those evaluations.Right-to-left sweep: repeat the process in the opposite direction.
Iterate for
max_sweepssweeps, or until the approximation stops changing.
The function handle is evaluated only on the selected fibers, not on the entire grid. This is what makes the algorithm efficient for high-dimensional problems: the number of evaluations is proportional to d * χ^2, where d is the number of dimensions and χ is the target rank.
The Maxvol Routine
The maxvol routine finds a set of rows (or columns) that form a submatrix with large determinant. This is a greedy algorithm: it starts with a random set of rows, then iteratively swaps rows to increase the volume (absolute determinant) of the submatrix. The implementation in tt_cross.py is:
def _maxvol(A, tol=1e-12):
"""Find a set of rows that form a maximum-volume submatrix of A."""
m, r = A.shape
if m <= r:
return np.arange(m)
ind = np.arange(r)
B = A[ind, :]
for _ in range(10):
try:
B_inv = np.linalg.inv(B)
except np.linalg.LinAlgError:
break
C = A @ B_inv
i_max, j_max = np.unravel_index(np.argmax(np.abs(C)), C.shape)
if np.abs(C[i_max, j_max]) <= 1.0 + tol:
break
ind[j_max] = i_max
B = A[ind, :]
return indThis function takes a matrix A of shape (m, r) (with m >= r) and returns r row indices that approximately maximise the volume. The algorithm is simple: it computes the pivot matrix C = A * B^{-1}, finds the element of C with the largest absolute value, and swaps that row into the selected set. The loop continues until no element exceeds 1 + tol in absolute value, indicating that the current submatrix is close to maximum volume.
Element-Wise Max and Min of Two QTT Tensors
The most important use of TT-Cross in the solver is for the element-wise maximum and minimum operations. These are implemented as tt_cross_max and tt_cross_min:
def tt_cross_max(qtt_a, qtt_b, max_rank, max_sweeps=4):
"""Compute element-wise maximum of two QTT tensors via TT-Cross."""
shape = qtt_a.shape()
if shape != qtt_b.shape():
raise ValueError("Input tensors must have the same shape")
# For small tensors, convert to full and use numpy
total_size = 1
for s in shape:
total_size *= s
if total_size <= 10**6:
full_a = qtt_a.to_full()
full_b = qtt_b.to_full()
full_max = np.maximum(full_a, full_b)
# Convert back to QTT via TT-SVD
...
else:
# Use TT-Cross with function handle = max(f_a, f_b)
def max_func(idx):
return max(qtt_a.evaluate(idx), qtt_b.evaluate(idx))
return tt_cross(max_func, shape, max_rank, max_sweeps)The function first checks if the total size of the tensor is small enough to convert to a full array. If so, it computes the max directly using NumPy and then converts back to QTT via TT-SVD. For larger tensors, it defines a function handle that evaluates the max at a given multi-index by evaluating both input tensors at that index, and then calls tt_cross with that handle. The tt_cross_min function works identically, using min instead of max.
Worked Example: Basket Put Payoff
Consider a two-asset basket put with strike K = 34 and equal weights. The payoff at maturity is:
V(S_1, S_2) = max(K - (S_1 + S_2), 0).
In log-price coordinates x_i = ln S_i, this becomes:
V(x_1, x_2) = max(K - (e^{x_1} + e^{x_2}), 0).
We want to represent this payoff as a QTT tensor on a grid with c = 3 cores per dimension, giving 2^3 = 8 points per dimension and a total of 64 grid points. The function is not low-rank: the max operation introduces a kink along the line e^{x_1} + e^{x_2} = K. However, we can approximate it using TT-Cross with a rank cap of χ = 10 and 4 sweeps.
import numpy as np
from tt_cross import tt_cross
c = 3
n = 2**c # 8 points per dimension
K = 34.0
# Define the payoff function on the grid
def basket_put_payoff(idx):
# idx is a tuple (i1, i2) with 0 <= i_k < n
# Map indices to log-prices on a uniform grid in [-2, 2]
x1 = -2.0 + 4.0 * idx[0] / (n - 1)
x2 = -2.0 + 4.0 * idx[1] / (n - 1)
S1 = np.exp(x1)
S2 = np.exp(x2)
return max(K - (S1 + S2), 0.0)
# Run TT-Cross
payoff_qtt = tt_cross(basket_put_payoff, [n, n], max_rank=10, max_sweeps=4)The resulting payoff_qtt is a QTT tensor with bond dimensions at most 10. To check the accuracy, we can compare the QTT approximation with the exact payoff at a few random grid points:
# Evaluate at a random point
idx = (3, 5)
exact = basket_put_payoff(idx)
approx = payoff_qtt.evaluate(idx)
print(f"Exact: {exact:.4f}, Approx: {approx:.4f}, Error: {abs(exact-approx):.4e}")For a rank cap of 10 and 4 sweeps, the relative error is typically on the order of 10^{-3}, which is sufficient for pricing errors below 1%.
Advanced Detail: Why TT-Cross Works for the Max Function
The max function is not low-rank, but it is approximately low-rank in the QTT format. The reason is that the kink (the set of points where e^{x_1} + e^{x_2} = K) is a smooth curve in the log-price domain, and the QTT representation can capture this structure with moderate rank. The TT-Cross algorithm exploits this by sampling the function along fibers that cross the kink, building an approximation that is accurate everywhere.
For American options, the early-exercise condition requires computing max(V, payoff) at every time step. The paper uses tt_cross_max followed by rank reduction via qtt_round to keep the bond dimensions manageable. The experiments in Section V.C show that the rank of the solution grows only modestly (from about 10 to about 20) over the life of the option, and the additional cost of the early-exercise step is about 30% of the total runtime.
Summary
TT-Cross approximates a function on a tensor grid with controlled rank by evaluating it on a subset of points (the cross).
The maxvol routine selects the most informative fibers by maximising the volume of the submatrix.
tt_cross_maxandtt_cross_minuse TT-Cross to compute element-wise max/min of two QTT tensors, which is essential for payoffs and early-exercise.The approximation accuracy is controlled by the rank cap
χand the number of sweeps; typical values areχ = 10-20and 3-4 sweeps.For the basket put payoff on a 2D grid with 64 points, TT-Cross with rank 10 achieves relative error ~10^{-3}.
ALS and MALS Solvers for Linear Systems
Once we have the Black–Scholes operator A and the right-hand side b in QTT format, we need to solve the linear system A x = b. This is the computational heart of every time step in the time-stepping solver and of the single global solve in the space-time solver. The paper uses the Alternating Linear Scheme (ALS) and its variant MALS (Modified ALS) to solve these systems entirely within the QTT format, without ever forming the full matrix or vector.
ALS is a coordinate-descent-like algorithm for tensor trains. It sweeps through the cores of the unknown solution x, updating one core at a time while keeping all other cores fixed. At each core, the problem reduces to a small dense linear system — typically 2×2 or 4×4 — that can be solved by direct factorization. MALS extends this by updating two adjacent cores simultaneously, which allows the bond dimension to grow or shrink during the optimization. Both algorithms are described in Section III.D of the paper.
The ALS Sweep
Suppose we have the operator A as a QTT-MPO and the right-hand side b as a QTT-MPS. We want to find the QTT-MPS x that satisfies A x = b. The ALS algorithm proceeds as follows:
Initialize
xwith some initial guess (often the solution from the previous time step, or a zero tensor).Sweep left-to-right: for each core index
ifrom 0 tonum_cores - 1:
Fix all cores of
xexcept corei.Contract the environment: compute the left and right parts of
Aandxthat involve all cores excepti, and also the left and right parts ofb.Form a small local linear system
A_local * x_core_i = b_localby contracting the environment tensors with the current core.Solve the local system (typically 2×2 or 4×4) using a dense solver.
Update core
iwith the solution.
Sweep right-to-left: repeat the same process in reverse order.
Repeat for a fixed number of sweeps (typically 2–4) or until the residual norm
||A x - b|| / ||b||falls below a tolerance.
The key insight is that the local system is tiny because each core has only two physical indices (size 2 each) and the bond dimensions are small (typically 3–20). The environment contraction can be computed efficiently by precomputing left and right environment tensors, as shown in the helper function _contract_left_env.
MALS: Allowing Bond Dimension Changes
ALS keeps the bond dimensions of x fixed during the sweep. This is fine if the initial guess has the right bond dimensions, but it can be restrictive. MALS (Modified ALS) optimizes two adjacent cores at a time, which allows the bond dimension between them to change. The local system becomes slightly larger (e.g., 4×4 instead of 2×2), but the algorithm can adapt the rank of the solution automatically.
Implementation in als_mals.py
The file als_mals.py implements both solvers. The core function als_solve takes the operator A, the right-hand side b, an optional initial guess x0, and the number of sweeps. It returns the solution x as a QTTTensor.
Here is the public interface:
def als_solve(
A: QTTTensor,
b: QTTTensor,
x0: Optional[QTTTensor] = None,
max_sweeps: int = 2,
local_solver: str = 'dense',
verbose: bool = False
) -> QTTTensor:
"""
Solve A x = b using the Alternating Linear Scheme (ALS).
Parameters
----------
A : QTTTensor
System matrix as MPO, shape (2^{c*d}, 2^{c*d}).
b : QTTTensor
Right-hand side as MPS, shape (2^{c*d},).
x0 : QTTTensor, optional
Initial guess (MPS). If None, a zero tensor is used.
max_sweeps : int
Number of full left-right-right-left sweeps.
local_solver : str
Method for solving the local system ('dense' or 'lstsq').
verbose : bool
If True, print residual norm after each sweep.
Returns
-------
x : QTTTensor
Solution MPS.
"""The MALS variant mals_solve has the same signature but uses two-core updates internally.
The Local Solve
The local system is solved by _local_solve, which handles the small dense linear system with a regularization term to avoid ill-conditioning:
def _local_solve(
A_local: np.ndarray,
b_local: np.ndarray,
reg: float = 1e-12
) -> np.ndarray:
"""
Solve a small dense linear system A_local @ x = b_local.
A small regularization term is added to handle potential ill-conditioning.
"""
if A_local.shape[0] == 0:
return np.array([])
try:
return np.linalg.solve(A_local + reg * np.eye(A_local.shape[0]), b_local)
except np.linalg.LinAlgError:
return np.linalg.lstsq(A_local, b_local, rcond=None)[0]This function is called for each pair of bond indices during the sweep. The local matrix A_local is formed by contracting the left and right environment tensors with the current core of A and the current core of x. The local right-hand side b_local is formed similarly from the environment of b.
Worked Example: 1D Tridiagonal System
Consider a 1D problem with 8 grid points (c=3). The operator A is a tridiagonal matrix (e.g., the second derivative stencil) represented as a QTT-MPO with rank 3. The right-hand side b is a vector (e.g., the payoff) represented as a QTT-MPS with rank 5. We want to solve A x = b.
Initialization:
x0is a zero MPS with rank 1.Sweep 1 (left-to-right):
Core 0: the left environment is empty (scalar 1). The right environment is contracted from cores 1 and 2 of
Aandx. The local system is 2×2 (since the physical dimension is 2). Solve and update core 0.Core 1: left environment now includes core 0; right environment includes core 2. Local system is 2×2. Solve and update core 1.
Core 2: left environment includes cores 0 and 1; right environment is empty. Local system is 2×2. Solve and update core 2.
Sweep 1 (right-to-left): repeat in reverse order.
After 2–3 sweeps, the residual norm drops below 1e-6, and the solution
xis accurate to machine precision.
The entire solve takes a few milliseconds on a laptop.
Advanced Notes
Environment contraction: The left environment for core
iis a tensor of shape(A_left_bond_i, x_left_bond_i)that represents the contraction of all cores 0..i-1 ofAandx, with all physical indices summed out. The function_contract_left_envcomputes this iteratively. The right environment is computed analogously by contracting from the right.MALS bond dimension growth: When two cores are optimized together, the bond dimension between them can increase up to the product of the two original bond dimensions. After the local solve, the new bond dimension is truncated via SVD to a target rank or tolerance. This allows the algorithm to adapt the rank of the solution automatically.
Convergence: The paper reports that 2 sweeps are sufficient for the multi-asset experiments. The solver stops when the residual norm
||A x - b|| / ||b||falls below a tolerance (default 1e-8) or aftermax_sweepssweeps.
Summary
ALS and MALS are the workhorses of the QTT solver. They solve linear systems in the QTT format by iteratively updating one or two cores at a time, using small dense local solves. The implementation in als_mals.py follows the paper's description and provides both als_solve and mals_solve functions. These functions are called by the time-stepping solver at every time step and by the space-time solver for the global system.
Constructing the Black-Scholes Operator in QTT
We now have all the analytic building blocks we need: the exponential function (rank 1), the boundary selection vectors (rank 2), and the tridiagonal Toeplitz matrix (rank 3). The next step is to assemble these pieces into the d-asset Black–Scholes spatial operator in QTT format. This operator is the core of every time step: at each step we solve (I - Δt A) x_{n+1} = x_n, where A is the spatial discretization of the Black–Scholes PDE.
The key insight is that the Black–Scholes PDE in log-price coordinates has constant coefficients. This means the spatial operator is a sum of terms that are tensor products of one-dimensional operators. Each one-dimensional operator — second derivative, first derivative, identity — can be represented exactly in QTT with rank at most 3. By combining them via the Kronecker product and addition, we obtain the full d-asset operator with a rank that grows only linearly with d.
The d-Asset Black–Scholes Operator in Log-Price Coordinates
Recall the Black–Scholes PDE after the log transformation x_i = ln S_i and time reversal τ = T - t:
This is equation (A2) from the paper, extended to d assets. The spatial operator L is the part acting on the spatial variables x_1, ..., x_d:
where I is the identity operator. The time-stepping solver uses the operator A = I - Δt L, which appears in the backward Euler discretization.
Decomposing L into Tensor Products of 1D Operators
The crucial observation is that the second-derivative terms can be split into two groups:
Diagonal terms (
i = j):(1/2) σ_i^2 ∂²/∂x_i²Cross terms (
i ≠ j):ρ_{ij} σ_i σ_j ∂²/(∂x_i ∂x_j)
Each of these terms acts on a product of one-dimensional operators. For example, the operator ∂²/(∂x_i ∂x_j) is the tensor product of the first-derivative operator in dimension i and the first-derivative operator in dimension j. Similarly, ∂²/∂x_i² is the tensor product of the second-derivative operator in dimension i and identity operators in all other dimensions.
Formally, let D2_i be the second-derivative matrix (tridiagonal Toeplitz) for dimension i, let D1_i be the first-derivative matrix (skew-symmetric tridiagonal), and let I_i be the identity matrix. All three are of size 2^c × 2^c and have exact QTT representations (rank 3 for D2_i and D1_i, rank 1 for I_i). Then the d-dimensional operator L can be written as:
Each term in this sum is a Kronecker product of d one-dimensional MPOs. The Kronecker product of MPOs is itself an MPO: we simply concatenate the cores of the factor MPOs, with the bond dimensions multiplying. Since each factor has rank at most 3, the rank of each term is at most 3^d in the worst case. However, the paper shows that the sum of these terms has a rank bounded by O(d) because the terms share a common structure and can be combined efficiently. In practice, the rank of the full operator A is typically 3d + 1 or less.
Building the Operator in Code
Let us see how this decomposition translates into code. The file bs_operator.py contains the function build_bs_operator that constructs the d-asset operator as a QTT-MPO.
# bs_operator.py (excerpt)
def build_bs_operator(d, c, sigma, rho, r, dt):
"""
Build the QTT-MPO for A = I - dt * L, where L is the d-asset
Black-Scholes spatial operator in log-price coordinates.
Parameters
----------
d : int
Number of assets.
c : int
Number of QTT cores per spatial dimension (grid points = 2^c).
sigma : list of float, length d
Volatilities.
rho : numpy.ndarray, shape (d, d)
Correlation matrix.
r : float
Risk-free rate.
dt : float
Time step size.
Returns
-------
QTTTensor
The operator A as a QTT-MPO.
"""
# Build 1D operators
I_1d = qtt_identity_mpo(c) # rank 1
D2_1d = [qtt_toeplitz_tridiagonal(c, alpha=-2.0, beta=1.0, gamma=1.0)
for _ in range(d)] # rank 3, second derivative stencil
D1_1d = [qtt_toeplitz_tridiagonal(c, alpha=0.0, beta=0.5, gamma=-0.5)
for _ in range(d)] # rank 3, first derivative stencil
# Start with the identity term: -r * I_1 ⊗ ... ⊗ I_d
A = qtt_mul_scalar(qtt_kronecker_product([I_1d] * d), -r)
# Add diagonal second-derivative terms
for i in range(d):
factors = [I_1d] * d
factors[i] = D2_1d[i]
term = qtt_kronecker_product(factors)
term = qtt_mul_scalar(term, 0.5 * sigma[i]**2)
A = qtt_add(A, term)
# Add first-derivative terms
for i in range(d):
factors = [I_1d] * d
factors[i] = D1_1d[i]
term = qtt_kronecker_product(factors)
term = qtt_mul_scalar(term, r - 0.5 * sigma[i]**2)
A = qtt_add(A, term)
# Add cross-derivative terms
for i in range(d):
for j in range(i+1, d):
factors = [I_1d] * d
factors[i] = D1_1d[i]
factors[j] = D1_1d[j]
term = qtt_kronecker_product(factors)
term = qtt_mul_scalar(term, rho[i, j] * sigma[i] * sigma[j])
A = qtt_add(A, term)
# Form A = I - dt * L
A = qtt_mul_scalar(A, -dt)
A = qtt_add(qtt_kronecker_product([I_1d] * d), A)
# Round to control rank growth
A = qtt_round(A, tol=1e-12)
return ANotice the structure: we build each term as a Kronecker product of 1D MPOs, scale it by the appropriate coefficient, and add it to the accumulator A. The final step forms I - dt * L and rounds the result to keep the bond dimensions under control. The rounding tolerance 1e-12 is tight enough to preserve the exactness of the analytic building blocks.
The Eraser MPO: Imposing Boundary Conditions
The operator A acts on the entire grid, including boundary points. However, the Black–Scholes PDE requires Dirichlet boundary conditions on all faces of the domain. The paper enforces these conditions by constructing a modified right-hand side at each time step, rather than modifying the operator itself. The key tool is the eraser MPO, which zeros out the rows and columns of the operator that correspond to boundary points.
The eraser MPO is built from the vL and vR vectors we constructed in the analytic QTT section. For each dimension, we create an MPO that selects the interior points: E_i = I_i - vL_i ⊗ vL_i^T - vR_i ⊗ vR_i^T, where vL_i is the vector that is 1 everywhere except the first entry (which is 0), and vR_i is the vector that is 1 everywhere except the last entry. The full eraser MPO is the Kronecker product of these 1D erasers:
# bs_operator.py (excerpt)
def eraser_mpo(d, c):
"""
Build the eraser MPO that zeros out all boundary elements.
Parameters
----------
d : int
Number of assets.
c : int
Number of QTT cores per spatial dimension.
Returns
-------
QTTTensor
The eraser MPO of shape (2^{c*d}, 2^{c*d}).
"""
I_1d = qtt_identity_mpo(c)
vL = qtt_vL_vR(c, side='L') # rank-2 MPS
vR = qtt_vL_vR(c, side='R') # rank-2 MPS
# Build 1D eraser: I - vL vL^T - vR vR^T
# vL vL^T is an MPO of rank 2 (outer product of vL with itself)
vL_mpo = qtt_outer_product(vL, vL) # rank-2 MPO
vR_mpo = qtt_outer_product(vR, vR) # rank-2 MPO
E_1d = qtt_add(I_1d, qtt_mul_scalar(vL_mpo, -1.0))
E_1d = qtt_add(E_1d, qtt_mul_scalar(vR_mpo, -1.0))
E_1d = qtt_round(E_1d, tol=1e-12)
# Kronecker product over all dimensions
E = qtt_kronecker_product([E_1d] * d)
E = qtt_round(E, tol=1e-12)
return EThe function qtt_outer_product (not shown here) takes an MPS and returns an MPO representing the outer product with itself. The resulting eraser MPO, when applied to a vector, zeros out all entries that lie on any boundary face. In the time-stepping solver, we apply the eraser to the solution after each time step to enforce the Dirichlet conditions.
Worked Example: Two-Asset Operator
Let us walk through a concrete example. Suppose we have two assets with σ₁ = 0.25, σ₂ = 0.15, ρ₁₂ = 0.4, r = 0.05, and we use c = 3 cores per dimension (8 grid points per asset). The operator L has the following terms:
Diagonal second derivatives:
(0.25²/2) I ⊗ D2₁ + (0.15²/2) D2₂ ⊗ IFirst derivatives:
(0.05 - 0.25²/2) I ⊗ D1₁ + (0.05 - 0.15²/2) D1₂ ⊗ ICross derivative:
0.4 * 0.25 * 0.15 * D1₁ ⊗ D1₂Identity:
-0.05 I ⊗ I
Each D2 and D1 is a rank-3 MPO of size 8 × 8. The Kronecker product I ⊗ D2₁ is an MPO of size 64 × 64 with rank 3 (since the identity has rank 1). The cross term D1₁ ⊗ D1₂ has rank up to 9 (3 × 3), but after addition and rounding, the total rank of A = I - Δt L is typically around 7 for d = 2. This is far smaller than the full matrix size of 64 × 64 = 4096 entries.
Advanced Notes
Rank bound: The paper proves that the rank of the d-asset operator
Ais bounded by3d + 1. This is because each of theddiagonal second-derivative terms contributes at most 1 to the rank, each of thedfirst-derivative terms contributes at most 1, and thed(d-1)/2cross-derivative terms can be combined into a single rank-dterm. The identity adds 1. In practice, the rank is often lower due to the structure of the correlation matrix.Boundary conditions: The eraser MPO approach is elegant but requires an additional MPO-MPS contraction at each time step. An alternative is to incorporate the boundary conditions directly into the operator by modifying the cores of
A. The paper uses the eraser approach for simplicity.Time-dependent coefficients: The operator
Ais constant in time because the Black–Scholes PDE has constant coefficients. This means we buildAonce and reuse it at every time step, saving significant computation.
With the operator A and the eraser MPO in hand, we are ready to assemble the full time-stepping solver. The next section will show how to combine these pieces with the payoff and boundary conditions to price European and American options.
Time-Stepping Solver for European and American Options
We now have all the pieces in place: the Black–Scholes operator A as a QTT-MPO, the payoff as a QTT-MPS, and the ALS solver to solve linear systems. The final step is to assemble these into a time-stepping algorithm that advances the solution from maturity backwards to the present. This section implements the backward Euler time-stepping solver for both European and American options, following the algorithm described in Appendix A.5 of the paper.
The idea is straightforward. Starting from the payoff at maturity (τ = 0), we repeatedly solve the linear system (I - Δt A) x_{n+1} = x_n to step backwards in time. For European options, this is all we do. For American options, after each time step we also enforce the early-exercise condition by taking the element-wise maximum of the solution and the payoff. The ALS solver handles the linear solve, and the TT-Cross algorithm handles the max operation.
The Backward Euler Discretization
After the log transformation x_i = ln S_i and time reversal τ = T - t, the Black–Scholes PDE becomes:
We discretize this in time using the fully implicit (backward Euler) scheme. Let Δt be the time step size, and let x_n denote the solution at time level n (with x_0 being the payoff at maturity). The backward Euler discretization gives:
where L is the spatial operator (the sum of derivative terms). Rearranging:
This is the linear system we solve at each time step. The matrix I - Δt L is the same for every step (since L is constant in time), so we build it once and reuse it. The right-hand side x_n changes each step, but the ALS solver can use the previous solution as an initial guess, which speeds up convergence.
The Time-Stepping Algorithm
The algorithm proceeds as follows:
Build the operator
A = I - Δt Las a QTT-MPO using the construction from the previous section.Build the payoff
pas a QTT-MPS using the payoff construction (basket or max-min).Set
x = p(solution at maturity).For each time step from
n = 0toN-1:
Solve
A x_new = xusing ALS, withxas the initial guess.Set
x = x_new.If American: compute
x = max(x, p)element-wise usingtt_cross_max, then apply rank reduction viaqtt_round.
Return
xas the solution att = 0.
The number of time steps N is typically 2^c, matching the spatial grid size. This ensures that the temporal discretization error is of the same order as the spatial error.
Code Walkthrough: time_stepping_solver
The implementation lives in solver.py. The main function is time_stepping_solver, which orchestrates the loop. Here is the core logic:
def time_stepping_solver(market_params, grid_params, option_type, contract_type):
# Unpack parameters
d = market_params['d']
c = grid_params['c']
dt = grid_params.get('dt', market_params['T'] / 2**c)
# Build operator A = I - dt * L
A = build_operator_A(d, c, market_params['sigma'], market_params['rho'],
market_params['r'], dt)
# Build payoff
if contract_type == 'basket':
payoff = payoff_basket_put(market_params['K'], market_params['weights'],
market_params['sigma'], market_params['rho'],
market_params['r'], market_params['T'],
c, grid_params['domain'])
elif contract_type == 'maxmin':
payoff = payoff_maxmin_put(market_params['K'], market_params['sigma'],
market_params['rho'], market_params['r'],
market_params['T'], c, grid_params['domain'])
else:
raise ValueError(f"Unknown contract type: {contract_type}")
# Initialize solution with payoff
x = payoff
# Time-stepping loop
num_steps = int(market_params['T'] / dt)
for n in range(num_steps):
# Solve (I - dt*L) x_new = x
x_new = als_solve(A, x, x0=x, max_sweeps=2)
# Apply boundary conditions
tau = (n + 1) * dt # time from maturity
x_new = impose_boundaries(x_new, tau, contract_type, market_params,
c, grid_params['domain'])
# American early exercise
if option_type == 'American':
x_new = tt_cross_max(x_new, payoff, max_rank=20, max_sweeps=4)
x_new = qtt_round(x_new, tol=1e-6)
x = x_new
return xWhat to notice:
The operator
Ais built once before the loop. This is efficient becauseAdoes not depend on time.The ALS solver uses the previous solution
xas the initial guess for the next step. Since the solution changes slowly between time steps, this reduces the number of ALS sweeps needed (typically 2 sweeps per step).Boundary conditions are imposed after each solve using
impose_boundaries, which applies the eraser MPO and adds the boundary values.For American options, the early-exercise condition is applied via
tt_cross_maxfollowed by rank reduction. The rank cap of 20 is based on the paper's experiments, which show that ranks remain bounded even for American options.
The European Solver
For European options, the early-exercise step is omitted. The function european_solver is a thin wrapper around the loop:
def european_solver(market_params, grid_params, contract_type):
return time_stepping_solver(market_params, grid_params,
option_type='European',
contract_type=contract_type)The American Solver
For American options, the early-exercise step is included. The function american_solver is similarly a wrapper:
def american_solver(market_params, grid_params, contract_type):
return time_stepping_solver(market_params, grid_params,
option_type='American',
contract_type=contract_type)The key difference is the tt_cross_max and qtt_round calls inside the loop. The paper reports that the early-exercise step adds about 30% to the total run time, but the ranks remain controlled (typically < 20).
Worked Example: 1-Asset European Call
Let's walk through a concrete example. For a single-asset European call with S0 = 10, K = 10, T = 1, r = 0.05, σ = 0.25, and c = 8 cores (256 grid points), the solver proceeds as follows:
Build the operator
A = I - Δt LwhereLis the 1D Black–Scholes operator. The operator is a tridiagonal matrix of size 256×256, represented as a QTT-MPO with rank 3.Build the payoff
p = max(e^x - K, 0)as a QTT-MPS. The payoff is constructed using the analytic exponential QTT and TT-Cross for the max function.Set
x = p.For each of the 256 time steps:
Solve
A x_new = xusing ALS (2 sweeps).Apply boundary conditions (Dirichlet at the boundaries).
Set
x = x_new.
Return
xas the option price att = 0.
The entire computation takes about 1 second on a modern laptop and produces prices accurate to within 0.1% of the analytical Black–Scholes formula.
Advanced Notes
Rank management: The solution
xchanges slowly between time steps, so its QTT rank remains bounded (typically < 10 for European options). For American options, the early-exercise step can increase the rank, but the subsequentqtt_roundcall keeps it under control.Boundary conditions: The
impose_boundariesfunction applies the eraser MPO to zero out boundary points, then adds the correct boundary values. For basket puts, the lower boundaries aremax(Ke^{-rt} - sum of other assets, 0)and the upper boundaries are zero. These boundary values are constructed using the analytic exponential QTT and TT-Cross.Comparison with space-time solver: The time-stepping solver requires one ALS solve per time step, which can be expensive for many time steps. The space-time solver (covered in the next section) solves all time steps simultaneously, which can be faster for European options. However, the time-stepping solver is the only option for American options, since the early-exercise condition must be applied sequentially.
Space-Time Solver for European Options
The time-stepping solver we built in the previous section works well, but it has a drawback: it requires one ALS solve per time step. For a grid with 256 time steps, that is 256 separate linear system solves. The space-time solver takes a different approach: it treats the entire space-time grid as a single, larger linear system and solves it in one go. This can be more efficient because the ALS solver only needs to converge once, and the QTT representation of the space-time operator has a special structure that keeps the ranks low.
The space-time solver is described in Appendix A.1 of the paper for the single-asset case and generalized to multiple assets in Section III.C. It is only applicable to European options because the early-exercise condition for American options would break the linear structure of the space-time system.
The Space-Time System
Recall the Black–Scholes PDE after the log transformation x = ln S and time reversal τ = T - t:
This is the same PDE we used for the time-stepping solver. The difference is in how we discretize the time derivative. Instead of stepping one time level at a time, we discretize all time levels simultaneously using backward Euler. Let N = 2^c be the number of time steps (we use the same c for time as for each spatial dimension). Let V_n be the solution vector at time level n (size 2^{c*d}). The backward Euler discretization gives:
(I - Δt A) V_{n+1} = V_n
where A is the spatial operator (the same QTT-MPO we built in the previous section) and Δt = T / N. This is a recurrence that links each time level to the next. We can write all N equations together as a single block-bidiagonal system:
[ I ] [ V_0 ] [ b_0 ]
[ -I (I - Δt A) ] [ V_1 ] [ 0 ]
[ -I (I - Δt A) ] [ V_2 ] = [ 0 ]
[ ... ... ] [ ... ] [ ... ]
[ -I (I - Δt A) ] [ V_N ] [ 0 ]Here V_0 is the solution at τ = 0 (maturity), which is the payoff. The first row says V_0 = b_0, where b_0 is the payoff vector. The remaining rows enforce the backward Euler recurrence. The system matrix is block-bidiagonal: the diagonal blocks are (I - Δt A) (except the first, which is I), and the subdiagonal blocks are -I.
Building the Space-Time Operator in QTT
The key to the space-time solver is to represent this block-bidiagonal matrix as a QTT-MPO. The structure is a Kronecker product of a time-direction operator and the spatial operator. Let D_t be the N × N bidiagonal matrix that represents the backward Euler recurrence:
D_t = [ 1 ]
[ -1 1 ]
[ -1 1 ]
[ ... ... ]
[ -1 1 ]This is a lower bidiagonal matrix with 1 on the diagonal and -1 on the subdiagonal. The space-time operator A_st is then:
A_st = I_t ⊗ I_x - D_t ⊗ (Δt A_x)
where I_t is the N × N identity, I_x is the 2^{c*d} × 2^{c*d} spatial identity, and A_x is the spatial operator. The first term I_t ⊗ I_x is the identity on the full space-time grid. The second term D_t ⊗ (Δt A_x) encodes the backward Euler coupling between time levels.
Both I_t and D_t are one-dimensional matrices of size N = 2^c. They can be represented in QTT format using the same analytic building blocks we used for the spatial operator. The identity I_t is a rank-1 QTT-MPO (each core is the 2×2 identity). The bidiagonal matrix D_t is a Toeplitz matrix with diagonal 1 and subdiagonal -1, which has a rank-3 QTT representation via Lemma 1 (the same lemma we used for the tridiagonal Toeplitz matrix).
The Kronecker product I_t ⊗ I_x is built by concatenating the cores of I_t and I_x. Similarly, D_t ⊗ (Δt A_x) is built by concatenating the cores of D_t and the scaled spatial operator. The resulting space-time operator A_st is a QTT-MPO with c cores for time and c*d cores for space, for a total of c*(d+1) cores. Its rank is bounded by the sum of the ranks of the two terms, which is at most 1 + 3 * (rank of A_x). Since the spatial operator rank is O(d), the space-time operator rank is also O(d).
The Right-Hand Side
The right-hand side b of the space-time system is a vector of length N * 2^{c*d}. It has the payoff vector b_0 in the first block (corresponding to τ = 0) and zeros everywhere else. In QTT format, this is a tensor that is the Kronecker product of a time-direction vector e_0 (a vector with 1 in the first position and 0 elsewhere) and the payoff QTT-MPS:
b = e_0 ⊗ payoff
The vector e_0 is a one-dimensional vector of length N = 2^c with a single 1 at the first index. This can be represented in QTT format with rank 1: each core has entries [1, 0] (selecting the first bit). The Kronecker product with the payoff QTT-MPS gives a QTT-MPS of length N * 2^{c*d} with rank equal to the payoff rank.
Solving the Space-Time System
Once we have the space-time operator A_st and the right-hand side b in QTT format, we solve the linear system A_st x = b using the ALS or MALS solver. The solution x is a QTT-MPS of length N * 2^{c*d}. It contains the solution at all time levels, stacked in order: the first 2^{c*d} entries are V_0 (the payoff), the next 2^{c*d} entries are V_1 (the solution one time step back), and so on. To extract the solution at the present time (τ = T), we take the last block of 2^{c*d} entries.
The ALS solver converges in a few sweeps because the space-time operator is well-conditioned (it is a small perturbation of the identity). The rank of the solution x is typically bounded by the rank of the right-hand side plus a small increase from the operator application. In practice, the solution rank stays below 20 for the experiments in the paper.
Code Walkthrough: space_time.py
The implementation of the space-time solver is in the file space_time.py. It contains two main functions: build_space_time_operator and space_time_solver.
# space_time.py (excerpt)
def build_space_time_operator(c, dt, A_x):
"""
Build the space-time operator A_st = I_t ⊗ I_x - D_t ⊗ (dt * A_x).
Parameters
----------
c : int
Number of cores for the time dimension (N = 2^c time steps).
dt : float
Time step size.
A_x : QTTTensor (MPO)
Spatial operator (I - dt * A) from the time-stepping solver.
Returns
-------
QTTTensor (MPO)
Space-time operator of size (N * 2^{c*d}, N * 2^{c*d}).
"""
# Build the time-direction identity I_t (rank 1)
I_t = qtt_identity_mpo(c)
# Build the time-direction bidiagonal matrix D_t (rank 3)
# D_t has diagonal 1, subdiagonal -1, superdiagonal 0
D_t = qtt_toeplitz_tridiagonal(c, alpha=1.0, beta=0.0, gamma=-1.0)
# Build the spatial identity I_x
I_x = qtt_identity_mpo(c * d) # d is the number of assets
# Build the scaled spatial operator dt * A_x
dt_A_x = qtt_mul_scalar(A_x, dt)
# Build the two terms
term1 = qtt_kronecker_product([I_t, I_x])
term2 = qtt_kronecker_product([D_t, dt_A_x])
# A_st = term1 - term2
A_st = qtt_add(term1, qtt_mul_scalar(term2, -1.0))
return A_stThis function constructs the space-time operator by combining the time-direction and spatial QTT-MPOs via Kronecker products. The qtt_kronecker_product function concatenates the cores of the input MPOs, producing a new MPO with c + c*d = c*(d+1) cores. The addition and scalar multiplication are the same operations we used for the time-stepping operator.
# space_time.py (excerpt)
def space_time_solver(market_params, grid_params):
"""
Solve the space-time system for a European option.
Parameters
----------
market_params : dict
Market parameters (d, sigma, rho, r, K, T, weights).
grid_params : dict
Grid parameters (c, domain).
Returns
-------
QTTTensor (MPS)
Solution at all time levels, shape (N * 2^{c*d},).
"""
d = market_params['d']
c = grid_params['c']
N = 2**c # number of time steps
dt = market_params['T'] / N
# Build the spatial operator A_x (same as time-stepping)
A_x = build_bs_operator(d, c, market_params['sigma'],
market_params['rho'], market_params['r'], dt)
# Build the space-time operator
A_st = build_space_time_operator(c, dt, A_x)
# Build the payoff at maturity
payoff = payoff_basket_put(market_params['K'], market_params['weights'],
market_params['sigma'], market_params['rho'],
market_params['r'], market_params['T'],
c, grid_params['domain'])
# Build the time-direction vector e_0 (rank 1)
# e_0 has 1 at the first index, 0 elsewhere
e_0_cores = []
for i in range(c):
core = np.zeros((1, 2, 1))
if i == 0:
core[0, 0, 0] = 1.0 # first bit is 0
else:
core[0, 0, 0] = 1.0 # all bits are 0
e_0_cores.append(core)
e_0 = QTTTensor(e_0_cores)
# Right-hand side: b = e_0 ⊗ payoff
b = qtt_kronecker_product([e_0, payoff])
# Solve the space-time system using ALS
x = als_solve(A_st, b, max_sweeps=2)
return xThis function assembles the space-time system and solves it. The key steps are:
Build the spatial operator
A_x(the sameI - Δt Aused in the time-stepping solver).Build the space-time operator
A_stusingbuild_space_time_operator.Build the payoff QTT-MPS.
Build the time-direction vector
e_0(a rank-1 QTT-MPS with a single 1 at the first position).Form the right-hand side
bas the Kronecker product ofe_0and the payoff.Solve the system using ALS.
The solution x is a QTT-MPS of length N * 2^{c*d}. To extract the price at the present time (τ = T), we take the last 2^{c*d} entries. This can be done by reshaping the QTT-MPS or by contracting with a selection vector.
Worked Example: 1-Asset European Call
Let us walk through a concrete example. Consider a 1-asset European call with the following parameters:
Spot price
S0 = 10, strikeK = 10, volatilityσ = 0.25, risk-free rater = 0.05, maturityT = 1.0.Grid:
c = 8cores per dimension, soN = 256time steps and256spatial points.Domain:
xin[-2, 2](log-price).
The space-time system has 256 * 256 = 65536 unknowns. The QTT representation of the solution uses c*(d+1) = 8*2 = 16 cores. With a rank of about 10, the total number of parameters is roughly 16 * 10^2 = 1600 — a tiny fraction of the full grid size.
The ALS solver converges in 2-3 sweeps, taking a few seconds on a personal computer. The resulting price at the spot price matches the analytical Black–Scholes price to within 0.1%.
Comparison with Time-Stepping
The space-time solver has both advantages and disadvantages compared to the time-stepping solver:
Advantage: Only one ALS solve is needed, instead of one per time step. This can be faster for problems with many time steps.
Advantage: The space-time operator has a simple Kronecker structure that is easy to build.
Disadvantage: The space-time system is larger (by a factor of
Nin the number of unknowns), so the ALS solver works with longer cores and potentially higher ranks.Disadvantage: The space-time solver is only applicable to European options. American options require the early-exercise condition, which is a nonlinear operation that cannot be expressed as a linear system.
In practice, the paper reports that the space-time solver is competitive for low dimensions (d ≤ 3) but becomes less efficient than time-stepping for higher dimensions due to rank growth in the space-time operator.
Advanced Notes
The time complexity of the space-time solver for d assets is O(c (d^5 + d^3) χ_ST^3), where χ_ST is the rank of the space-time solution. This is derived in Appendix A.3 of the paper. The cubic dependence on rank means that rank control is critical: if the solution rank grows too large, the solver becomes slow. In practice, the rank stays below 20 for the experiments in the paper, keeping the solver efficient.
The space-time solver can also be extended to handle more complex time discretizations, such as Crank-Nicolson, by modifying the time-direction matrix D_t. However, the paper uses backward Euler for its simplicity and low-rank preservation.
Computing Greeks from the QTT Grid
A full-grid solver gives you more than just the option price at a single point. Once you have the entire solution surface in QTT format, you can compute the Greeks — Delta and Gamma — over the whole grid with negligible additional cost. The key insight is that the derivative operators in log-price space have exact, low-rank QTT representations. Applying them to the solution MPS produces the Greeks in log-price coordinates, and a simple chain-rule conversion gives the familiar Delta and Gamma in the original asset-price coordinates.
This section implements the Greek computation as described in Appendix A.2 of the paper. The generated file greeks.py provides two public functions: compute_delta and compute_gamma. Each can return either a QTTTensor representing the Greek over the entire grid or a single interpolated value at a given spot price.
Derivative MPOs from the Tridiagonal Building Block
The first and second derivatives in log-price space are discretized using central differences. The first derivative uses the stencil [-1/(2dx), 0, 1/(2dx)], and the second derivative uses [1/dx^2, -2/dx^2, 1/dx^2]. Both are tridiagonal Toeplitz matrices, so they can be built exactly using the same qtt_toeplitz_tridiagonal function from qtt_analytic.py (Lemma 1, rank 3).
def _derivative_mpo_first(c, dx):
"""
Build a 1D first-derivative MPO using central differences.
Implements the tridiagonal Toeplitz matrix with stencil [-1/(2dx), 0, 1/(2dx)]
(Lemma 1 / Appendix C.3).
Parameters
----------
c : int
Number of QTT cores (grid size = 2^c).
dx : float
Spatial step size in the log-price domain.
Returns
-------
QTTTensor
MPO of shape (2^c, 2^c) with rank 3.
"""
# Central difference: subdiagonal = -1/(2dx), diagonal = 0, superdiagonal = 1/(2dx)
beta = 0.5 / dx
gamma = -0.5 / dx
return qtt_toeplitz_tridiagonal(c, alpha=0.0, beta=beta, gamma=gamma)
def _derivative_mpo_second(c, dx):
"""
Build a 1D second-derivative MPO using central differences.
Stencil [1/dx^2, -2/dx^2, 1/dx^2] (Lemma 1 / Appendix C.3).
Parameters
----------
c : int
Number of QTT cores (grid size = 2^c).
dx : float
Spatial step size in the log-price domain.
Returns
-------
QTTTensor
MPO of shape (2^c, 2^c) with rank 3.
"""
alpha = -2.0 / (dx * dx)
beta = 1.0 / (dx * dx)
gamma = 1.0 / (dx * dx)
return qtt_toeplitz_tridiagonal(c, alpha=alpha, beta=beta, gamma=gamma)Notice that the first-derivative MPO has a zero diagonal (alpha=0.0), while the second-derivative MPO has alpha = -2/dx^2. Both have rank 3, matching the tridiagonal Toeplitz construction from Lemma 1.
Building a d-Dimensional Derivative MPO
To differentiate with respect to a single asset in a d-asset problem, we need an MPO that applies the 1D derivative in the target dimension and the identity in all other dimensions. The function _build_derivative_mpo_d constructs this by concatenating the cores of the 1D derivative MPO for the target dimension and the cores of the 1D identity MPO for every other dimension.
def _build_derivative_mpo_d(d, c, asset_index, dx):
"""
Build a d-dimensional MPO that applies the first derivative in the
`asset_index`-th dimension and identity in all other dimensions.
The MPO is constructed by concatenating the 1D cores for each dimension:
- For the target dimension: use the 1D derivative MPO cores.
- For all other dimensions: use the 1D identity MPO cores.
Parameters
----------
d : int
Number of assets (dimensions).
c : int
Number of QTT cores per dimension.
asset_index : int
Index of the dimension along which to differentiate (0-based).
dx : float
Spatial step size in the log-price domain.
Returns
-------
QTTTensor
MPO of shape (2^{c*d}, 2^{c*d}).
"""
derivative_1d = _derivative_mpo_first(c, dx)
identity_1d = qtt_identity_mpo(c)
cores = []
for dim in range(d):
if dim == asset_index:
cores.extend(derivative_1d.cores)
else:
cores.extend(identity_1d.cores)
return QTTTensor(cores, is_mpo=True)This function is used for the first derivative. For the second derivative, the same pattern applies but with _derivative_mpo_second in place of _derivative_mpo_first. The resulting MPO has rank 3 in the target dimension and rank 1 in all others, so the overall rank is bounded by 3.
Computing Delta and Gamma
With the derivative MPOs in hand, computing the Greeks is a matter of applying them to the solution MPS and converting from log-price coordinates to asset-price coordinates.
Delta in log-price coordinates is ∂V/∂x. The chain rule gives the familiar Delta in asset-price coordinates:
Gamma requires both the first and second derivatives in log-price:
These formulas follow from the change of variables x = ln S and the chain rule. The paper derives them in Appendix A.2.
Here is the compute_delta function from greeks.py:
def compute_delta(solution_qtt, asset_index, market_params, grid_params, spot=None):
"""
Compute the Delta (∂V/∂S) for a given asset from the QTT solution.
The chain rule is used: Δ = (1/S) * ∂V/∂x, where x = log S.
If `spot` is provided, returns the interpolated Delta at that spot price
for the specified asset, with other assets fixed at their initial spots.
Otherwise, returns a QTTTensor representing the Delta grid over all assets.
Parameters
----------
solution_qtt : QTTTensor
MPS of the option price at the final time (shape 2^{c*d}).
asset_index : int
Index of the asset for which to compute Delta (0-based).
market_params : dict
Must contain 'S0' (list of initial spot prices) and other parameters.
grid_params : dict
Must contain 'c' (int) and 'domain' (tuple (xmin, xmax)).
spot : float, optional
If provided, return the interpolated Delta at this spot.
Returns
-------
QTTTensor or float
The Delta grid (QTTTensor) or interpolated value (float).
"""
c = grid_params['c']
xmin, xmax = grid_params['domain']
dx = (xmax - xmin) / (2**c - 1)
d = len(market_params['S0'])
# Build derivative MPO for the specified asset dimension
derivative_mpo = _build_derivative_mpo_d(d, c, asset_index, dx)
# Apply to get ∂V/∂x
dV_dx = qtt_contract_mpo_mps(derivative_mpo, solution_qtt)
# Round to control rank growth (optional, but good practice)
dV_dx = qtt_round(dV_dx, tol=1e-8)
if spot is not None:
# Interpolate at the given spot
val = _interpolate_at_spot(dV_dx, spot, asset_index, market_params, grid_params)
# Δ = (1/S) * dV/dx
return val / spot
else:
# Build the d-dimensional QTT for e^{-x} on the asset_index dimension
# and all-ones on other dimensions.
# First, build 1D exponential QTT: e^{-x} on the whole domain.
exp_minus_x = analytic_qtt_exponential(alpha=-1.0, c=c, interval=(xmin, xmax))
# Build list of 1D QTTs for Kronecker product
ones_qtt = analytic_qtt_exponential(alpha=0.0, c=c, interval=(xmin, xmax)) # all ones
qtt_list = []
for i in range(d):
if i == asset_index:
qtt_list.append(exp_minus_x)
else:
qtt_list.append(ones_qtt)
# Use qtt_kronecker_product to combine (only works for rank-1 MPS)
from qtt_core import qtt_kronecker_product
inv_s_qtt = qtt_kronecker_product(qtt_list)
# Element-wise multiply: Δ = dV_dx * e^{-x}
delta_qtt = _elementwise_mul_mps(dV_dx, inv_s_qtt)
# Round to control rank
delta_qtt = qtt_round(delta_qtt, tol=1e-8)
return delta_qttThe function first builds the derivative MPO for the specified asset dimension and applies it to the solution to get ∂V/∂x. If a specific spot price is requested, it interpolates ∂V/∂x at that point and divides by spot to get Delta. Otherwise, it builds a QTT representation of 1/S = e^{-x} using the analytic exponential function (rank 1) and multiplies element-wise with ∂V/∂x to get the full Delta grid.
The compute_gamma function follows the same pattern but computes both ∂V/∂x and ∂²V/∂x², subtracts them, and multiplies by e^{-2x}:
def compute_gamma(solution_qtt, asset_index, market_params, grid_params, spot=None):
"""
Compute the Gamma (∂²V/∂S²) for a given asset from the QTT solution.
The chain rule is used: Γ = (1/S²) * (∂²V/∂x² - ∂V/∂x), where x = log S.
If `spot` is provided, returns the interpolated Gamma at that spot price
for the specified asset. Otherwise, returns a QTTTensor representing the
Gamma grid over all assets.
Parameters
----------
solution_qtt : QTTTensor
MPS of the option price at the final time (shape 2^{c*d}).
asset_index : int
Index of the asset for which to compute Gamma (0-based).
market_params : dict
Must contain 'S0' (list of initial spot prices).
grid_params : dict
Must contain 'c' (int) and 'domain' (tuple (xmin, xmax)).
spot : float, optional
If provided, return the interpolated Gamma at this spot.
Returns
-------
QTTTensor or float
The Gamma grid (QTTTensor) or interpolated value (float).
"""
c = grid_params['c']
xmin, xmax = grid_params['domain']
dx = (xmax - xmin) / (2**c - 1)
d = len(market_params['S0'])
# ---- Compute ∂V/∂x ----
deriv1_mpo = _build_derivative_mpo_d(d, c, asset_index, dx)
dV_dx = qtt_contract_mpo_mps(deriv1_mpo, solution_qtt)
dV_dx = qtt_round(dV_dx, tol=1e-8)
# ---- Compute ∂²V/∂x² ----
# Build second derivative MPO for the target dimension
# We need a d-dimensional MPO for the second derivative.
# Reuse the same construction as first derivative but with second derivative MPO.
deriv2_1d = _derivative_mpo_second(c, dx)
identity_1d = qtt_identity_mpo(c)
cores = []
for dim in range(d):
if dim == asset_index:
cores.extend(deriv2_1d.cores)
else:
cores.extend(identity_1d.cores)
deriv2_mpo = QTTTensor(cores, is_mpo=True)
d2V_dx2 = qtt_contract_mpo_mps(deriv2_mpo, solution_qtt)
d2V_dx2 = qtt_round(d2V_dx2, tol=1e-8)
# ---- Compute (∂²V/∂x² - ∂V/∂x) ----
# Element-wise subtraction
# We need a subtraction function; we can use qtt_add with negative scalar
from qtt_core import qtt_add, qtt_mul_scalar
neg_dV_dx = qtt_mul_scalar(dV_dx, -1.0)
diff = qtt_add(d2V_dx2, neg_dV_dx)
diff = qtt_round(diff, tol=1e-8)
if spot is not None:
# Interpolate the difference at the spot
val = _interpolate_at_spot(diff, spot, asset_index, market_params, grid_params)
# Γ = (1/S²) * val
return val / (spot * spot)
else:
# Build the d-dimensional QTT for e^{-2x} on the asset_index dimension
exp_minus_2x = analytic_qtt_exponential(alpha=-2.0, c=c, interval=(xmin, xmax))
ones_qtt = analytic_qtt_exponential(alpha=0.0, c=c, interval=(xmin, xmax))
qtt_list = []
for i in range(d):
if i == asset_index:
qtt_list.append(exp_minus_2x)
else:
qtt_list.append(ones_qtt)
from qtt_core import qtt_kronecker_product
inv_s2_qtt = qtt_kronecker_product(qtt_list)
# Element-wise multiply: Γ = diff * e^{-2x}
gamma_qtt = _elementwise_mul_mps(diff, inv_s2_qtt)
gamma_qtt = qtt_round(gamma_qtt, tol=1e-8)
return gamma_qttWorked Example: Delta and Gamma for a 1-Asset European Call
Let us walk through a concrete example. Suppose we have solved the Black–Scholes PDE for a 1-asset European call with c = 8 cores (256 grid points), sigma = 0.25, r = 0.05, T = 1.0, and K = 10. The solution MPS V has shape (256,) and rank typically around 5–10.
To compute Delta at the spot price S0 = 10:
Build the 1D first-derivative MPO
D_xwithdx = (xmax - xmin) / (2^c - 1).Apply
D_xtoVto getdV_dx(a new MPS of shape(256,)).Interpolate
dV_dxatx0 = ln(10)to getdV_dx_at_spot.Divide by
S0 = 10to getDelta = dV_dx_at_spot / 10.
For Gamma, we additionally compute d2V_dx2 by applying the second-derivative MPO, then compute diff = d2V_dx2 - dV_dx, interpolate at x0, and divide by S0^2 = 100.
The paper reports in Table IV that for a 1-asset European call, the grid-wide Delta and Gamma errors are on the order of 1e-4 to 1e-3, and the computation adds only a few percent to the total solver runtime. This is because the derivative MPOs have rank at most 3, and the element-wise multiplication with e^{-x} or e^{-2x} uses rank-1 MPS, so the overall rank of the Greek tensors remains small.
Advanced Detail: Interpolation on the Dyadic Grid
The function _interpolate_at_spot uses nearest-neighbor interpolation on the dyadic log-price grid. For a given spot price S, it computes x = ln(S), clamps it to the domain [xmin, xmax], finds the nearest grid index, and evaluates the MPS at that multi-index using _evaluate_mps_at_indices. The evaluation contracts the MPS cores with one-hot vectors corresponding to the binary representation of each index, which is efficient because the MPS is never expanded to its full vector form.
For higher accuracy, one could replace nearest-neighbor with linear or cubic interpolation, but the paper shows that nearest-neighbor is sufficient for the reported error levels.
Summary
Derivative MPOs for first and second derivatives in log-price space are built from the analytic tridiagonal Toeplitz representation (Lemma 1, rank 3).
A d-dimensional derivative MPO is constructed by concatenating 1D derivative cores for the target dimension and 1D identity cores for all other dimensions.
Delta and Gamma are computed by applying the derivative MPOs to the solution MPS and converting from log-price to asset-price coordinates using the chain rule.
The conversion uses the analytic rank-1 QTT for
e^{-x}(for Delta) ore^{-2x}(for Gamma), multiplied element-wise with the derivative MPS.The entire computation is efficient because all operations stay within the QTT format, and the ranks of the derivative MPOs and conversion factors are small (3 or 1).
Numerical Experiments and Results
We now bring together all the components built in the previous sections to reproduce the key numerical results from the paper. The experiments cover European and American basket puts and max-min puts in dimensions d = 3, 4, 5, using the time-stepping solver as the primary tool. For European options we also test the space-time solver and confirm that it delivers the same prices (within truncation error) with a different efficiency profile. All results are compared against reference prices obtained by independent methods: Gauss–Hermite quadrature for European contracts and Longstaff–Schwartz Monte Carlo for American contracts.
Experimental Setup
The market parameters are taken from Appendix B.2 of the paper. A single set of parameters is used across all experiments, with the number of assets d determining which entries are active:
Spot prices (log‑price coordinates):
S = (10, 11, 12, 13, 14)for assets 1 through 5.Volatilities:
sigma = (0.25, 0.15, 0.20, 0.10, 0.15).Correlation matrix
rhoof size 5 x 5 with off‑diagonal entries(0.4, 0.3, 0.2, 0.1, 0.4, ...)(see paper for the full matrix).Risk‑free rate
r = 0.05(5% per annum).Maturity
T = 1.0year.Strikes: for basket puts,
K = 34, 47, 62for d = 3, 4, 5 respectively; for max‑min puts,K = 10for all d.
The grid parameters are:
Number of spatial cores
cper dimension, typically 7, 8, or 9. The grid per dimension therefore hasN = 2^cpoints.Number of time steps in the time‑stepping solver is also
2^c(the paper uses a fully implicit scheme withtheta = 1).The spatial domain
[x_min, x_max]^dis determined adaptively by a pilot coarse run (see the functionadaptive_domaininutils.py).
Reference prices are computed as follows:
European basket put: Gauss–Hermite quadrature with 100 points per dimension (using the function
gauss_hermite_quadrature_basketinutils.py). This gives a highly accurate reference value.European max‑min put: The same Gauss–Hermite method applied to the payoff
max(K - (min_i e^{x_i}), 0).American basket and max‑min puts: Longstaff–Schwartz least‑squares Monte Carlo with 100,000 paths and 100 time steps (function
longstaff_schwartz_american). The Monte Carlo reference itself has a standard deviation of roughly 1–2% of the price.
Error metrics used in all tables:
|price_error| = |price_QTT - price_ref|(absolute error).relative_error = |price_error| / price_ref.For Greeks, the same metrics are computed over the whole grid (Delta and Gamma are compared to the reference grid, which is obtained by finite differences on the Gauss–Hermite solution).
Running a Basket Experiment
The central experiment for a given d and option type is orchestrated by the function run_basket_experiment in experiments.py. A simplified version of its body is shown below; it relies on the solver and utilities built earlier.
def run_basket_experiment(d, c, option_type, ref_price):
# Market parameters for d assets
params = default_market_params(d) # from utils.py
# Adaptive grid domain (pilot coarse run with c=5)
domain = adaptive_domain(params, c_coarse=5)
grid_params = {'c': c, 'domain': domain}
# Solve using the time-stepping solver
solution_qtt = time_stepping_solver(
params, grid_params,
option_type=option_type,
contract_type='basket'
)
# Interpolate at the spot price (first element of params['spot'])
spot = params['spot'][0] # for d=3, S0 = 10
price = solution_qtt.interpolate_at_spot(spot)
# Compute error
error = abs(price - ref_price) / ref_price
return price, errorFor American options the same function is called with option_type='American'; the solver internally applies the early‑exercise condition (via tt_cross_max) after every time step.
Reproducing the Paper’s Tables
The paper presents its main results in Tables V through VIII. The following paragraphs summarise those results and explain how the code reproduces them. Because no computation was executed during this writing, the values quoted here are taken directly from the paper; running the code with the same parameters should yield numbers matching the published tables within the reported tolerances.
Table V: European basket put (d = 3, 4, 5). The time‑stepping solver with c = 8 cores per dimension (256 grid points per asset) achieves absolute pricing errors below 0.05 for all dimensions. The run time grows modestly with d: for d = 3 the solver completes in under 10 seconds on the reference hardware (Apple M3 Pro); for d = 5 the runtime is still under one minute. The rank of the solution never exceeds 15. The space‑time solver, when applied to the same problems, yields identical prices (within 1e‑6) and is faster for d = 3 but requires more memory because the global solution surface has rank that is 2–3 times larger than the time‑stepping representation.
Table VI: European max‑min put. The max‑min payoff introduces the min function, which is approximated by TT‑Cross with rank cap chi = 12. The resulting pricing errors are slightly larger than for the basket contract — typically 0.1–0.3 absolute error — but still well within the 1–2% range needed for practical applications. The run times are comparable to the basket case because the extra TT‑Cross cost appears only during the initial payoff construction.
Table VII: American basket put. The early‑exercise condition adds about 30% to the total runtime. The paper reports that the rank of the solution grows by at most 2–3 additional bond dimensions compared to the European case, confirming that the QTT representation remains efficient. Pricing errors are slightly larger than for European options because of the additional approximation introduced by the TT‑Cross max operation at each time step.
Table VIII: American max‑min put. This is the most challenging contract because both the payoff and the early‑exercise condition involve the min and max functions. The observed ranks are 2–3 times higher than for the basket put, and the runtime increases accordingly. Nevertheless, for d = 3 the solver still finishes in under 30 seconds, and for d = 5 it remains under 5 minutes on a personal computer.
Worked Example: d = 3 Basket Put, c = 8
To illustrate the typical workflow, consider the 3‑asset European basket put. The steps are:
Load the market parameters:
sigma = [0.25, 0.15, 0.20],rho(the 3x3 block of the full matrix),r = 0.05,K = 34,T = 1.0.Determine the domain via a coarse pilot run (e.g.,
c = 5, domain expands to cover the region where the payoff is non‑zero). The result might be[x_min, x_max] = [-2.1, 2.1].Build the operator MPO
A(using the analytic QTT building blocks) and the payoff QTT (usingpayoff_basket_put).Run the time‑stepping loop: 256 time steps, each solving
(I - dt A) x_{n+1} = x_nwith the ALS solver (2 sweeps, initial guess from previous step).Extract the price at the spot point
(10, 11, 12)by interpolating the final time‑slice QTT.
The computed price should match the reference (obtained by Gauss–Hermite quadrature) to within 0.02–0.03 absolute error, corresponding to a relative error of less than 1%. The entire computation takes about 8 seconds on the reference machine.
Advanced Details
Rank behaviour: For all experiments reported in the paper, the bond dimensions of the solution QTT remained below 20 even for the most demanding 5‑asset American max‑min case. This is crucial for the polynomial scaling of the method.
Space‑time vs. time‑stepping: The space‑time solver is approximately twice as fast as the time‑stepping solver for d = 3 European options, but its memory footprint is larger because it stores the entire space‑time surface. For d > 3 the time‑stepping solver is preferred because the space‑time operator rank grows more rapidly.
Greek computation: The functions
compute_deltaandcompute_gamma(fromgreeks.py) can be called on the final solution QTT to obtain full‑grid Greeks with negligible extra cost. The paper reports that the grid‑wide Delta and Gamma errors are comparable to the price errors (Table IV for the 1D case; similar behaviour holds for d > 1).
Code Organisation
The experimental workflow is encapsulated in experiments.py. Its two main entry points are:
def run_basket_experiment(d, c, option_type, ref_price=None):
"""
Run the time-stepping solver for a basket put.
Returns the QTT solution, the interpolated spot price, and the relative error.
If ref_price is given, error is computed; otherwise it returns None.
"""
...
def run_maxmin_experiment(d, c, option_type, ref_price=None):
"""
Similar for max-min put.
"""
...These functions use the utility functions default_market_params, adaptive_domain, and the reference‑price calculators from utils.py. The paper’s tables can be reproduced by looping over d, c, and option_type and recording the outputs.
No code execution is performed as part of this tutorial; the reader is encouraged to run the solver and compare the results with the published Tables V–VIII.
Use the Button below to download the source code.