Tensor Cross Interpolation

Tensor Cross Interpolation (TCI) approximates a high-dimensional function as a tensor train by sampling only a small fraction of its entries. This crate provides two levels of API:

• tensor4all-tensorci for low-level TCI directly on integer indices.
• tensor4all-quanticstci for quantics interpolation on discrete or continuous grids.

Low-Level TCI (tensor4all-tensorci)

Use crossinterpolate2 when you already know the local dimensions and want direct control over the algorithm.

Defining the function

The function must accept a &Vec<usize> of 0-indexed multi-indices and return a scalar value. Here is a simple 2D example where f(i, j) = i + j + 1:

#![allow(unused)]
fn main() {
use tensor4all_simplett::prelude::*;
use tensor4all_tensorci::prelude::*;

let f = |idx: &Vec<usize>| (idx[0] + idx[1] + 1) as f64;
let local_dims = vec![4, 4];
let initial_pivots = vec![vec![3, 3]]; // pick where |f| is large

let result = crossinterpolate2::<f64, _, fn(&[Vec<usize>]) -> Vec<f64>>(
    f,
    None,
    local_dims,
    initial_pivots,
    TCI2Options {
        tolerance: 1e-10,
        seed: Some(42),
        ..Default::default()
    },
).unwrap();

assert_eq!(result.termination, TCI2Termination::Converged);

let tt = result.tci.to_tensor_train().unwrap();
let value = tt.evaluate(&[2, 3]).unwrap();
assert!((value - 6.0).abs() < 1e-10);
assert!(tt.rank() >= 1);
}

Choosing TCI2Options

The most important parameters:

Interpreting the results

crossinterpolate2 returns a TCI2OptimizationResult:

Do not infer convergence from errors.last() alone. Full convergence also requires no recent global pivots and a stable rank history.

The default global pivot finder walks the residual |f - tt| from random starts, retaining each maximizing coordinate for the next coordinate scan. It repeats full sweeps until the residual stops increasing, exceeds ten times the acceptance threshold, or reaches 100 sweeps. The acceptance threshold is abs_tol * tol_margin_global_search, where abs_tol includes the configured error normalization. Tensor-train evaluations use a shared cache and batches of coordinate candidates.

Custom GlobalPivotFinder implementations return Result<Vec<MultiIndex>> and require T: Scalar + TTScalar + EinsumScalar; use Ok(pivots) for a successful search. optimize_with_finder propagates search failures. The EinsumScalar bound also applies to crossinterpolate2, integrate, and the adaptive interpolation entry points that call them.

Convergence diagnostics

The errors vector tracks the normalized bond error after each half-sweep. The algorithm converges when:

1. The last ncheck_history (default: 3) entries are all below tolerance.
2. No global pivots were added in those iterations.
3. The rank has stabilized.

If the errors plateau above your tolerance, try:

• Increasing max_bond_dim (the function may need higher rank).
• Increasing max_iter (more sweeps may be needed).
• Choosing better initial pivots (where |f| is large).

Convert to a tensor train for further manipulation:

#![allow(unused)]
fn main() {
use tensor4all_simplett::prelude::*;
use tensor4all_tensorci::prelude::*;
let f = |idx: &Vec<usize>| (idx[0] + idx[1] + 1) as f64;
let tensor4all_tensorci::TCI2OptimizationResult { tci, ranks: _ranks, errors: _errors, .. } = crossinterpolate2::<f64, _, fn(&[Vec<usize>]) -> Vec<f64>>(
    f, None, vec![4, 4], vec![vec![3, 3]],
    TCI2Options { seed: Some(42), ..Default::default() },
).unwrap();
let tt = tci.to_tensor_train().unwrap();
assert!(tt.rank() >= 1);

// Evaluate the tensor train at specific indices
let val = tt.evaluate(&[1, 2]).unwrap();
assert!((val - 4.0).abs() < 1e-10); // f(1,2) = 1+2+1 = 4
}

Continuous vs discrete

crossinterpolate2 works on discrete integer indices. For functions on continuous domains, use the quantics representation provided by tensor4all-quanticstci (see below), which maps floating-point coordinates to binary tensor-train indices.

High-Level Quantics TCI (tensor4all-quanticstci)

The quantics representation encodes each grid index in binary and arranges the bits across tensor-train sites. This often yields much lower bond dimensions than a naive encoding.

Important conventions

• Indexing differs between the two APIs:
  • crossinterpolate2 (low-level): indices and pivots are 0-indexed (0..local_dim)
  • quanticscrossinterpolate_discrete_batch (high-level): grid indices are 0-indexed (0..grid_size)
• Equal dimensions: quanticscrossinterpolate_discrete_batch requires all dimensions to have the same number of points.
• Power-of-2 grid sizes: all grid dimensions must be powers of 2 (4, 8, 16, 32, …).

Choosing between discrete and continuous APIs

Tip on n_random_init_pivot: The n_random_init_pivot option (default: 5) controls how many random initial pivot points are used to seed the TCI algorithm. For functions with multiple separated features or high-dimensional problems, increase this to 10–20 to improve robustness.

Discrete grid interpolation

Use quanticscrossinterpolate_discrete_batch when your function is naturally defined on an integer grid. Indices are passed as &[usize] and are 0-indexed.

fn main() -> anyhow::Result<()> {
use tensor4all_quanticstci::prelude::*;

let f = |idx: &[usize]| (idx[0] + idx[1]) as f64;
let sizes = vec![16, 16];

let (qtci, ranks, errors) = quanticscrossinterpolate_discrete_batch::<f64, _>(
    &sizes, pointwise_index_batch(f),
    None,
    QtciOptions::default().with_tolerance(1e-10))?;

assert!(*errors.last().unwrap() < 1e-10);
assert!(!ranks.is_empty());

let value = qtci.evaluate(&[5, 10])?;
assert!((value - 15.0).abs() < 1e-8);

// Sum of (i + j) for i, j in 0..16 = 2 * 16 * (15 * 16 / 2) = 3840
let total = qtci.sum()?;
assert!((total - 3840.0).abs() < 1e-6);
Ok(())
}

Continuous grid interpolation with DiscretizedGrid

For functions on continuous domains, build a DiscretizedGrid that maps grid indices to physical coordinates. The number of quantics bits per dimension is set via the builder.

fn main() -> anyhow::Result<()> {
use tensor4all_quanticstci::prelude::*;

// 2^4 = 16 grid points on [0, 1)
let grid = DiscretizedGrid::builder(&[4])
    .with_lower_bound(&[0.0])
    .with_upper_bound(&[1.0])
    .build()
    .unwrap();

let f = |x: &[f64]| x[0] * x[0];

let (qtci, _ranks, errors) = quanticscrossinterpolate_batch::<f64, _>(
    &grid, pointwise_coordinate_batch(f),
    None,
    QtciOptions::default())?;

assert!(*errors.last().unwrap() < 1e-8);

// Verify the interpolation at grid point 0 (x = 0.0)
let value = qtci.evaluate(&[0])?;
assert!((value - 0.0).abs() < 1e-10);
Ok(())
}

Integration

integral() computes a left Riemann sum:

integral ≈ Σ f(xᵢ) × Δx

This has O(h) convergence where h is the grid spacing. The result depends on whether include_endpoint is set on the DiscretizedGrid.

fn main() -> anyhow::Result<()> {
use tensor4all_quanticstci::prelude::*;
let grid = DiscretizedGrid::builder(&[4])
    .with_lower_bound(&[0.0])
    .with_upper_bound(&[1.0])
    .build()
    .unwrap();
let f = |x: &[f64]| x[0] * x[0];
let (qtci, _, _) = quanticscrossinterpolate_batch::<f64, _>(
    &grid, pointwise_coordinate_batch(f), None, QtciOptions::default(),
)?;
let integral = qtci.integral()?;
// Left Riemann sum of x^2 over [0, 1) with 16 points
assert!((integral - 1.0 / 3.0).abs() < 5e-2);

let sum = qtci.sum()?;
assert!((sum * grid.grid_step()[0] - integral).abs() < 1e-12);
Ok(())
}

For discrete grids created without a continuous domain, integral() returns the plain sum identical to sum().

Practical Example: Multi-scale 1D Function

This section corresponds to the Julia Quantics TCI of univariate function notebook.

The function below mixes several length scales:

f(x) = cos(x/B) * cos(x/(4*sqrt(5)*B)) * exp(-x^2) + 2*exp(-x)

with B = 2^-30 on [0, ln(20)) using R = 40 quantics bits.

fn main() -> anyhow::Result<()> {
use tensor4all_quanticstci::prelude::*;

// Use R = 10 bits (2^10 = 1024 points) for a fast doctest.
// In practice, set R = 40 for the full multi-scale resolution.
let r = 10;
let x_max = 20.0_f64.ln();
let grid = DiscretizedGrid::builder(&[r])
    .with_lower_bound(&[0.0])
    .with_upper_bound(&[x_max])
    .include_endpoint(false)
    .build()
    .unwrap();

// A smooth multi-scale function: f(x) = cos(10x) * exp(-x^2) + 2*exp(-x)
let f = move |coords: &[f64]| {
    let x = coords[0];
    (10.0 * x).cos() * (-x * x).exp() + 2.0 * (-x).exp()
};

let tol = 1e-8;
let (qtci, _ranks, errors) = quanticscrossinterpolate_batch::<f64, _>(
    &grid, pointwise_coordinate_batch(f),
    None,
    QtciOptions::default()
        .with_tolerance(tol)
        .with_max_bond_dim(64)
        .with_nrandominitpivot(8))?;

assert!(*errors.last().unwrap() < tol);

// Verify at several grid points across the domain
for &grid_idx in &[0usize, 99, 511, 1023] {
    let got = qtci.evaluate(&[grid_idx])?;
    assert!(got.is_finite(), "value at grid index {} should be finite", grid_idx);
}

// Check the integral is reasonable (positive, since f > 0 near x = 0)
let integral = qtci.integral()?;
assert!(integral > 1.0);
Ok(())
}

Multivariate (2D) Example

This section corresponds to the Julia Quantics TCI of multivariate function notebook.

fn main() -> anyhow::Result<()> {
use tensor4all_quanticstci::prelude::*;

// Use 64x64 for a faster doctest (256x256 in practice)
let sizes = vec![64, 64];
let f = |idx: &[usize]| {
    let x = idx[0] as f64;
    let y = idx[1] as f64;
    (x / 24.0).cos() + (y / 17.0).cos() + 0.1 * ((x + y) / 13.0).sin()
};

let (qtci, _ranks, errors) = quanticscrossinterpolate_discrete_batch::<f64, _>(
    &sizes, pointwise_index_batch(f),
    None,
    QtciOptions::default()
        .with_tolerance(1e-10)
        .with_max_bond_dim(64)
        .with_nrandominitpivot(8))?;

assert!(*errors.last().unwrap() < 1e-8);

for &(i, j) in &[(0usize, 0usize), (16, 32), (31, 49), (63, 63)] {
    let exact = (i as f64 / 24.0).cos()
        + (j as f64 / 17.0).cos()
        + 0.1 * (((i + j) as f64) / 13.0).sin();
    let got = qtci.evaluate(&[i, j])?;
    assert!((got - exact).abs() < 1e-6);
}
Ok(())
}