# Chebyshev approximation and composition of functions in matrix product
states for quantum-inspired numerical analysis

Juan José Rodríguez-Aldavero  Institute of Fundamental Physics IFF-CSIC, Calle
Serrano 113b, Madrid 28006, Spain  Paula García-Molina  Institute of
Fundamental Physics IFF-CSIC, Calle Serrano 113b, Madrid 28006, Spain  Luca
Tagliacozzo  Institute of Fundamental Physics IFF-CSIC, Calle Serrano 113b,
Madrid 28006, Spain  Juan José García-Ripoll  Institute of Fundamental Physics
IFF-CSIC, Calle Serrano 113b, Madrid 28006, Spain

(Jul 2024)

###### Abstract

This work explores the representation of univariate and multivariate functions
as matrix product states (MPS), also known as quantized tensor-trains (QTT).
It proposes an algorithm that employs iterative Chebyshev expansions and
Clenshaw evaluations to represent analytic and highly differentiable functions
as MPS Chebyshev interpolants. It demonstrates rapid convergence for highly-
differentiable functions, aligning with theoretical predictions, and
generalizes efficiently to multidimensional scenarios. The performance of the
algorithm is compared with that of tensor cross-interpolation (TCI) and
multiscale interpolative constructions through a comprehensive comparative
study. When function evaluation is inexpensive or when the function is not
analytical, TCI is generally more efficient for function loading. However, the
proposed method shows competitive performance, outperforming TCI in certain
multivariate scenarios. Moreover, it shows advantageous scaling rates and
generalizes to a wider range of tasks by providing a framework for function
composition in MPS, which is useful for non-linear problems and many-body
statistical physics.

##  I Introduction

Tensor networks [1] are decompositions of tensors in networks of mutually
contracted smaller tensors. On some occasions, these decompositions can result
in accurate, exponentially more efficient approximations. Matrix product
states (MPS), originally discovered in the study of strongly correlated
quantum systems and quantum Markov chains [2, 3], are a succesful example of
these decompositions. MPS underlie the structure of the DMRG algorithm [4, 5],
a well-known and powerful method for diagonalizing many-body Hamiltonians [6].
They were rediscovered under the name of tensor-trains (TT) [7] and quantized
tensor-trains (QTT) [8] as an efficient decomposition of high-dimensional
tensors, featuring approximation tools similar to DMRG. Currently, MPS stand
as a platform for quantum-inspired numerical methods [9], with applications in
Fourier analysis and interpolation [10, 9, 11, 12], extrapolation [13],
solving differential equations [14, 15, 16, 17] and multivariate integration
[18, 19, 20], among others. These methods potentially yield exponential
advantages in time and memory usage over conventional numerical techniques
[15, 20, 16]. Furthermore, the interpretation of tensor network formalisms as
quantum circuits [21] makes them an efficient method for discovering quantum
algorithms to load specific functions in quantum registers [22], thus
interfacing with quantum computing applications.

A central problem in using tensor-network decompositions for numerical
analysis is loading generic functions into such representations. The
conventional approach is based on discretizing such functions over a
multidimensional grid and using singular value decompositions (SVD). However,
this method scales polynomially in time and memory with the volume of the
grid, thus exponentially with its dimension, becoming susceptible to the so-
called “curse of dimensionality”. There are essentially two types of methods
that tackle this curse. Some methods rely on a clever sampling of the function
to gather an optimal tensor network structure. This family of methods includes
multiple variations of tensor cross-interpolation (TCI) [23] and generative
modeling techniques [24]. Other approaches rely on algebraic representations
of the function that enable constructive approximations of the optimal tensor
structure. These include methods based on polynomial interpolation [25, 26]
and algorithmic constructs that map concrete algebraic representations of
functions to tensor networks [19].

In this work, we propose a combination of the latter two approaches, relying
on Chebyshev interpolation [27, 28] to create polynomial approximations that
can be efficiently constructed as MPS, using a finite-precision non-linear
algebra framework [9, 16]. This extends the algebraic circuits from Ref. [19]
to more generic function representations. The resulting algorithm, heavily
inspired by previous uses of quantum filters [29], provides approximations
that converge exponentially fast with the circuit depth for analytic
functions. Furthermore, this method not only allows us to express a generic
function $f(x)$ but also to combine approximations, composing functions
$f(g(x))$ for which a Chebyshev interpolant has been developed, as well as
their derivatives and primitives, providing constructive approximations for
higher-dimensional composite functions—e.g., $f(g(x)h(y))$,
$f(g^{\prime\prime}(x)\int h(y))$, etc.

The performance of the proposed algorithm is tested on a collection of
univariate and multivariate functions. In the univariate scenario, the method
demonstrates a stable convergence that matches the theoretical convergence
rates. At the same time, it exhibits an efficient performance, requriring
polynomial resources and avoiding the curse of dimensionality. In particular,
its runtime scales close to linearly with the approximation order and sub-
linearly with the number of qubits, hence sub-logarithmically with the density
of the discretization. Moreover, the method generalizes to the multivariate
scenario maintaining its accurate convergence rates and polynomial resource
requirements.

Despite this, in the univariate scenario the algorithm generally does not
outperform tensor cross-interpolation or multiscale interpolative
constructions. This is particularly true when function evaluations are
inexpensive or when the function is not highly-differentiable. However, it
outperforms TCI in some multivariate scenarios, such as those that benefit
from the interleaved qubit order. Moreover, the method generalizes for the
composition of functions on generic MPS/QTT structures, a task where the
remaining algorithms are less applicable. This capability is useful in
scenarios such as non-linear problems [30] and many-body statistical physics
[31]. In addition, the method can applied in a straightforward manner to
compose functions on generic matrix product operators (MPO).

This work is organized as follows: Section II presents the mapping between
functions and MPS and the finite-precision non-linear algebra framework
required to approximate the MPS operations. Section III introduces the
Chebyshev approximation algorithm and its formulation in terms of MPS. Section
IV benchmarks the Chebyshev expansion over a collection of univariate and
multivariate functions, comparing the results with TCI and multiscale
interpolative constructions. Finally, Section V discusses this work’s results
with a focus on applications and future extensions.

##  II Quantum-inspired numerical analysis

Low-dimensional, weakly entangled quantum states admit an efficient
representation as matrix product states [32]. In this form, the exponentially
large tensor of amplitudes of a quantum system with $n$ components

$\ket{\psi}=\sum_{s_{1}=1}^{d_{1}}\sum_{s_{2}=1}^{d_{2}}\cdots\sum_{s_{n}=1}^{d_{n}}\psi_{s_{1}s_{2}\ldots
s_{n}}\ket{s_{1},s_{2},\ldots,s_{n}}$ (1)

becomes a collection of mutually contracted tensors

$\psi_{s_{1}s_{2}\ldots
s_{n}}=\sum_{\alpha_{1}=1}^{\chi_{1}}\sum_{\alpha_{2}=1}^{\chi_{2}}\cdots\sum_{\alpha_{n}=1}^{\chi_{n}}A^{s_{1}}_{\alpha_{1}}A^{s_{2}}_{\alpha_{1}\alpha_{2}}\cdots
A^{s_{n}}_{\alpha_{n}},$ (2)

where a product of matrices determines the individual amplitudes.

Figure 1: Matrix product state composed of $n$ mutually contracted rank-3
tensors $A_{i}$, with virtual indices $\alpha_{i}$ and physical indices
$s_{i}$.

In Penrose’s graphical notation [c.f. Figure 1], this decomposition is
depicted as a one-dimensional train of boxes, where each box correspond to a
tensor, the lines that connect them represent contracted indices,
$\alpha_{i}$, and the free indices become freely standing “legs”, $s_{i}$
where $i=1,\ldots,n$ .

By construction, the number of elements in a MPS grows polynomially with the
number of components $\mathcal{O}(ns\chi_{\max}^{2})$, where
$\chi_{\max}=\max\chi_{i}$ is the largest “bond dimension”. This bond
dimension does not depend on the system size for weakly entangled states. In
those situations, the MPS representation is exponentially more efficient than
the tensor representation $\psi_{s_{1},s_{2}\ldots s_{n}}$ of the quantum
state. Let us now discuss how these improvements apply to the representation
of functions.

###  II.1 MPS/QTT representation of functions

Besides aiding in the study of quantum many-body systems, MPS are useful for
encoding multidimensional data. In the functional tensor train (FTT)
formalism, a multivariate function is encoded as a tensor with one index per
variable, and the result is decomposed into MPS form. This technique was
independently discovered in quantum mechanics [33] and applied mathematics
[34]. In the “quantized” or “quantics” tensor train (QTT) formalism, which is
the focus of this work, each dimension uses multiple indices, enabling greater
flexibility and efficiency. This technique has been rediscovered in multiple
works from mathematics [8, 35] and quantum-inspired algorithms [36, 30, 9,
15].

The QTT transforms a discretized function into a tensor whose indices label
the grid points. For a one-dimensional continuous variable $x$ defined in an
interval $[0,1)$, a set of $n$ digits (or indices) running over $k$ values,
$s_{i}\in\\{0,1,\ldots,k-1\\}$, can label up to $k^{n}$ points. In particular,
when the grid is uniform, a $k$-nary representation of the points is given by

$x_{s_{1}s_{2}\ldots s_{n}}=0.s_{1}s_{2}\ldots
s_{n}=\sum_{i=1}^{n}s_{i}k^{-i}.$ (3)

These $k^{n}$ values may be regarded as a tensor with $n$ indices, $s_{1}$ to
$s_{n}$, each addressing progressively smaller length scales. The tensor
encoding of the $x$ coordinate generalizes to arbitrary univariate functions

$f_{s_{1}s_{2}\ldots s_{n}}:=f(x_{s_{1}s_{2}\ldots s_{n}}),$ (4)

as well as extend to multidimensional functions, with the caveat that we now
have multiple choices for the ordering of indices [9]. In this work we use
either a sequential or serial order

$g_{s_{1}^{1}s_{2}^{1}\ldots s_{n}^{1}\ldots s_{1}^{d}s_{2}^{d}\ldots
s_{n}^{d}}=g(x_{s_{1}^{1}s_{2}^{1}\ldots
s_{n}^{1}},\ldots,x_{s_{1}^{d}s_{2}^{d}\ldots s^{d}_{n}}),$ (5)

and an interleaved order

$g_{s_{1}^{1}\ldots s_{1}^{d},s_{2}^{1}\ldots s_{2}^{d}\ldots s_{n}^{1}\ldots
s_{n}^{d}}=g(x_{s_{1}^{1}s_{2}^{1}\ldots
s_{n}^{1}},\ldots,x_{s_{1}^{d}s_{2}^{d}\ldots s^{d}_{n}}).$ (6)

All these tensors are amenable to an MPS/QTT decomposition (2) with an
efficiency that is determined by their bond dimension structure. Often, their
“compressibility” depends on factors such as the algebraic structure,
smoothness, or displacement structure of the function [37]. Moreover, the bond
dimensions also depend on the value of $k$, which in this work is taken to be
$k=2$, operating with two-dimensional “quantum” components or “qubits”.

Sufficiently smooth univariate functions are known to be highly compressible,
as they yield a MPS/QTT with bounded entanglement entropy [9, 38]. In
particular, polynomials of finite degree $d$ admit an exact representation
with $\chi_{\text{max}}\leq d+1$ [25] which can be constructed analytically as
proven in Appendix A—unfortunately, this construction is ill-conditioned for
constructing Chebyshev interpolants. Interestingly, Chebyshev polynomials
discretized over the Chebyshev-Gauss nodes require a constant bond dimension
$\chi_{\text{max}}=2$ [8]. Other canonical examples include the exponential
function $(\chi_{\text{max}}=1)$, the sine and cosine functions
$(\chi_{\text{max}}=2)$ and periodic vectors of period $\ell$, which require
$\chi_{\text{max}}=\ell$.

In general, loading a function $f(x)$ in an MPS requires a specialized
algorithm. A naïve approach discretizes the function on a tensor
$f_{s_{1}s_{2}\ldots s_{n}}$ and applies a sequence of singular value
decompositions (SVD). This creates an MPS with an optimal error in Frobenius
norm [23] using a time and memory that scale with the volume of the tensor.
Fortunately, other algorithms exist with a logarithmic cost relative to the
size of the discretization. As mentioned in the introduction, these methods
are divided into those that require some randomized or induced sampling—e.g.,
tensor cross-interpolation methods (TCI) [23, 39, 40, 41] for both tensor-
trains and for QTT [40, 42], and sketching methods [24]—and those that are
based on some algebraic structure—such as rank-revealing interpolative
constructions [25] and arithmetic circuit tensor networks [19].

This work proposes a different approximation technique that combines the
benefits of interpolative formulas [25] with the algebraic structures of
arithmetic circuit tensor networks [19]. It is based on a finite-precision
non-linear algebra supporting function multiplication and addition in the MPS
formalism, using it to implement a Chebyshev approximation scheme. The
resulting algorithm creates algebraic structures for individual functions of
one $f(x)$ or more variables $f(x,y,z)$, which can be composed in arbitrary
ways, such as $f(g(x))$ or $f(g(x)+h(y))$, etc. This approach has the benefit
of generalizing to multiple dimensions, going beyond examples that result from
tensorization of univariate functions [43, 37].

###  II.2 Finite-precision non-linear algebra

Figure 2: Illustration of the tensor product of two MPS $\ket{A}$ and
$\ket{B}$ in serial (a) and interleaved (b) orders.

Let us endow the space of matrix product states with three non-linear
operations. The fundamental one is the approximation of a linear combination
of MPS by another with potential limits on the size and complexity of its
tensors. Formally, a state

$\ket{\phi}=\operatorname{argmin}_{\phi\in\mathrm{MPS}}\left\|\ket{\phi}-\sum_{i=1}^{L}\alpha_{i}\ket{\psi_{i}}\right\|_{2}^{2},$
(7)

is searched over the set $\mathrm{MPS}$ of matrix product states with maximum
bond dimension $\chi$. This work uses a two-site variational truncation scheme
parameterized by a tolerance $\epsilon$ that controls the maximum error, as
detailed in Refs. [9, 16]. The second operation is the element-wise product of
MPS, resulting from the the non-linear multiplication of the vector elements

$\ket{u}\cdot\ket{v}:=\sum_{s_{1},s_{2},\ldots,s_{n}}u_{s_{1}s_{2}\ldots
s_{n}}v_{s_{1}s_{2}\ldots s_{n}}\ket{s_{1},s_{2},\ldots,s_{n}}.$ (8)

The third operation is the tensor product of two MPS. This is a
straightforward operation based on tensor concatenation that combines the two
Hilbert spaces. These tensors can be concatenated following a serial or
interleaved order [c.f. Figure 2].

$\displaystyle\ket{x}\otimes\ket{y}$
$\displaystyle:=\sum_{\vec{s},\vec{r}}x_{s_{1}\ldots s_{n}}y_{r_{1}\ldots
r_{n}}\ket{s_{1},\ldots,s_{n},r_{1},\ldots,r_{n}},$
$\displaystyle\ket{x}\otimes\ket{y}$
$\displaystyle:=\sum_{\vec{r},\vec{s}}x_{s_{1}\ldots s_{n}}y_{r_{1}\ldots
r_{n}}\ket{s_{1},r_{1},\ldots,s_{n},r_{n}}.$ (9)

The complexity of the resulting state depends on the structure of its
correlations [9].

##  III Chebyshev approximation

Any continuous univariate function can be approximated by a polynomial up to a
given order. The Chebyshev polynomials form a convenient basis[27, 28] that
results in a close to optimal approximation—the minimax polynomial—hence
providing excellent convergence guarantees, i.e., exponential and algebraic
decay of errors for analytical and differentiable functions. This technique is
commonly used within pseudospectral methods for the numerical solution of
partial differential equations [44]. In this section, this method is adapted
to approximate MPS/QTT representations of functions using the finite-precision
non-linear algebra framework from Section IIa. This adaptation involves
understanding the convergence guarantees, MPS complexity estimates, and
optimal evaluation orders.

###  III.1 Expansion on the Chebyshev basis

Figure 3: Representation of the first five Chebyshev polynomials.

The Chebyshev polynomials of the first kind [c.f. Fig. 3] are monic
polynomials $T_{k}(x):[-1,1]\rightarrow[-1,1]$ defined by the trigonometric
relation $T_{k}(\cos\theta)=\cos(k\theta),\quad k\in\mathbb{N}$. They are
orthogonal with respect to the weight function $\omega(x)=1/\sqrt{1-x^{2}}$

$\int_{-1}^{1}T_{k}(x)T_{l}(x)\omega(x)\mathrm{d}x=\frac{\pi}{2}(1+\delta_{k0})\delta_{kl}$
(10)

and they satisfy a three-term recurrence relation

$T_{k+1}(x)=2xT_{k}(x)-T_{k-1}(x),\;k\geq 1$ (11)

with $T_{0}(x)=1$, $T_{1}(x)=x$.

Any continuous function $f(x)$ defined in $[-1,1]$ shows a convergent series
in the Chebyshev basis, as

$\displaystyle f(x)$ $\displaystyle=\sum_{k=0}^{\infty}a_{k}T_{k}(x),\mbox{
with}$ (12) $\displaystyle a_{k}$
$\displaystyle=\frac{2-\delta_{k0}}{\pi}\int_{-1}^{1}f(x)T_{k}(x)\omega(x)\mathrm{d}x.$
(13)

The Chebyshev projection of order $d$ is the $d$-th order polynomial that
results from truncating (12) up to the first $d+1$ terms
$f_{d}(x)=\sum_{k=0}^{d}a_{k}T_{k}(x).$ These expansions can be trivially
generalized to functions in other intervals $g(u):[a,b]\to\mathbb{C}$, via the
affine transformation

$u=\alpha x+\beta,$ (14)

with $\alpha=(a-b)/2$, and $\beta=(a+b)/2$.

Alternatively, a function can be approximated by the interpolation on either
the Chebyshev-Gauss nodes, given by the $d$ zeros of the $d$-th order
Chebyshev polynomials $r_{k}\in(-1,1)$,

$r_{k}=\cos\left(\frac{\pi(2k-1)}{2d}\right),\quad k=1,2,3,\ldots,d,$ (15)

or the Chebyshev-Lobatto nodes, given by the $d+1$ extrema of $T_{d}(x)$ in
the closed set $[-1,1]$

$\ell_{k}=\cos\left(\frac{k\pi}{d}\right),\quad k=0,1,2,\ldots,d.$ (16)

The new interpolant is the polynomial series that results from evaluating the
interpolated function on the chosen set of nodes [28]. For the case of
Chebyshev-Gauss nodes,

$p_{d}(x)=\sum_{k=0}^{d}c_{k}T_{k}(x),\mbox{ with
}c_{k}=\frac{2}{d}\sum_{i=1}^{d}{}^{\prime}f(r_{i})T_{k}(r_{i}),$ (17)

where the primed sum involves halving the term $i=1$.

The coefficients $c_{k}$ can be computed from a discrete cosine transform
(DCT)

$c_{k}=\mathrm{DCT}\\{f(r_{k})\\},$ (18)

which may be implemented with $\mathcal{O}(d\log d)$ cost using the fast
Fourier transform.

These expansions may be generalized to the multidimensional scenario. However,
in such case the coefficients become a tensor $c_{k_{1}\ldots k_{m}}$ that
scales exponentially with the dimension and the contraction

$f(x_{1},\ldots,x_{m})=c_{k_{1}\ldots k_{m}}T_{k_{1}}(x_{1})\ldots
T_{k_{m}}(x_{m}),$ (19)

becomes computationally intractable. This complexity may be avoided for
functions that result from the composition of other functions, such as
$f(g(x)h(y))$. Also, the derivatives or primitives of a Chebyshev
approximation follow from simple recurrence relations of its coefficients
[27]. This enables the approximation of functions involving primitives or
derivatives, such as $f(g^{\prime\prime}(x)\int h(y))$, that need not have a
closed analytical form.

###  III.2 Convergence properties

The rate of convergence of polynomial approximations is usually measured in
terms of the distance in $L^{\infty}$ norm, also known as the Chebyshev norm,
between the function and its approximant. This distance is defined as the
limit of

$\displaystyle\|f(x)-g(x)\|_{p}$
$\displaystyle=\left(\int_{-\infty}^{\infty}|f(x)-g(x)|^{p}dx\right)^{1/p}$
(20)

when $p\rightarrow\infty$, and corresponds to the maximum pointwise error

$\|f(x)-g(x)\|_{\infty}=\max_{x}|f(x)-g(x)|.$ (21)

In practice, the functions are discretized on a grid with $N$ points,
resulting in the estimation

$\|f(x)-g(x)\|^{(N)}_{p}=\left(\frac{1}{N}\sum_{i=1}^{N}|f(x_{i})-g(x_{i})|^{p}\right)^{1/p}.$
(22)

The convergence of a Chebyshev interpolant (17) is determined by the decay
rate of its Gibbs oscillations in between the interpolated nodes. These depend
on the smoothness properties of the approximated function. Models with a
finite number of derivatives and analytical functions converge algebraically
and exponentially fast with the expansion order, respectively [28].

###### Theorem 1

For an integer $\nu\geq 0$, let $f$ and its derivatives through $f^{(\nu-1)}$
be absolutely continuous on $[-1,1]$ and suppose the $\nu$-th derivative
$f^{(\nu)}$ is of bounded variation $V$. Then for any interpolation order
$d>\nu$, its Chebyshev interpolants $p_{d}$ satisfy

$\left\|f-p_{d}\right\|_{\infty}\leq\frac{4V}{\pi\nu(d-\nu)^{\nu}}.$

###### Theorem 2

Let a function $f$ analytic in $[-1,1]$ be analytically continuable to the
open Bernstein ellipse $E_{\rho}$, where it satisfies $|f(x)|\leq M$ for some
$M$. Then for each interpolation order $d\geq 0$ its Chebyshev interpolants
$p_{d}$ satisfy

$\left\|f-p_{d}\right\|_{\infty}\leq\frac{4M\rho^{-d}}{\rho-1}.$

###  III.3 Clenshaw’s evaluation method

The evaluation of the partial sum (17) can pose numerical instabilities [45,
46, 27] and requires the additional computation of $d$ Chebyshev polynomials.
Both issues are simultaneously addressed by Clenshaw’s evaluation method [45].
This method allows for the efficient evaluation of partial sums of orthogonal
polynomials by expanding their three-term recurrence relation and keeping only
the non-vanishing terms. For Chebyshev polynomials, it reduces the evaluation
of the partial sum (17) to a linear combination $f(x)=y_{0}-xy_{1}$, where the
quantities $y_{j}$ are defined according to the following inverse recurrence
relation

$\displaystyle y_{d+1}$ $\displaystyle=y_{d}=0,$ (23) $\displaystyle y_{n-1}$
$\displaystyle=c_{n-1}-y_{n+1}+2xy_{n},\mbox{ for }2<n\leq d$ $\displaystyle
y_{1}$ $\displaystyle=c_{1}-y_{3}+2xy_{2},\mbox{ and }$ $\displaystyle y_{0}$
$\displaystyle=c_{0}-y_{2}+2xy_{1}.$

###  III.4 Chebyshev approximation with MPS/QTT

The partial sum (17) and its equivalent formulation using Clenshaw’s method
(23) require the identity and an initial condition $x$. These primitives can
be efficiently loaded in MPS using the recipe shown in Appendix A. Therefore,
the partial sum can be reformulated as

$\ket{p_{d}(x)}=\sum_{k=0}^{d}c_{k}\ket{T_{k}(x)},$ (24)

where the MPS Chebyshev polynomials are computed following their recurrence
relation

$\ket{T_{k+1}(x)}=2\ket{x}\cdot\ket{T_{k}(x)}-\ket{T_{k-1}(x)}.$ (25)

These expressions depend only on the identity state $\ket{1}$ and the initial
condition $\ket{x}$, as well as on linear combinations and product operations.
To avoid an unbounded growth of the bond dimension, the intermediate states of
the approximation shall be periodically simplified. Equivalently, the Clenshaw
evaluation (23) can be reformulated in terms of a sequence of MPS
$\ket{y_{n}}$ and a linear combination
$\ket{f}=\ket{y_{0}}-\ket{x}\cdot\ket{y_{1}}$. As shown in Appendix B, this
approach proves to be more efficient and robust than the standard evaluation.

Algorithm 1 MPS/QTT Chebyshev interpolation of $f(g(x))$ in the Chebyshev-
Gauss nodes

Univariate scalar function $f$, MPS initial condition $\ket{g}$ with support
in $[a,b]$, interpolation order $d$, tolerance $\epsilon$.

Compute the Chebyshev coefficients for $f$

$\tilde{r}_{k}\leftarrow\cos\left(\frac{\pi(2k-1)}{2d}\right),\quad
k=1,\ldots,d$

$r_{k}\leftarrow\textsc{affine}(\tilde{r}_{k},[-1,1]\rightarrow[a,b])$

$c_{k}\leftarrow\textsc{DCT}(f(r_{k}))$

Rescale the MPS argument

$\ket{\tilde{g}}\leftarrow\textsc{affine}(\ket{g},[a,b]\rightarrow[-1,1])$

Perform the Clenshaw evaluation

$n\leftarrow\mbox{number of qubits in}\ket{g}$

$\ket{\mathbf{I}}\leftarrow\mbox{Identity MPS with $n$ qubits}$

$\ket{y_{d+1}},\;\ket{y_{d}}\leftarrow\mbox{Zero MPS with $n$ qubits}$

for $k=d$ to 0 do

$\ket{y_{k}}\leftarrow\textsc{simp}(c_{k}\ket{\mathbf{I}}-\ket{y_{k+1}}+2\ket{\tilde{g}}\cdot\ket{y_{k+1}};\epsilon)$

end for

return $\textsc{simp}(\ket{y_{0}}-\ket{\tilde{g}}\cdot\ket{y_{1}};\epsilon)$

The steps required to compute the MPS Chebyshev approximation of a function
composition $f(g(x))$ on the Chebyshev-Gauss nodes (15) are summarized in
Algorithm 1. There, $g(x)$ is given as a prior MPS representation $\ket{g}$.
This algorithm specializes in loading functions $f(x)$ by making $g(x)=x$, and
generalizes as a method for function composition, extending earlier MPS
algebraic frameworks[19]. This approach has the caveat of needing to know the
support of $g(x)$ in order to rescale it to fit in the domain of the Chebyshev
expansion $[-1,1]$ using an efficient affine transformation. The algorithm
relies on the simplification routine, refered as simp and parametrized by a
tolerance $\epsilon$, to perform the linear combinations of MPS and their
subsequent simplification. Note that this approach may encounter numerical
instabilities for large MPS.

The cost of Chebyshev expansions and function composition is dominated by the
variational simplification routine, which scales as $\mathcal{O}(\chi^{3})$
for every tensor of the simplified state. Since the intermediate states are
polynomials of degree $k\leq d$, this sets a complexity upper bound of
$\mathcal{O}(\sum_{k=0}^{d}(k+1)^{3})=\mathcal{O}(d^{4})$. However, in
practice the bond dimension structure of the intermediate states allows for a
much better performance.

##  IV Numerical benchmarks

This section evaluates the performance of the MPS/QTT Chebyshev approximation
for loading collection of univariate and multivariate functions. The figures
of merit are the error in discretized $L^{\infty}$ norm (22), the maximum bond
dimension $\chi_{\text{max}}$ of the representation—a proxy for the memory
requirements given by the $\mathcal{O}(2\chi_{\text{max}}^{2})$ scaling—and
the runtime of the algorithm. These metrics are computed for different
parameter sets, defined by the number $n$ of qubits of the representation, the
polynomial order $d$ and the tolerance $\epsilon$ that parameterizes the
intermediate simplification steps. A comprehensive study requires nine
independent experiments. These are organized in three columns of mutually
shared axes.

The performance of MPS Chebyshev expansions will be compared with two state-
of-the-art techniques: multiscale interpolative constructions, as devised by
Lindsey [25], and tensor cross-interpolation (TCI) [23]. The description of
the TCI algorithm employed in this work, as well as the implementation
details, are presented in Appendix C.

###  IV.1 Univariate functions

Figure 4: Four qualitatively different functions proposed for the univariate
performance evaluation of the different methods. The color code labeling each
function is preserved along the subsequent studies. In blue, a Gaussian
probability distribution $f_{G}$ (26). In red, an oscillating function $f_{O}$
(27). In orange, the absolute value function $f_{A}$ (28). Lastly, in green,
the Heaviside step function $f_{S}$ (29). Figure 5: Results for the MPS/QTT
Chebyshev approximation of the four univariate functions $f_{G}$ (26), $f_{O}$
(27), $f_{A}$ (28) and $f_{S}$ (29). This study considers nine independent
experiments given by the combination of the three parameter settings with the
three figures of merit. All axes are mutually shared and shown in logarithmic
scale. The second row contains a thin dotted curve representing the maximum
bond dimension of a SVD decomposition. The results are asymptotically fitted
with thick dashed curves and the scalings are shown in the legend. The main
messages of the results can be understood by columns. The first column
suggests that MPS polynomial approximations are of low entanglement, as their
bond dimensions saturate to a value that is independent of the grid and only
function-dependent. The second column illustrates that the method converges
according to the theoretical rates for Chebyshev expansions, but with a
complexity that scales algebraically below the bound for MPS polynomials. The
third column demonstrates that the local truncation error of the MPS, given by
the tolerance $\epsilon$, provides a good proxy for the error and complexity
entailed in the approximations. Figure 6: Bond dimension structure of the MPS
corresponding to the oscillating function $f_{O}$ (27) for each of the 1300
intermediate Clenshaw evaluation steps. A similar behavior is noted and
omitted for the remaining test functions. The bond dimension varies noticeable
with the position in the chain, meaning that we never have to deal with a
square MPS tensor, and as a result the computational cost of the full
algorithm is much lower than the expected $\chi_{\textrm{max}}^{3}$ required
to perform an SVD decomposition of a
$\chi_{\textrm{max}}\times\chi_{\textrm{max}}$ matrix.  Figure 7: Results for
the three variants of multiscale interpolative constructions for the
oscillating function (27). The results of MPS/QTT Chebyshev approximation are
incorporated as a red solid baseline. The variation based on local
interpolation is repeated for local orders $\ell=1,10,30$. The results follow
a similar structure than the previous study for the MPS/QTT Chebyshev
approximation algorithm, shown in Figure 5. Figure 8: Results for the tensor
cross-interpolation of the four univariate functions $f_{G}$ (26), $f_{O}$
(27), $f_{A}$ (28) and $f_{S}$ (29). The results follow a similar structure
than the previous studies for the MPS/QTT Chebyshev approximation and
multiscale interpolative constructions, shown respectively in Figure 5 and
Figure 7. Figure 9: Benchmark for univariate function loading using MPS
Chebyshev approximations, multiscale interpolative constructions and TCI. The
results compare the required runtime and number of function evaluations
required to achieve a given approximation error. On the one hand, Figure (a)
compares the performance of MPS Chebyshev approximation with the three
variants of multiscale interpolative constructions for the oscillating
function (27). On the other hand, Figures (b) and (c) compare the performance
of MPS Chebyshev approximations (solid line) and TCI (dashed line) for the
four univariate functions, with respect to both the required runtime and
number of function evaluations.

All algorithms are tested on a setof qualitatively different one-dimensional
functions, shown in Figure 4. The first two are the Gaussian distribution

$f_{G}(x;\sigma)=\frac{1}{\sigma\sqrt{2\pi}}e^{-\left(\frac{x}{2\sigma}\right)^{2}}$
(26)

with $\sigma=1/3$, and a highly-oscillatory function

$f_{O}(x;\varepsilon)=\cos\left(\frac{1}{x^{2}+\varepsilon}\right)$ (27)

with $\varepsilon=10^{-2}$. Both are analytical functions and are expected to
exhibit an exponential convergence rate with the interpolation order $d$ after
a sufficiently large plateau of errors. The other two functions are non-
smooth, allowing to explore the limitations of polynomial approximations
compared to TCI. These include the absolute value function

$f_{A}(x;x_{c})=|x-x_{c}|,$ (28)

and the Heaviside step function

$f_{S}(x;x_{c})=\Theta(x-x_{c}).$ (29)

These functions are centered on an offset grid point
$x_{c}=\frac{1}{2}+\frac{1}{2^{5}}$ to avoid considering some artificially
simple representations. The Heaviside step function is discontinuous, so its
Chebyshev expansion exhibits constant Gibbs oscillations that lead to a finite
error. The absolute value function has only first-order derivatives of bounded
variation, leading to an inefficient algebraic attenuation of its Gibbs
oscillations.

####  IV.1.1 MPS/QTT Chebyshev approximation

The study explores grids with up to $2^{25}$ points—i.e., $n=25$ qubits—and a
numerical tolerance of up to $\epsilon=10^{-14}$, leading to a truncation
error close to machine precision. The interpolation order $d$ runs up to 1300
for all functions, except for the Gaussian function, which requires a maximum
of $d=50$ in order to prevent a degradation of the accuracy of the Clenshaw
method, as detailed in Appendix B. The results are grouped in Figure 5, where
they are addressed and interpreted column-wise.

The first column addresses the representability of the functions as MPS by
setting the order of the expansion and tolerance $\epsilon$ to the most
precise values, and evaluating the effects of the grid size on the
$L^{\infty}$ norm, maximal bond dimension and CPU time, as respectively shown
in Figures 5(a), 5(b) and 5(c). The results depict that the $L^{\infty}$ norm
is a smooth function of the grid spacing and converges for grids involving a
few thousand points, for which the maximal bond dimension saturates to a value
that only depends on the specific function—a manifestation of the bounded
entanglement content of polynomials—. Also, it matches those of the optimal
decomposition of the tensor $f_{s_{1}s_{2}\ldots s_{n}}$ using the iterative
SVD for the two analytical functions. The computational cost scales
sublinearly, with an averaged rate of $\mathcal{O}(n^{0.64\pm 0.16})$, as the
tensors corresponding to the less significant bits asymptotically converge to
a reduced bond dimension [25]. This decay can be interpreted to mean that, at
those scales, the polynomial approximations of sufficiently well-behaved
functions can be interpolated almost linearly, even if we are approximating
non-smooth functions. This decay is demonstrated for every intermediate step
in Figure 6.

The second column in Figure 5 evaluates the effects of the interpolation order
for the finest grid with $2^{25}$ nodes. Figure 5(d) demonstrates that the
approximations converge according to the theoretical expectations of Theorems
1 and 2. Moreover, Figure 5(e) shows that the maximum bond dimension scales
sub-linearly with the interpolation order as $\mathcal{O}(d^{0.28})$ for all
test functions. This curve lays systematically below the polynomial upper
bound of $\chi_{\text{max}}=d+1$ known for $d$-th order polynomials [c.f.
Appendix A]. At convergence, this bond dimension is consistent with the
optimal SVD value, as shown for the Gaussian and oscillatory functions for
$d_{G}=50$ and $d=1300$, respectively. Finally, Figure 5(h) shows how the
computational cost scales at a close-to-linear rate of $\mathcal{O}(d^{1.37})$
with the interpolation order, well below the expected worst case scenario of
$\mathcal{O}(d^{4})$. Again, this is attributed to the intermediate bond
dimension structure of the MPS approximations shown in Figure 6. As the MPS
tensors are highly rectangular, the expected $\mathcal{O}(d^{3})$ cost for SVD
is subleading.

Finally, the third column analyzes whether the truncation error at each step
acts as an indicator of the global approximation error. Figures 5(g), 5(h) and
5(i) illustrate that the global approximation error, maximal bond dimension
and CPU time are parameterized by the tolerance parameter $\epsilon$. In
particular, the global error is a monothonically decreasing function of
$\epsilon$ until it is dominated by the error of the expansion. Furthermore,
the maximal bond dimension increases mildly—logarithmically—with the tolerance
parameter for all functions. This parametrization allows exploring the
tradeoffs between accuracy and computational cost when a good Chebyshev
expansion is warranted.

Based on these results, the MPS/QTT Chebyshev approximation proves to be an
efficient method for loading highly-differentiable univariate functions. It
offers a reliable performance for a large numbers of qubits, attributed to the
advantageous bond dimension structure of the intermediate steps and
demonstrated by its sub-logarithmic scaling with discretization size.
Moreover, it exhibits excellent error guarantees for smooth and highly-
differentiable functions that align with the theoretical bounds. This
motivates its generalization to multivariate problems. However, before
addressing this scenario, it is of interest to compare it against other state-
of-the-art methods for loading functions, such as multiscale interpolative
constructions [25] and TCI.

####  IV.1.2 Multiscale interpolative constructions

In the one dimensional case there exist alternative methods to construct the
MPS representation of a Chebyshev interpolant. Notably, M. Lindsey has
developed several variants of a multiscale interpolative construct based on
Lagrange interpolation [25], which shall be compared with the previous
iterative algorithm.

In the following, we briefly summarize the spirit of this technique, and refer
the interested reader to the original paper. Essentially, the basic construct
builds an MPS interpolant in a single pass. It begins by assembling a left-
most tensor core that captures the information of the values of a function
$f(x)$ on two Chebyshev-Lobatto grids (16), each with $d+1$ points and sharing
the intersection. This procedure constructs an effective interpolant by
joining two interpolants of order $d+1$ and using $2d+1$ function values.
Then, two additional tensor cores are assembled belonging to the bulk and
right-most edges of the MPS. These cores propagate the function information
and implement the Lagrange interpolating polynomial. Once they are assembled,
these cores can be reused for subsequent interpolations. Hence, this
interpolation procedure is universal, given that different functions only
differ in the first tensor while all the subsequent ones are the same. The
resulting state has all bond dimensions set to $\chi=d+1$ at maximum, which is
a large overestimation and benefits from a simplification at the end of the
algorithm. However, the same fundamental algorithm may be implemented in a
rank-revealing form, where the tensor cores from the basic construct are
simplified along the way using the SVD with a given tolerance $\epsilon$,
avoiding the need for a large-scale final simplification. Additionally, this
algorithm can be reformulated using local Lagrange interpolation with a
specified local order $\ell$. Then, the tensor cores become sparse improving
the efficiency of tensor decompositions and contractions.

The parameters of all these variants are similar to those of the previous
iterative algorithm, which enables a similar study. However, for simplicity
and in order to analyze all algorithmic variants, the study focuses only on
the oscillatory function (27), which is the most complex example given for the
MPS Chebyshev expansion. The results are presented in Figure 7, where the
outcome of the previous study is incorporated as a reference line in red color
with the label “Iterative”.

As shown in Figure 7(a), global interpolative constructions converge
accurately and stably with the number of qubits. However, local constructions
require a large enough local order of approximately above $\ell=10$ to do so,
given that small enough local orders lead to a large deterioration of the
convergence and maximum bond dimension as shown in Figures 7(a) and 7(b).
Interestingly, Figures 7(d) and 7(e) demonstrate that these constructions
converge with half the polynomial order than MPS Chebyshev expansions. This is
due to the non-linear convergence rates of the two interpolants at each
Chebyshev-Lobatto grid, spanning each half of the interval with
$\lfloor\frac{d}{2}\rfloor-1$ nodes.

The complexity of the interpolative construct is dominated by the assembly of
the left-most, bulk and right-most cores. Hence, the algorithm’s performance
almost does not scale with the system size, as shown in Figure 7(c). However,
the resources required for the assembly scale with the interpolation order $d$
with a rate that depends on the interpolation variant. As illustrated in
Figure 7(f), the rate for the rank-revealing construction is close to qubic,
as $\mathcal{O}(d^{2.73})$, while local constructions exhibit a slower rate,
concretely of $\mathcal{O}(d^{1.25})$ for $\ell=30$. The former rate surpasses
the close-to-linear rate of $\mathcal{O}(d^{1.37})$ of the previous algorithm
shown in Figure 5(h). Nevertheless, we note that these rates may depend on
specific implementation choices. In any case, Figures 7(c) and 7(i) suggest
that the interpolative constructions outperform the iterative expansion for
this function for $n\geq 10$, exhibiting lower CPU times within one order of
magnitude. Finally, Figures 7(g), 7(h) and 7(i) illustrate how these
algorithms can be similarly parameterized by the tolerance parameter of the
simplification routine.

####  IV.1.3 Tensor cross-interpolation (TCI)

Tensor cross-interpolation (TCI) is a state-of-the-art technique that
constructs the tensor-train representation of a black-box function by sampling
its elements. As opposed to the previous algorithms, TCI does not create a
polynomial or algebraic structure that can be reused, but rather samples the
target function (or tensor) to get the smallest tensor cores that minimize a
given cost function. The details of TCI as well as the particular
implementation used in this work are detailed in Appendix C.

Essentially, TCI iteratively optimizes an initial MPS approximation of the
function—which may start from a random guess. This optimization is done
following “sweeps” that increase its bond dimension until convergence is
achieved. In our suboptimal but performant implementation, the bond dimension
grows with each sweep by 2, and is bounded by a predefined threshold value
$\chi_{\text{thr}}$. The relation between the order of polynomials and the
maximum bond dimension [c.f. Appendix A] suggests that the threshold
$\chi_{\mathrm{thr}}$ may be considered analogous to the polynomial order $d$
of the previous methods. Despite this upper limit, the algorithm is designed
to halt when the approximation error in $L^{\infty}$ norm (21) falls below a
predefined value. Interestingly, as TCI does not perform a polynomial
approximation, it is not constrained by the smoothness of the target function
and can in fact load non-differentiable targets efficiently.

This study replicates the analysis from Sect. IV.1.1 up to a maximum bond
dimension $\chi_{\text{thr}}=30$ and a halting approximation error of
$\epsilon$. By default, the resulting MPS is not simplified with a tolerance
$\epsilon$. However, a final simplification with tolerance $\epsilon$ is
considered in order to replicate the results shown in the third column of the
previous studies. As the outcome now depend on a random initial condition,
they are evaluated 10 times to extract the mean and standard deviation.

The results are presented in Figure 8. As shown in Figures 8(a) and 8(d), TCI
achieves machine precision convergence for all test functions, even for non-
differentiable targets, exhibiting an exponential convergence with greater
rates than the previous algorithms. Figures 8(b) and 8(e) show how the
complexities resemble those computed with SVD, while Figures 8(c) and 8(f)
illustrate that asymptotic time scaling of TCI is close to linear with both
the number of qubits and the threshold bond dimension. Hence, the required
computational time shows an averaged scaling rate of $\mathcal{O}(n^{1.11\pm
0.1})$ with the number of qubits, which is less favorable than the sublinear
rate of the MPS/QTT Chebyshev approximation and the constant rate of
multiscale interpolative constructions.

####  IV.1.4 Performance benchmark

The previous studies suggest that the three algorithms are capable of loading
smooth univariate functions in MPS accurately and efficiently, avoiding the
curse of dimensionality. However, TCI provides an exponential advantage in the
representation of non-smooth functions. Moreover, the algorithms show
different scaling rates. This section compares their performance explicitly.

The results are presented in Figure 9. Figure 9(a) compares the runtime
required to achieve a given error between the iterative MPS/QTT Chebyshev
approximation and the multiscale interpolative constructions for the
oscillating function (27) for different values of the order $d$. Similarly,
Figures 9(b) and 9(c) benchmark MPS Chebyshev expansions and TCI for all test
functions, also considering the number of required function evaluations.

As expected, all three methods converge exponentially for the analytical
functions but differ in the prefactors of the overall scaling. Figure 9(a)
illustrates how the constructive methods outperform the iterative Chebyshev
approximation for the oscillating function within one order of magnitude.
Similarly, Figure 9(b) shows that TCI outperforms the iterative method within
more than one order of magnitude, due to its superior convergence rate. Also,
it shows an exponentially better performance for the non-smooth functions.
However, it underperforms for the Gaussian function, which requires a much
lower interpolation order to converge. The performance gaps can be
qualitatively attributed to the total information brought into the sampling.
As shown in Figure 9(c), TCI uses approximately two orders of magnitude more
evaluations than the iterative expansion.

Despite these performance gaps, all three approximation techniques load smooth
univariate functions with comparable accuracy. In this context, the function
composition capabilities of the iterative Chebyshev algorithm enable tasks
beyond function loading. These capabilities have no straightforward analogue
using interpolative constructions, and require expensive state samples
involving costly tensor contractions for the case of TCI. In contrast, the
iterative Chebyshev approximation requires fewer inexpensive function
evaluations.

###  IV.2 Multivariate functions

Figure 10: Results for multivariate function loading using MPS Chebyshev
approximation and TCI. Columns 1 and 2 illustrate the results for the MPS
Chebyshev approximation of the multivariate functions $f_{\text{product}}$
(30) and $f_{\text{squeezed}}$ (31). Column 3 shows the results for the tensor
cross-interpolation of the two smooth Gaussian functions together with the
multivariate non-smooth functions $f_{A,m}$ (32) and $f_{S,m}$ (33). The
circle and triangle markers respectively correspond to the serial and
interleaved qubit orders.

The performance of the MPS/QTT Chebyshev approximation for loading univariate
functions motivates its generalization to the multivariate scenario. In
general, this requires creating an $m$-dimensional tensor of coefficients with
an exponential cost (19). However, this can be avoided if the functions have a
simple algebraic structure that can be derived from the composition of
univariate functions. The goal of this study is to confirm the applicability
of the MPS function algebra from Sect. II and the Chebyshev expansion to
achieve a complete algebraic framework for function loading [19] that
generalizes to more complex arbitrary functions and MPS/QTT encodings, while
at the same time comparing the performance with the state-of-the-art.

This hypothesis is explored using two target functions: an ordinary Gaussian

$f_{\text{product}}(\mathbf{x})=\exp\left(\sum_{i=1}^{m}x_{i}^{2}\right),$
(30)

and a squeezed Gaussian

$f_{\text{squeezed}}(\mathbf{x})=\exp\left(\left[\sum_{i=1}^{m}x_{i}\right]^{2}\right).$
(31)

Both functions can be respectively understood as the composition of the
exponential $h(x)=\exp(-x)$ with two polynomials, which are
$\sum_{i}x_{i}^{2}$ and $\left(\sum_{i}x_{i}\right)^{2}$. Hence, both can be
obtained by evaluating the MPS/QTT Chebyshev expansion of $h(x)$ on either of
these polynomials—which have an efficient MPS representation combining the
recipe of Appendix A with tensor products.

The study explores the performance of this composition strategy on problems
with growing dimensionality, following a similar structure to previous studies
with univariate functions. The dimension of the function, $m$, ranges between
1 and 10, and the tolerance of the intermediate simplifications is varied
between $10^{-1}$ and $10^{-14}$. The results are compared against TCI for the
same multivariate functions. However, since TCI enables loading non-smooth
functions accurately, we also benchmark this method against the multivariate
generalization of $f_{A}$ (28) and $f_{S}$ (29),

$\displaystyle f_{A,m}(\mathbf{x};x_{c})$
$\displaystyle=\left|\sum_{i=1}^{m}x_{i}-x_{c}\right|,\mbox{ and}$ (32)
$\displaystyle f_{S,m}(\mathbf{x};x_{c})$
$\displaystyle=\Theta\left(\sum_{i=1}^{m}x_{i}-x_{c}\right).$ (33)

Note that, as explained in Sect. II, there are many different orderings of
qubits in a multidimensional tensorized encoding. The study analyzes
sequential and interleaved orders from Fig. 2.

The largest problem involves 10 variables with 20 qubits each, resulting in a
Hilbert space with $\sim 10^{60}$ elements. This makes the computation of the
error in discretized $L^{\infty}$ norm intractable. Hence, this error is
estimated using standard Monte Carlo sampling techniques. This estimation is
assessed to converge to the full-tensor error for a moderate amount of
function samples in Appendix D. In the following, the approximation error is
estimated as the average of 10 measurements with 1000 random samples each. The
results are presented in Figure 10.

The first column addresses the performance of the iterative Chebyshev
approximation for the two multivariate Gaussian functions with respect to
their dimension, in both the serial and interleaved qubit orders. As suggested
by Figures 10(a) and 10(b), the accuracy and complexity of the method now
largely depends on the qubit order. The product Gaussian distribution (30) is
stably loaded in the serial order with a reduced bond dimension, but has an
exponentially growing bond dimension in the interleaved order as
$e^{1.57m}\approx 5^{m}$, preventing any approximation beyond $m=4$. On the
other hand, the squeezed Gaussian distribution (31) is loaded efficiently in
both the serial and interleaved orders, although more stably and accurately in
the latter, exhibiting a sublinear bond dimension growth of
$\mathcal{O}(m^{0.49})$ and implying an exponential advantage over its tensor
counterpart. In all cases, the approximation error seems to be largely
correlated to the bond dimension. Figure 10(c) shows that the CPU time scales
polynomially cases except for the Gaussian with exponentially growing bond
dimension. Interestingly, the squeezed function in interleaved order shows a
superior performance than in serial order. This is in part due to its enhanced
stability, which allows for a more aggresive simplification of the initial
condition.

The second column addresses the parametrization of the results with the local
truncation error induced by the tolerance $\epsilon$, for the most stringent
parameter setting of $n=20$ and $m=10$. The inefficient interleaved product
Gaussian function is omitted. Figures 10(d) and 10(e) show how the tolerance
parameterizes the error and maximum bond dimension similarly than in the
univariate scenario. However, now the stability of the maximum bond dimension
is enhanced in the interleaved qubit order. Figure 10(f) shows how this
parameterization reduces the CPU time by more than one order of magnitude at
the cost of a larger approximation error.

Figure 11: CPU time required to reach a given approximation error for the
MPS/QTT Chebyshev approximation method and TCI. The results for the Chebyshev
approximation are parameterized with the tolerance $\epsilon$ between
$10^{-1}$ and $10^{-14}$. The large isolated markers correspond to TCI
results, while the curves correspond to Chebyshev results. The functions are
labeled following the legend of Figure 10(a), which is omitted for
conciseness. Figure 12: Number of evaluations required by TCI for the four
multivariate functions with respect to their dimension. The functions are
labeled acording to the legend shown in Figure 10(a), which is omitted for
conciseness.

Finally, the third column addresses the representability of these functions,
including the two additional non-smooth functions (32) and (33) for the TCI
algorithm. Figures 10(g) and 10(h) confirm that TCI is able to load accurately
both smooth and non-smooth functions, given the appropriate qubit order. In
particular, it shows how the non-smooth functions are exponentially more
accurate and efficient in the interleaved order, while the remaining functions
follow a similar behavior as before. However, the convergence of the
multivariate absolute value function starts to degrade for a dimension larger
than 5. Also, we note that the sampled estimation for the approximation error
may not be as reliable for non-smooth functions, whose errors are expected to
be larger close to the discontinuity. In the appropriate qubit order, the
maximum bond dimension scales with an approximately linear rate with the
dimension. This induces a polynomial scaling of the required computational
resources, as depicted in Figure 10(i), overcoming the curse of
dimensionality.

The tradeoff between performance and accuracy exhibited by MPS Chebyshev
approximations enable them to beat the performance of TCI in some multivariate
scenarios. Both methods are benchmarked in Figure 11 comparing the CPU time
required to reach a given approximation error. These results illustrate how
the Chebyshev approximation exhibits a better performance for the squeezed
Gaussian (31) in interleaved order, and can achieve a competitive performance
for the remaining scenarios for a sufficiently large local truncation
tolerance—while still achieving values that are useful in practice, below
$10^{-6}$ in all remaining cases.

Moreover, multivariate MPS Chebyshev approximations require many orders of
magnitude less function evaluations than TCI. This is illustrated in Figure
12, showing that the amount of function evaluations required by TCI reaches
several million, while Chebyshev approximations require of the order of
$\mathcal{O}(d_{\text{max}})$ evaluations. These results suggest that, as
opposed to MPS Chebyshev approximations, TCI may be less suitable for
scenarios where the function evaluations pose a significant overhead. A
relevant example is given by the composition of functions in generic states,
for which function evaluation involves costly tensor contractions.

##  V Discussion

This work has explored the problem of loading univariate and multivariate
functions in matrix product states or quantized tensor trains, motivated by
their application in quantum-inspired numerical methods. To this end, an
adaptation of the Chebyshev approximation algorithm has been proposed for
MPS/QTT structures. This method is based on an non-linear algebraic framework
for MPS/QTT polynomial addition and multiplication, which enables the
translation of Chebyshev expansions and Clenshaw evaluation schemes.

The numerical studies illustrate that the Chebyshev MPS/QTT iterative
algorithm is stable, achieves the theoretical convergence rates and requires
polynomial resources both in time and memory. In particular, in the univariate
case the method scales linearly with the interpolation order and sublinearly
with the number of qubits. This scaling generalizes to multivariate functions,
which also can be encoded using polynomial resources.

Furthermore, MPS Chebyshev approximations provide a generic framework for
function composition on MPS/QTT that enables non-linear manipulation of
generic states. These capabilities are useful for the solution of non-linear
equations with quantum-inspired methods [30] and the implementation of filters
in Hamiltonian problems [31].

The following caveats are identified. First, the method is limited to the
composition of low-dimensional, highly-differentiable functions, as the volume
of the tensor of Chebyshev coefficients (19) scales exponentially with
dimension. Second, the method generally exhibits slower convergence rates than
multiscale interpolative constructions—requiring an interpolation order twice
as large—and TCI. Finally, the required intermediate simplifications and
linear combinations may be susceptible to numerical instabilities for large
enough states.

Lastly, the following extensions are considered. On the one hand, the
iterative scheme proposed for MPS Chebyshev expansions may be applicable for
more elaborate algebraic schemes, such as rational approximation. These
schemes may be able to extend the domain of applicability to piecewise smooth
functions [47]. On the other hand, the exploration of function composition on
irregular discretizations would be of interest as it enables efficient
integration strategies, such as those based on Clenshaw-Curtis [48] or Gauss-
Kronrod [49] quadrature rules.

The algorithms and methods shown in this work are freely available in the
SeeMPS library [50], an open-source Python framework that aims to ease the
prototyping and experimentation with MPS/QTT quantum-inspired algorithms. This
library has supported other related applications [9, 16]. All the simulations
described in this study were conducted on a workstation equipped with an
18-core Intel(R) Xeon(R) W-2295 CPU and 64GB of RAM. The source files required
to replicate the simulations are accessible in the public repository [51].

###### Acknowledgements.

This work has been supported by Spanish Projects No. PID2021-127968NB-I00 and
No. PDC2022-133486-I00, funded by MCIN/AEI/10.13039/501100011033 and by the
European Union “NextGenerationEU”/PRTR”1. The authors acknowledge the
discussions on related topics with A. Bou Comas and E. Staffetti. The authors
also gratefully acknowledge the Scientific Computing Area (AIC), SGAI-CSIC,
for their assistance while using the DRAGO Supercomputer for performing the
simulations, and Centro de Supercomputación de Galicia (CESGA) who provided
access to the supercomputer FinisTerrae. The authors thank the Institut Henri
Poincaré (UAR 839 CNRS-Sorbonne Université) and the LabEx CARMIN
(ANR-10-LABX-59-01) for their support. PGM acknowledges the funds given by
“FSE invierte en tu futuro” through an FPU Grant FPU19/03590 and by
MCIN/AEI/10.13039/501100011033.

##  Appendix A Analytical encoding of polynomials in MPS/QTT

Figure 13: Decomposition of a parametrized collection of monomials $x^{i}$
contracted with a vector of coefficients $c_{i}$ in a MPS.

This Appendix provides a recipe to construct the MPS/QTT representation of the
polynomials $f_{d}(x)$ of degree $d$

$f_{d}(x)=\sum_{i=0}^{d}p_{i}x^{i}$ (34)

over a finite domain $x\in[a,b]$. The constructive proof creates a tensor
representation for the discretization of $f$ on a grid of arbitrary density,
with tensors $A_{i,r}^{s}\in\mathbb{C}^{\ell\times(d+1)^{2}}$ whose index
sizes are determined by polynomial degree $\chi=d+1$ and the quantization
$\ell$ ($\ell=2$ in our simulations).

The proof is done by induction. Starting from $\ket{x^{i}_{n}}$ , the tensor
decomposition of a monomial $x^{i}$ over a grid with $\ell^{n}$ components,
the construct derives a recipe for $\ket{x^{i}_{n+1}}$ by adding one more
variable:

$\ket{x^{i}_{n+1}}=\sum_{r=0}^{d}A_{ri}^{s_{n+1}}\ket{x^{r}_{n}}\ket{s_{n+1}}.$
(35)

The new discretization can be expressed in terms of the original, coarser-
grained one $x(s_{1},s_{2},\ldots,s_{n})$ and a displacement that depends on
the new quantization variable $s_{n+1}\in\\{0,1,\ldots,\ell-1\\}$

$x(\mathbf{s}_{n+1})=x(\mathbf{s}_{n})+\frac{d}{\ell^{n+1}}s_{n+1}.$ (36)

Hence, it follows

$\displaystyle x(\mathbf{s}_{n+1})^{i}$
$\displaystyle=\sum_{r=0}^{i}x(\mathbf{s}_{n})^{r}\left(\frac{s_{n+1}}{\ell^{n+1}}\right)^{i-r}{{i}\choose{r}}$
(37)
$\displaystyle=\sum_{r=0}^{d}x(\mathbf{s}_{n})^{r}\left(\frac{s_{n+1}}{\ell^{n+1}}\right)^{i-r}{{i}\choose{r}}\Theta(i-r)$
(38) $\displaystyle=:\sum_{r=0}^{d}x(\mathbf{s}_{n})^{r}A_{r,i}^{s_{n+1}}.$
(39)

This expression uses the asymmetric Heaviside function, which is nonzero
$\Theta(i-r)=1$ only when $i\geq r$.

Equation (39) is a recursive formula to construct the tensor representation of
$\ket{x^{i}_{n}}$

$\ket{x^{i}_{n}}=A_{0\alpha_{1}}^{s_{1}}A_{\alpha_{1}\alpha_{2}}^{s_{2}}\cdots
A_{\alpha_{N-1}i}^{s_{n}}\ket{s_{1},s_{2},\ldots,s_{n}},$ (40)

with tensors

$A_{r,k}^{s_{n}}=\left(\frac{s_{n}}{\ell^{n}}\right)^{k-r}{{k}\choose{r}}\Theta(k-r),$
(41)

where the indices $r,k,\alpha_{j}$ run between 0 and $d$. This construct is a
tensor network with one free index, which is the degree of the monomial, as
illustrated in Fig. 13a.

A similar proof allows building a reversed construct

$\ket{x^{i}_{n}}=B_{i\alpha_{1}}^{s_{1}}B_{\alpha_{1}\alpha_{2}}^{s_{2}}\cdots
B_{\alpha_{N-1}i}^{s_{n}}\ket{s_{1},s_{2},\ldots,s_{n}},$ (42)

with flipped tensors

$B_{k,r}^{s_{n}}=\left(\frac{s_{n}}{\ell^{n}}\right)^{k-r}{{k}\choose{r}}\Theta(k-r).$
(43)

The interest of these constructs is that they extend to arbitrary polynomials
of a fixed degree. Thus, for any combination of monomials of degree $d$ or
lower we acquire a tensor network decomposition with bonds of size $d+1$

$p(x)=\sum_{i}c_{i}x^{i}\;\Rightarrow\;\ket{p(x)_{n}}=\sum_{i}c_{i}\ket{x^{i}_{n}},$
(44)

from either Eq. (40) or (42). This construct only modifies one of the tensors
in the expansion [c.f. Fig. 13b].

##  Appendix B Performance of the Clenshaw evaluation method

Figure 14: Comparison between the Clenshaw evaluation (23) (in red) and the
direct evaluation (in brown) of the partial sum (24) for the MPS Chebyshev
approximation of the oscillating function (27). This study follows the same
structure as the numerical studies performed in Section IV, showing individual
experiments with mutually shared axes in logarithmic scale. Figure 15: Maximum
bond dimension for the $d$ intermediate states of a MPS approximation for the
Gaussian (26) with $n=25$, $d=200$ and $\epsilon=10^{-14}$. The dotted line
represents the polynomial bound $\chi\propto d+1$.

This Appendix compares the two possible ways to evaluate the partial sums (12)
and (17) for MPS/QTT Chebyshev approximation. The former is based on
evaluating the Chebyshev polynomials $T_{k}(\ket{g})$ using the recurrence
relation (25) and then directly evaluating the partial sums. The latter is
based on the Clenshaw summation method, based on evaluating (23).

In most occasions, both evaluation methods have similar accuracy, but the
Clenshaw method is more efficient, with smaller runtimes and maximum bond
dimensions, provided the expansion uses the right number of polynomials. If
the expansion goes beyond the required order, introducing contributions whose
weight is below machine precision, the algorithm will lead to an unbounded
growth of the intermediate bond dimensions, resulting in a significantly
longer runtime. However, since the optimal order can be estimated beforehand
through the decay of the Chebyshev coefficients, this should not pose a
practical problem.

A suitable benchmark to compare both evaluation orders is the Chebyshev
expansion of the oscillating function (27), using the most stringent
parameters $n=25$, $d=1300$ and $\epsilon=10^{-14}$. Figure 14 shows how the
Clenshaw method leads to significantly smaller bond maximum bond dimension and
runtimes. Similar accuracies between both methods are noted and omitted for
conciseness.

Figure 15, shows the maximum bond dimension of the intermediate steps,
$\chi_{\text{i, max}}$, for the interpolation of the Gaussian function (26) up
to an order $d=200$. The Clenshaw method overestimates the intermediate bond
dimensions when the interpolation order is also overestimated—i.e., the
expansion included terms with a contribution below machine precision—.
Interestingly, the direct evaluation does not exhibit this problem, which in
this case would imply a better performance. However, this advantage ceases to
exist when the required expansion order is properly estimated beforehand.

##  Appendix C Tensor cross-interpolation

Tensor cross-interpolation (TCI), also known as “TT-cross”, is a well-known
algorithm that computes the tensor-train representation of a black-box
multivariate function from the sampling of its elements along specific
patterns known as crosses. This function can implicitly represent a tensor, so
TCI specializes as an efficient tensor decomposition method. Other potential
applications of TCI include function composition and loading functions in MPO
structures.

TCI exhibits a polynomial time complexity and offers an exponential advantage
over the standard tensor decomposition method based on iterated SVD. This
establishes it as one of the preferred method for loading highly multivariate
functions. It has demonstrated a remarkable performance in multivariate
integration [39, 52, 20, 53], with applications in quantum field theory [20]
and condensed matter physics [42, 54], showing exponential improvements over
state-of-the-art integration techniques [39, 53]. Moreover, it has proven
useful for other tasks including multivariate optimization [55, 56],
multivariate distribution sampling [57] or quantum state tomography [58].

TCI was introduced in a seminal work by Oseledets [23] and conceived as a
generalization of a low-rank matrix decomposition, known as the skeleton
decomposition [59], to tensors. This early version of TCI, however, was not
rank adaptive and required prior knowledge of the tensor ranks. Since its
introduction, numerous advancements have been made to improve the efficiency
and rank adaptivity of TCI. One such advancement is the incorporation of two-
site optimizations, similar to those used in the Density Matrix
Renormalization Group (DMRG) algorithm [40]. Other approaches have introduced
greedy optimization strategies that involve the use of “pivots” [41, 39, 20].
Additionally, rank-adaptive one-site optimizations utilizing the rectangular
skeleton decomposition [60] have been proposed.

The implementation of TCI used in this work follows this latter method, and is
based on a former tensor-train implementation employed in Refs. [61, 55].

###  C.1 The skeleton decomposition

The skeleton decomposition is a low-rank matrix decomposition of a matrix $A$
in a column submatrix $C$, a row submatrix $R$ and a submatrix $\widehat{A}$
formed by the intersection of the chosen columns and rows, as

$A=C\widehat{A}^{-1}R.$ (45)

This decomposition is exact for rank-$r$ matrices if $r$ linearly independent
rows and columns are chosen. In the approximate case, its accuracy is
dependent on the chosen set of rows and columns. A notation that includes the
sets of indices explicitly is adopted for convenience. Given the sets
$\mathcal{I}=\\{i_{1},\ldots,i_{r}\\}$ and
$\mathcal{J}=\\{j_{1},\ldots,j_{r}\\}$, each containing $r$ indices for rows
and columns respectively, the skeleton decomposition can be reformulated as

$A=A(\mathbb{I},\mathcal{J})\left[A(\mathcal{I},\mathcal{J})\right]^{-1}A(\mathcal{I},\mathbb{J}),$
(46)

as depicted in Figure 16. Here, $\mathbb{I}$ and $\mathbb{J}$ are the sets of
all rows and columns of $A$.

Figure 16: Skeleton decomposition of the matrix $A$, using tensor network
notation (a) and explicitly showing the matrix elements (b). The selected rows
and columns (blue) are contained in the index sets $\mathcal{I}$ and
$\mathcal{J}$, resulting in intersections (pink) and pivots (red).

.

Finding the optimal set of $r$ rows and columns for the skeleton decomposition
is challenging [62]. However, there are effective heuristics for achieving a
good approximation [59]. These heuristics follow the “maxvol principle”, which
suggests that an accurate approximation is achieved by maximizing the
“volume”—defined as the absolute value of the determinant—of the intersection
submatrix $\widehat{A}$, as

$\operatorname{vol}(\widehat{A}):=|\operatorname{det}\widehat{A}|.$ (47)

This principle is generalizable to cases where the intersection submatrices
are rectangular, allowing for the selection of rectangular submatrices of
maximum volume [60].

###  C.2 Generalization to tensors

The skeleton decomposition can be applied on the “unfolding” matrices of a
tensor. These matrices are given by the reorganization of the tensor’s indices
in two distinct multi-indices. Given a tensor $A(s_{1},s_{2},\ldots,s_{n})$,
where each index $s_{j}$ has a corresponding size $\mathbb{I}_{k}$, its
indices before and after any chosen position $k$ can be grouped into two
multi-indices,

$\displaystyle s_{<k}:=s_{1}s_{2}\ldots s_{k-1}$
$\displaystyle\in\mathbb{I}_{1}\times\mathbb{I}_{2}\times\ldots\times\mathbb{I}_{k-1},$
(48) $\displaystyle s_{>k}:=s_{k+1}s_{k+2}\ldots s_{n}$
$\displaystyle\in\mathbb{I}_{k+1}\times\mathbb{I}_{k+2}\times\ldots\times\mathbb{I}_{n}.$

Then, by grouping the index $s_{k}$ with the multi-index $s_{<k}$, it can be
unfolded in a matrix

$A(s_{1},s_{2},\ldots,s_{n}):=A(s_{<k}s_{k},s_{>k}).$ (49)

The skeleton decomposition of this matrix computes the sets multi-indices
$\mathcal{I}_{\leqslant k}$ and $\mathcal{J}_{>k}$, that are respectively in
$s_{<k}s_{k}$ and $s_{>k}$, which maximize the volume of the intersection
matrix:

$\displaystyle A\left(s_{<k}s_{k},s_{>k}\right)\approx$ (50)
$\displaystyle\approx
A\left(s_{<k}s_{k},\mathcal{J}_{>k}\right)\left[A\left(\mathcal{I}_{\leqslant
k},\mathcal{J}_{>k}\right)\right]^{-1}A\left(\mathcal{I}_{\leqslant
k},s_{>k}\right).$

The tensor-train decomposition of the tensor $A$ is achieved by sequentially
applying the skeleton decomposition to each of its unfolding matrices.
Initially, the decomposition targets the matrix formed by the multi-indices
$\mathcal{I}_{\leqslant 1}$ and $\mathcal{J}_{>1}$, calculating a matrix
$A\left(s_{<1}s_{1},\mathcal{J}_{>1}\right)\left[A\left(\mathcal{I}_{\leqslant
1},\mathcal{J}_{>1}\right)\right]^{-1}$ which is then unfolded to represent
the first tensor of the tensor-train. The subsequent structure, represented by
$A\left(\mathcal{I}_{\leqslant 1},s_{>1}\right)$, is then unfolded into a
matrix $A\left(\mathcal{I}_{\leqslant 1}s_{2},s_{>2}\right)$ by combining the
index $s_{2}$ with the set $\mathcal{I}_{\leqslant 1}$, and then similarly
decomposed. This process is sequentially applied across all unfolding matrices
of $A$, systematically constructing the tensor-train

$\displaystyle A(s_{1},s_{2},\ldots,s_{n})\approx$ (51) $\displaystyle\approx$
$\displaystyle
A(s_{1},\mathcal{J}_{>1})\left[A(\mathcal{I}_{<2},\mathcal{J}_{>1})\right]^{-1}\times$
$\displaystyle\times$ $\displaystyle
A(\mathcal{I}_{<2},s_{2},\mathcal{J}_{>2})\left[A(\mathcal{I}_{<3},\mathcal{J}_{>2})\right]^{-1}\times\ldots$
$\displaystyle\ldots\times$ $\displaystyle
A(\mathcal{I}_{<k},s_{k},\mathcal{J}_{>k})\left[A(\mathcal{I}_{<k+1},\mathcal{J}_{>k})\right]^{-1}\times\ldots$
$\displaystyle\ldots\times$ $\displaystyle A(\mathcal{I}_{<n},s_{n}).$

The elements $A(\mathcal{I}_{<k},s_{k},\mathcal{J}_{>k})$ are referred to as
“fibers”. These fibers become transform into the cores of the tensor train
following their contraction with the inverse intersection matrices
$\left[A(\mathcal{I}_{<k+1},\mathcal{J}_{>k})\right]^{-1}$. By construction,
the multi-indices $\mathcal{I}_{<k}$ and $\mathcal{J}_{>k}$ are nested,
satisfying the properties

$\mathcal{I}_{k+1}\subset\mathcal{I}_{k}\times
s_{k}\quad\text{and}\quad\mathcal{J}_{k}\subset
s_{k+1}\times\mathcal{J}_{k+1}.$ (52)

In that case, the tensor-train approximation preserves the elements of $A$
that are included in the fibers, while it interpolates the remaining elements
[39].

This approach is inefficient, as it acts on the explicit tensor representation
of $A$. In practice, TCI provides an exponential improvement by acting on a
black-box representation of $A$. It begins with an initial tensor-train
approximation of rank-$r$, as in (51), that includes a collection of nested
multi-indices $\mathcal{J}_{k}$. Then, it starts a recursive procedure. The
first fiber of the approximation, $A(s_{1},\mathcal{J}_{>1})$, is sampled from
the black-box function and decomposed. This produces a set of multi-indices
$\mathcal{I}_{<2}$, which facilitates the sampling of the intersection
submatrix $A(\mathcal{I}_{<2},\mathcal{J}_{>1})$ and the subsequent fiber
$A(\mathcal{I}_{<2},s_{2},\mathcal{J}_{>2})$. This process is recursively
continued through the final tensor of the approximation, after which the
procedure reverses from the last tensor back to the first. Each of these
“sweeps” updates every tensor core of the approximation, repeating in a left-
to-right and right-to-left iterative scheme until a convergence criterion is
satisfied.

This formulation is not rank-adaptive, as it does not adjust the tensor ranks
dinamically. This limitation can be addressed by replacing the standard
skeleton decomposition with its generalization for rectangular matrices [60].
This allows for the progressive expansion of the index sets $\mathcal{I}_{<k}$
and $\mathcal{J}_{>k}$, increasing the ranks or bond dimensions of the
approximation following each sweep. This avoids the need to estimate the ranks
of the solution in advance, and allows to start with a rank-1 approximation
which can be taken at random. The expansion rate of the index sets is
controlled by a hyperparameter. In this work, it is set to perform a rank-1
increase at every decomposition, hence increasing the ranks by two at every
complete sweep.

###  C.3 Extension to MPS/QTT

TCI can be extended to compute MPS/QTT structures by imposing the quantization
of the degrees of freedom of the black-box function. This involves sampling a
different set of fibers $A(\mathcal{I}_{<k},s_{k},\mathcal{J}_{>k})$ from the
sets of MPS/QTT indices. In turn, these indices are related to the tensor-
train indices—that is, the indices labeling the discretization nodes of the
function degrees of freedom—by a simple linear map which arises naturally from
the $k$-nary representation shown in Eq. (3). This map depends on the order of
the MPS tensors, as well as on the number of tensors assigned to each
dimension. Instead of deriving a general expression, we illustrate this map
for a simple example. Consider a bivariate function $f(x_{1},x_{2})$ with
$n_{1}$ and $n_{2}$ qubits per dimension. This function yields the linear maps

$M_{ij}^{\text{serial}}=\begin{pmatrix}2^{n_{1}-1}&0\\\ 2^{n_{1}-2}&0\\\
\vdots&\vdots\\\ 2^{0}&0\\\ 0&2^{n_{2}-1}\\\ 0&2^{n_{2}-2}\\\
\vdots&\vdots\end{pmatrix},\quad
M^{\text{intl}}_{ij}=\begin{pmatrix}2^{n_{1}-1}&0\\\ 0&2^{n_{2}-1}\\\
2^{n_{1}-2}&0\\\ 0&2^{n_{2}-2}\\\ \vdots&\vdots\end{pmatrix}$ (53)

for the serial and interleaved orders respectively. These maps facilitate the
identification and computation of the appropriate fibers for MPS/QTT
decomposition, independent of the specific details of the underlying TCI
implementation.

Crucially, often the pivot matrices $A(\mathcal{I}_{<k+1},\mathcal{J}_{>k})$
of MPS/QTT structures are ill-conditioned for inversion, leading to
significant numerical instabilities. This can be avoided by orthogonalizing
the fibers prior to performing the skeleton decomposition using the QR
decomposition, enhancing the stability and accuracy of the method [23].

##  Appendix D Evaluation of the MPS/QTT error estimation using random
sampling

Figure 17: (a) Error analysis between a product Gaussian distribution and its
SVD approximation. Comparison between the error in $L^{\infty}$ norm of the
full discretization tensor and of randomly sampled subtensors. (b) Error
analysis between a larger product Gaussian distribution computed with SVD and
with iterated Chebyshev expansions. Comparison between the error in norm-2
calculated using the scalar product and using samples. The legends represent
the configurations $(m,n)$ of dimension $m$ and qubits per dimension $n$. The
error bars represent the standard deviation for 100 iterations.

Evaluating the distance in norm $L^{p}$ of a multivariate function with its
approximation requires assessing their exponentially large discretization
tensors. However, it can be estimated following a Monte Carlo approach from a
collection of random samples. This Appendix assesses this estimation for the
product Gaussian distribution up to 200 qubits in two experiments.

The first experiment, shown in Figure 17(a), consists in loading the product
Gaussian function 30 using SVD with a tolerance of $\epsilon=10^{-8}$, from 1
to 6 dimensions and up to 24 qubits. Then, its error in $L^{\infty}$ norm is
estimated from a collection of random samples. Apparently, the estimation
converges monotonically to the error of the full discretization tensor with
the number of samples.

The second experiment, shown in Figure 17(b) consists in loading the product
Gaussian function with two different methods, up to 200 qubits, and comparing
their error. The first method is given by the tensorization of the univariate
function loaded with SVD. The second method is given by an interative MPS
Chebyshev expansion, following Section III. Then, the sampled error in $L^{2}$
norm is compared with the error computed as the MPS norm of their difference,
given by their scalar product. Apparently, the sampled errors converge to the
quasi-exact errors in $L^{2}$ norm given by their scalar product for more than
$10^{4}$ samples.

## References

  * Orús [2014] R. Orús, Annals of Physics 349, 117 (2014).
  * Affleck  _et al._ [1987] I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Phys. Rev. Lett. 59, 799 (1987).
  * Fannes  _et al._ [1989] M. Fannes, B. Nachtergaele, and R. F. Werner, Europhysics Letters 10, 633 (1989).
  * Östlund and Rommer [1995] S. Östlund and S. Rommer, Phys. Rev. Lett. 75, 3537 (1995).
  * Dukelsky  _et al._ [1998] J. Dukelsky, M. A. Martín-Delgado, T. Nishino, and G. Sierra, Europhysics Letters 43, 457 (1998).
  * White [1992] S. R. White, Physical Review Letters 69, 2863 (1992).
  * Oseledets [2011] I. V. Oseledets, SIAM Journal on Scientific Computing 33, 2295 (2011).
  * Khoromskij [2011] B. N. Khoromskij, Constructive Approximation 34, 257 (2011).
  * García-Ripoll [2021] J. J. García-Ripoll, Quantum 5, 431 (2021), arxiv:1909.06619 [quant-ph] .
  * Dolgov  _et al._ [2012a] S. V. Dolgov, B. N. Khoromskij, and I. V. Oseledets, SIAM Journal on Scientific Computing 34, A3016 (2012a).
  * Chen  _et al._ [2023] J. Chen, E. Stoudenmire, and S. R. White, PRX Quantum 4, 040318 (2023).
  * García-Molina  _et al._ [2022] P. García-Molina, J. Rodríguez-Mediavilla, and J. J. García-Ripoll, Physical Review A 105, 012433 (2022), arxiv:2104.02668 [quant-ph] .
  * Lin and White [2023] L. Lin and S. R. White, A Novel Method of Function Extrapolation Inspired by Techniques in Low-entangled Many-body Physics (2023), arxiv:2308.09001 [quant-ph] .
  * Dolgov  _et al._ [2012b] S. Dolgov, B. Khoromskij, and D. Savostyanov, Journal of Fourier Analysis and Applications 18, 915 (2012b).
  * Gourianov  _et al._ [2022] N. Gourianov, M. Lubasch, S. Dolgov, Q. Y. van den Berg, H. Babaee, P. Givi, M. Kiffner, and D. Jaksch, Nature Computational Science 2, 30 (2022), arxiv:2106.05782 [physics, physics:quant-ph] .
  * García-Molina  _et al._ [2023] P. García-Molina, L. Tagliacozzo, and J. J. García-Ripoll, Global optimization of MPS in quantum-inspired numerical analysis (2023), arxiv:2303.09430 [quant-ph] .
  * Ye and Loureiro [2024] E. Ye and N. Loureiro, Quantized tensor networks for solving the Vlasov-Maxwell equations (2024), arxiv:2311.07756 [physics, physics:quant-ph] .
  * Khoromskij and Veit [2016] B. Khoromskij and A. Veit, Computational Methods in Applied Mathematics 16, 145 (2016).
  * Peng  _et al._ [2023] R. Peng, J. Gray, and G. K.-L. Chan, Physical Review Research 5, 013156 (2023).
  * Núñez Fernández  _et al._ [2022] Y. Núñez Fernández, M. Jeannin, P. T. Dumitrescu, T. Kloss, J. Kaye, O. Parcollet, and X. Waintal, Physical Review X 12, 041018 (2022).
  * Schön  _et al._ [2005] C. Schön, E. Solano, F. Verstraete, J. I. Cirac, and M. M. Wolf, Phys. Rev. Lett. 95, 110503 (2005).
  * Melnikov  _et al._ [2023] A. A. Melnikov, A. A. Termanova, S. V. Dolgov, F. Neukart, and M. R. Perelshtein, Quantum Science and Technology 8, 035027 (2023).
  * Oseledets and Tyrtyshnikov [2010] I. Oseledets and E. Tyrtyshnikov, Linear Algebra and its Applications 432, 70 (2010).
  * Hur  _et al._ [2023] Y. Hur, J. G. Hoskins, M. Lindsey, E. Stoudenmire, and Y. Khoo, Applied and Computational Harmonic Analysis 67, 101575 (2023).
  * Lindsey [2023] M. Lindsey, Multiscale interpolative construction of quantized tensor trains (2023), arxiv:2311.12554 [cs, math] .
  * Chen and Lindsey [2024] J. Chen and M. Lindsey, Direct interpolative construction of the discrete Fourier transform as a matrix product operator (2024), arXiv:2404.03182 [quant-ph].
  * Mason and Handscomb [2003] J. C. Mason and D. C. Handscomb,  _Chebyshev Polynomials_ (CRC Press, Boca Raton; London New York, 2003).
  * Trefethen [2019] L. N. Trefethen,  _Approximation Theory and Approximation Practice, Extended Edition_ (Society for Industrial and Applied Mathematics, Philadelphia, PA, 2019).
  * Bañuls  _et al._ [2020] M. C. Bañuls, D. A. Huse, and J. I. Cirac, Phys. Rev. B 101, 144305 (2020).
  * Lubasch  _et al._ [2018] M. Lubasch, P. Moinier, and D. Jaksch, Journal of Computational Physics 372, 587 (2018).
  * Yang  _et al._ [2020] Y. Yang, S. Iblisdir, J. I. Cirac, and M. C. Bañuls, Phys. Rev. Lett. 124, 100602 (2020).
  * Cirac  _et al._ [2021] J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Rev. Mod. Phys. 93, 045003 (2021).
  * Iblisdir  _et al._ [2007] S. Iblisdir, J. I. Latorre, and R. Orús, Phys. Rev. Lett. 98, 060402 (2007).
  * Gorodetsky  _et al._ [2019] A. Gorodetsky, S. Karaman, and Y. Marzouk, Computer methods in applied mechanics and engineering 347, 59 (2019).
  * Dolgov  _et al._ [2012c] S. V. Dolgov, B. N. Khoromskij, and I. V. Oseledets, SIAM Journal on Scientific Computing 34, A3016 (2012c), ttps://doi.org/10.1137/120864210 .
  * Latorre [2005] J. I. Latorre, Image compression and entanglement (2005), arXiv:quant-ph/0510031 [quant-ph] .
  * Shi and Townsend [2021] T. Shi and A. Townsend, SIAM Journal on Matrix Analysis and Applications 42, 275 (2021).
  * Jobst  _et al._ [2023] B. Jobst, K. Shen, C. A. Riofrío, E. Shishenina, and F. Pollmann, arXiv preprint arXiv:2311.07666 (2023).
  * Dolgov and Savostyanov [2020] S. Dolgov and D. Savostyanov, Computer Physics Communications 246, 106869 (2020), arxiv:1903.11554 [math] .
  * Savostyanov and Oseledets [2011] D. Savostyanov and I. Oseledets, in  _The 2011 International Workshop on Multidimensional (nD) Systems_ (IEEE, Poitiers, France, 2011) pp. 1–8.
  * Savostyanov [2014] D. V. Savostyanov, Linear Algebra and its Applications 458, 217 (2014), arxiv:1305.1818 [math] .
  * Ritter  _et al._ [2024] M. K. Ritter, Y. Núñez Fernández, M. Wallerberger, J. von Delft, H. Shinaoka, and X. Waintal, Physical Review Letters 132, 056501 (2024).
  * Oseledets [2013] I. V. Oseledets, Constructive Approximation 37, 1 (2013).
  * Boyd [1989] J. P. Boyd,  _Chebyshev & Fourier Spectral Methods_, edited by C. A. Brebbia, S. A. Orszag, J. Argyris, K.-J. Bathe, A. S. Cakmak, J. Connor, R. McCrory, C. S. Desai, K.-P. Holz, F. A. Leckie, G. Pinder, A. R. S. Ponter, J. H. Seinfeld, P. Silvester, P. Spanos, W. Wunderlich, and S. Yip, Lecture Notes in Engineering, Vol. 49 (Springer Berlin Heidelberg, Berlin, Heidelberg, 1989).
  * Clenshaw [1955] C. W. Clenshaw, Mathematics of Computation 9, 118 (1955).
  * Grant [1970] J. A. Grant, The Mathematical Gazette 54, 96 (1970).
  * Driscoll and Fornberg [2001] T. A. Driscoll and B. Fornberg, Numerical Algorithms 26, 77 (2001).
  * Waldvogel [2006] J. Waldvogel, BIT Numerical Mathematics 46, 195 (2006).
  * Golub and Welsch [1969] G. H. Golub and J. H. Welsch, Mathematics of computation 23, 221 (1969).
  * García-Ripoll [2024] J. J. García-Ripoll, Minimal matrix-product state algorithms library., https://github.com/juanjosegarciaripoll/seemps2 (2024).
  * Rodríguez-Aldavero [2024] J. J. Rodríguez-Aldavero, Simulations for the figures in the paper., https://github.com/jjrodriguezaldavero/seemps_chebyshev_approximation (2024).
  * Vysotsky  _et al._ [2021] L. I. Vysotsky, A. V. Smirnov, and E. E. Tyrtyshnikov, Lobachevskii Journal of Mathematics 42, 1608 (2021), arxiv:2103.12129 [hep-ph] .
  * Alexandrov  _et al._ [2023] B. Alexandrov, G. Manzini, E. W. Skau, P. M. D. Truong, and R. G. Vuchov, Mathematics 11, 534 (2023).
  * Erpenbeck  _et al._ [2023] A. Erpenbeck, W.-T. Lin, T. Blommel, L. Zhang, S. Iskakov, L. Bernheimer, Y. Núñez-Fernández, G. Cohen, O. Parcollet, X. Waintal, and E. Gull, Physical Review B 107, 245135 (2023), arxiv:2303.11199 [cond-mat] .
  * Chertkov  _et al._ [2022] A. Chertkov, G. Ryzhakov, G. Novikov, and I. Oseledets, Optimization of Functions Given in the Tensor Train Format (2022), arxiv:2209.14808 [cs, math] .
  * Sozykin  _et al._ [2022] K. Sozykin, A. Chertkov, R. Schutski, A.-H. Phan, A. S. Cichocki, and I. Oseledets, Advances in Neural Information Processing Systems 35, 26052 (2022).
  * Dolgov  _et al._ [2020] S. Dolgov, K. Anaya-Izquierdo, C. Fox, and R. Scheichl, Statistics and Computing 30, 603 (2020).
  * Lidiak  _et al._ [2022] A. Lidiak, C. Jameson, Z. Qin, G. Tang, M. B. Wakin, Z. Zhu, and Z. Gong, Quantum state tomography with tensor train cross approximation (2022), arXiv:2207.06397 [quant-ph].
  * Goreinov  _et al._ [2010] S. A. Goreinov, I. V. Oseledets, D. V. Savostyanov, E. E. Tyrtyshnikov, and N. L. Zamarashkin, How to Find a Good Submatrix, in  _Matrix Methods: Theory, Algorithms and Applications_ (WORLD SCIENTIFIC, 2010) pp. 247–256.
  * Mikhalev and Oseledets [2018] A. Mikhalev and I. V. Oseledets, Linear Algebra and its Applications 538, 187 (2018), arxiv:1502.07838 [cs, math] .
  * Chertkov  _et al._ [2023] A. Chertkov, G. Ryzhakov, and I. Oseledets, SIAM Journal on Scientific Computing 45, A2101 (2023).
  * Bartholdi [1982] J. J. Bartholdi, Operations Research Letters 1, 190 (1982).