# A hybrid classical-quantum approach to solve the heat equation using quantum
annealers

Giovani G. Pollachini  QuanBy Computação Quântica, Florianópolis, Brazil.
Juan P. L. C. Salazar  Engenharia Aeroespacial, Universidade Federal de Santa
Catarina, Joinville, SC, 89.219-600, Brazil  Caio B. D. Góes  Departamento de
Física, Universidade Federal de Santa Catarina, Florianópolis, SC, 88040-900,
Brazil  Thiago O. Maciel  Departamento de Física, Universidade Federal de
Santa Catarina, Florianópolis, SC, 88040-900, Brazil  Eduardo I. Duzzioni
QuanBy Computação Quântica, Florianópolis, Brazil.  Departamento de Física,
Universidade Federal de Santa Catarina, Florianópolis, SC, 88040-900, Brazil
<EMAIL_ADDRESS>
###### Abstract

The numerical solution of partial differential equations by discretization
techniques is ubiquitous in computational physics. In this work we benchmark
this approach in the quantum realm by solving the heat equation for a square
plate subject to fixed temperatures at the edges and random heat sources and
sinks within the domain. The hybrid classical-quantum approach consists in the
solution on a quantum computer of the coupled linear system of equations that
result from the discretization step. Owing to the limitations in the number of
qubits and their connectivity, we use the Gauss-Seidel method to divide the
full system of linear equations into subsystems, which are solved iteratively
in block fashion. Each of the linear subsystems were solved using 2000Q and
Advantage quantum computers developed by D-Wave Systems Inc. By comparing
classical numerical and quantum solutions, we observe that the errors and
chain break fraction are, on average, greater on the 2000Q system. Unlike the
classical Gauss-Seidel method, the errors of the quantum solutions level off
after a few iterations of our algorithm. This is partly a result of the span
of the real number line available from the mapping of the chosen size of the
set of qubit states. We verified this by using techniques to progressively
shrink the range mapped by the set of qubit states at each iteration
(increasing floating-point accuracy). As a result, no leveling off is
observed. However, an increase in qubits does not translate to an overall
lower error. This is believed to be indicative of the increasing length of
chains required for the mapping to real numbers and the ensuing limitations of
hardware.

##  I Introduction

In fields of computational physics, the mathematical description of the
problem of interest relies on the solution of partial differential equations
restricted to boundary conditions and initial conditions. Most of the
techniques are based on discretization methods, where the physical domain is
divided into cells or volumes. When implicit methods are chosen, the discrete
set of coupled equations to be solved are cast as a system of linear equations
of the form $A\mathbf{x}=\mathbf{b}$, where $A$ is generally a $N\times N$
sparse matrix and the vector $\mathbf{b}$ is known, while the vector
$\mathbf{x}$ is the sought solution to the problem. The largest fraction of
computational time is spent on the solution of the linear systems.

After the advent of the Harrow, Hassidim, and Lloyd (HHL) algorithm [1] to
solve the linear quantum problem, subsequent work has used it as subroutine to
solve ordinary or partial differential equations [2, 3, 4, 5, 6, 7, 8, 9].
However, the solutions for the differential equations were presented in form
of a vector state, restricted to the calculation of global properties of the
system of interest, otherwise the cost to obtain each amplitude of the quantum
state would render the quantum approach too expensive. Clader et al. [3] argue
that an exponential speedup can be achieved by using the HHL algorithm to
solve the electromagnetic scattering cross-section problem via the finite
element method. However, Montanaro and Pallister [6] showed that the relation
between the size of the system of linear equations and the accuracy of the
solution was missed, which implies that the speedup can be reduced or even
eliminated. A thorough analysis of the resources needed to solve the
electromagnetic scattering cross-section of a 2D target was made by Scherer et
al. [10], concluding that the number of resources must be reduced by many
orders of magnitude for the algorithm to become practical. There are other
issues concerning quantum state preparation and the necessity of an oracle to
encode the coefficients of the matrix $A$. For reasons like these, it is not
yet clear if exponential speedups have been achieved by quantum algorithms
using the linear quantum problem as a subroutine. Recently, Linden et al. [11]
showed that HHL-like algorithms used to solve the heat equation are never
faster than the best classical algorithms. However, a quantum alternative
which uses amplitude amplification can accelerate the process, achieving at
most a quadratic speedup. Although this a very encouraging result, the gate
model of quantum computing needs advance to approach large problems. Following
this reasoning, quantum annealers can be more attractive at the current stage
of quantum computing due to their large number of qubits and increase
connectivity. For instance, a graph-coloring-based methodology to solve
differential equations through quantum annealers was proposed in [12].

An alternative method for solving the linear system of equations by quantum
annealing was presented in Refs. [13] and Rogers and Singleton Jr [14]. An
analysis of the resources employed for this task was made in Ref. [15],
showing that in some circumstances it is possible to have advantage when
compared to traditional algorithms to solve least squares problems. The linear
system is encoded in a QUBO (quadratic unconstrained binary optimization)
problem [16] through a well-defined protocol, whose ground state is the
solution to our problem of interest. As opposed to a HHL-type algorithm, the
final solution of the linear system of equations is a classical state
codifying the vector $\mathbf{x}$. In this work we solve the steady-state 2D
heat equation by a hybrid classical-quantum algorithm. Applying the finite
difference method to the partial differential equation yields a system of
coupled linear equations. In order to enable the solution of larger linear
systems of equations, we use the block Gauss-Seidel method and solve the full
system iteratively, where the smaller subsystems of equations are solved in
the quantum processing units (QPU). The quantum algorithm is implemented in
the 2000Q and Advantage systems of D-Wave Systems Inc.

The paper is divided as follows. In §II we introduce the D-Wave quantum
annealing methodology, review the approach of Rogers and Singleton Jr [14] to
solve the linear system of equations, present the block Gauss-Seidel method,
and elaborate about quantum resources required to implement the algorithm. The
results are presented in §III and are followed by conclusions and perspectives
in §IV.

##  II Methods

As described above, the basic idea of our hybrid approach is to replace a
classical subroutine to solve the system of linear equations that result from
the discretization of partial differential equations by a quantum algorithm.
In the proposed method the system of linear equations is solved using quantum
annealing, where the final Hamiltonian encodes a QUBO problem [14]. In the
following we briefly review the quantum annealing followed by the QUBO
problem.

###  II.1 D-Wave quantum annealing

D-Wave quantum annealing is based on the adiabatic quantum evolution of the
transverse Ising Hamiltonian with the system in thermal equilibrium at a very
low temperature [17]. The goal of the computation is to minimize an objective
function (energy) for which the solution to the desired problem is
approximately encoded in its ground state [18]. In order to achieve the ground
state of the system with a high success rate, the quantum fluctuations must be
dominant over thermal fluctuations [17].

The idea behind this model is the preparation of the system as the ground
state of an easy-to-prepare Hamiltonian, for instance ($\hbar\omega=1$),

$\displaystyle H(0)=-\sum_{i=1}^{n}\sigma_{x}^{i},$ (1)

where $\sigma_{x}^{i}$ is the Pauli matrix in the $x$ direction for the
$i$-particle and its ground state is $|\psi(0)\rangle=|+\rangle^{\otimes n}$,
which contains the superposition of all states in the computational basis
equally distributed. By choosing the final Hamiltonian $H(1)$ as the operator
whose ground state contains the solution to the desired problem, it is
possible to interpolate between these two Hamiltonians with a schedule
function $f(s)$, subject to $f(0)=0$ and $f(1)=1$, resulting in,

$\displaystyle H(s)=\left(1-f(s)\right)H(0)+f(s)H(1),$ (2)

where $s=t/T$ ($0\leq s\leq 1$) is the dimensionless time, $t$ is the current
time and $T$ is the total time of the computation. The final Hamiltonian
implemented by D-Wave is the Ising Hamiltonian,

$\displaystyle
H(1)=H_{Ising}=\sum_{i}\alpha_{i}\sigma_{z}^{i}+\sum_{i,j}\beta_{ij}\sigma_{z}^{i}\sigma_{z}^{j}\
,$ (3)

with the parameters $\alpha_{i}$ and $\beta_{ij}$ being chosen to encode the
solution to the problem and $\sigma_{z}^{i}$ is the Pauli matrix in the $z$
direction for the particle $i$.

According to the adiabatic theorem [19], if the evolution is performed
adiabatically, then the ground state of the problem Hamiltonian $H(1)$ is
found with high probability, provided that the evolution time $T$ from the
initial Hamiltonian $H(0)$ to $H(1)$ is proportional to $g_{min}^{-2}$, where
$g_{min}=E_{1}-E_{0}$ is the minimum energy gap between the ground and first
excited states [19, 20]. It is important to observe that a suitable choice for
the schedule function may result in totally different evolution times [21, 20,
22]. Therefore, the final Hamiltonian $H(1)$ must contain the solution to the
linear system of equations [14]. Next, we show how to build $H(1)$, thus
casting the system of linear equations into a QUBO problem.

###  II.2 QUBO algorithm to solve linear systems of equations

O’Malley and Vesselinov [13] and Rogers and Singleton Jr [14] showed how to
map the problem of solving a linear system of equations into a QUBO problem
[16]. Also, in Ref. [23], the linear system of equations can be obtained as a
particular case of a nonlinear system of equations. Given a linear system
$A\mathbf{x}=\mathbf{b}$, the solution $\mathbf{x}$ can be found by minimizing
the function

$\displaystyle
H(\mathbf{x})=(A\mathbf{x}-\mathbf{b})^{{\dagger}}(A\mathbf{x}-\mathbf{b})\ .$
(4)

We can approximate each component $x_{i}$ ($i=1\dots N$) of the vector
$\mathbf{x}$ by its binary representation in a given interval
$[-d_{i},2c_{i}-d_{i})$, with $R$ bits, as:

$\displaystyle x_{i}=c_{i}\sum_{r=0}^{R-1}q_{r}^{i}2^{-r}-d_{i}\ ,$ (5)

where each $q_{r}^{i}\in\\{0,1\\}$ is the $r$-th bit in the binary
representation of $x_{i}$. If $d_{i}>0$ and $c_{i}>d_{i}/2$, the domain of
$x_{i}$ contains positive and negative real values. This can easily be
generalized to complex variables [14]. After substituting (5) into (4) and
dropping the positive constant term, the equation becomes,

$\displaystyle
H(q)=\sum_{r=0}^{R-1}\sum_{i=0}^{N-1}a_{r}^{i}q_{r}^{i}+\sum_{r,s=0}^{R-1}\sum_{i,j=0}^{N-1}b_{rs}^{ij}q_{r}^{i}q_{s}^{j}\
,$ (6)

where the coefficients are given by,

$\displaystyle a_{r}^{i}$
$\displaystyle=-2\left(\sum_{j,k=0}^{N-1}A_{ki}A_{kj}c_{i}d_{j}+\sum_{j=0}^{N-1}A_{ji}c_{i}b_{j}\right)2^{-r}\
,$ (7) $\displaystyle b_{rs}^{ij}$
$\displaystyle=\left(\sum_{k=0}^{N-1}A_{ki}A_{kj}c_{i}c_{j}\right)2^{-(r+s)}.$
(8)

Here, $A_{ij}$ are coefficients of the matrix $A$ and $b_{i}$ are components
of the vector $\mathbf{b}$. It is important not to confuse $b_{i}$ with
$b^{ij}_{rs}$.

In quantum annealers, the physical qubits are indexed by a 1-dimensional
linear index $\ell$ (see Eq. (3)), while the logical qubits are indexed by
2-dimensional indices $i$ and $r$, where $i=0,1,...,N-1$ and $r=0,1,...,R-1$,
as shown in Eq. (6). The mapping between the logical an physical indices is
given by,

$\displaystyle\ell(i,r)$ $\displaystyle=i.R+r,\,\,\ell=0,...,NR-1,\,\,$ (9)

such that the inverse mapping is,

$\displaystyle i_{\ell}$ $\displaystyle=\lfloor\ell/R\rfloor\ ,$
$\displaystyle r_{\ell}$ $\displaystyle=\ell\mod R.$ (10)

The mapping from $q_{\ell}\rightarrow\sigma_{\ell}$ is
$\sigma_{\ell}=2q_{\ell}-1$ with $\sigma_{\ell}=\\{+1,-1\\}$.

Finally, the quantization of the QUBO problem is made through
$\hat{q}_{\ell}|q\rangle=q_{\ell}|q\rangle$, where the operators are such that
$\hat{q}^{2}_{\ell}=\hat{q}_{\ell}$, with eigenvalues $q_{\ell}\in\\{0,1\\}$
and $|q\rangle=|q_{0}\rangle\otimes...\otimes|q_{NR-1}\rangle$. In order to
solve the system of linear equations, the coefficients $a^{i}_{r}\rightarrow
a_{\ell}$ and $b^{ij}_{rs}\rightarrow b_{\ell\kappa}$, described in Eqs. (7)
and (8), must be provided to the D-Wave Ocean Software [24], the interface
that communicates with the quantum hardware.

The readout of the quantum algorithm is a bit string $q_{0}q_{1}...q_{NR-1}$
of $N.R$ qbits that returns to the classical algorithm for further processing.
In principle, the qubits must be fully connected. However, in situations where
the matrix $A$ is sparse, several coefficients $b_{rs}^{ij}$ are set to zero,
hence not all connections between the qubits are necessary. With knowledge of
the above, if the full linear system of equations is solved directly, the
number of required qubits $N.R$ for equivalent floating-point accuracy of
classical algorithms makes this implementation prohibitive at the current
stage of technological development. However, the linear system can be divided
into subsystems, and the full linear system solved iteratively, in a hybrid
quantum-classical approach, as detailed next.

###  II.3 Iterative QUBO method for solving linear systems

The block Gauss-Seidel method consists in solving the linear system
$A\mathbf{x}=\mathbf{b}$ iteratively, by applying the Gauss-Seidel method on
the original system divided in blocks [25]. Suppose, for simplicity, that the
size of the system is an even number $N$ and divide it into two subsystems,

$\displaystyle\bmqty{undef}\bmqty{undef}=\bmqty{undef}\ .$ (11)

This translates to

$\displaystyle A_{11}\mathbf{x}_{1}+A_{12}\mathbf{x}_{2}$
$\displaystyle=\mathbf{b}_{1},\ $ (12) $\displaystyle
A_{21}\mathbf{x}_{1}+A_{22}\mathbf{x}_{2}$ $\displaystyle=\mathbf{b}_{2}\ .$
(13)

In the iterative method we make an initial guess for $\mathbf{x}_{2}$, such as
$\mathbf{x}_{2}^{(0)}=0$, where the upper index of the block-vector
$\mathbf{x}_{2}^{(0)}$ defines the step of the iteration process.
Rearrangement of equations (12) and (13) allow us to set up an iterative
process where the initial guess is updated by solving two linear systems of
size $N/2$.

$\displaystyle A_{11}\mathbf{x}_{1}^{(1)}$
$\displaystyle=\mathbf{b}_{1}-A_{12}\mathbf{x}_{2}^{(0)}\ ,$ (14)
$\displaystyle A_{22}\mathbf{x}_{2}^{(1)}$
$\displaystyle=\mathbf{b}_{2}-A_{21}\mathbf{x}_{1}^{(1)}\ .$ (15)

At the $k$-th step,

$\displaystyle A_{11}\mathbf{x}_{1}^{(k)}$
$\displaystyle=\mathbf{b}_{1}-A_{12}\mathbf{x}_{2}^{(k-1)}\ ,$ (16)
$\displaystyle A_{22}\mathbf{x}_{2}^{(k)}$
$\displaystyle=\mathbf{b}_{2}-A_{21}\mathbf{x}_{1}^{(k)}\ .$ (17)

If the sequences $\mathbf{x}_{1}^{(k)}$ and $\mathbf{x}_{2}^{(k)}$ converge,
then $\mathbf{x}_{1}^{(k)}\simeq\mathbf{x}_{1}^{(k-1)}$,
$\mathbf{x}_{2}^{(k)}\simeq\mathbf{x}_{2}^{(k-1)}$ and the original linear
system is satisfied. Supposing $A_{11}$ and $A_{22}$ are invertible, it can be
shown that,

$\displaystyle\mathbf{x}_{1}^{(k)}-\mathbf{x}_{1}$
$\displaystyle=\left(A_{22}^{\,-1}A_{21}A_{11}^{\,-1}A_{12}\right)^{k}\left(\mathbf{x}_{1}^{(0)}-\mathbf{x}_{1}\right)\
,$ (18) $\displaystyle\mathbf{x}_{2}^{(k)}-\mathbf{x}_{2}$
$\displaystyle=\left(A_{11}^{\,-1}A_{12}A_{22}^{\,-1}A_{21}\right)^{k}\left(\mathbf{x}_{2}^{(0)}-\mathbf{x}_{2}\right)\
.$ (19)

Hence, a sufficient condition for the convergence of this method is:

$\displaystyle\norm{A_{22}^{\,-1}A_{21}A_{11}^{\,-1}A_{12}}$ $\displaystyle<1\
,$ (20) $\displaystyle\norm{A_{11}^{\,-1}A_{12}A_{22}^{\,-1}A_{21}}$
$\displaystyle<1\ .$ (21)

The generalization of the method is possible by dividing the original linear
system into several blocks whose size is adjustable to the particular
interest. If the system is divided into $D$ blocks, each smaller system has
size $N/D\times N/D$. The number of required iterations is not known
beforehand, but defined by a convergence criterion. A typical criterion
depends on the normalized residual $r^{(k)}$, defined by

$\displaystyle
r^{(k)}\equiv\frac{\norm{A\mathbf{x}^{(k)}-\mathbf{b}}}{\norm{\mathbf{b}}}$
(22)

where $\mathbf{x}^{(k)}$ is the solution of the system of linear equations
after $k$ iteration steps. A solution is considered converged when $r^{(k)}$
meets a specified tolerance $\tau$, such that $r^{(k)}\leq\tau$. Here, we
instead monitor the relative error of the solution vector at the $k$-th
iteration $\mathbf{x}^{(k)}$, given as,

$\displaystyle
e^{(k)}\equiv\frac{\norm{\mathbf{x}^{(k)}-\mathbf{x}}}{\norm{\mathbf{x}}}\,.$
(23)

It can be shown that $e^{(k)}$ and $r^{(k)}$ are related according to,

$\displaystyle e^{(k)}\leq\kappa(A)r^{(k)}\,,$ (24)

where $\kappa(A)$ is the condition number of the coefficient matrix $A$.
Preconditioning techniques should be used when $A$ is ill-conditioned. In this
work we do not pursue this issue further, since it would be relevant for both
classical and quantum algorithms.

#### Scaling of the method

The representation of a vector $\mathbf{x}$ with $N$ components requires
$N\cdot R$ qubits, where $R$ is the number of qubits that gives the desired
accuracy of floating-point numerical representation, $\epsilon$, of each
variable $x_{i}$. However, as stated above, it is possible to solve smaller
subsystems of size $N/D$ by partitioning the full system in $D$ parts. In this
case the number of qubits scales as $(N/D)\cdot R$. In order to reduce the
number of qubits required in the computation for a given floating-point
accuracy of $x_{i}$, one possibility is to shrink the interval
$[-d_{i},2c_{i}-d_{i}]$ at each iteration, which we adopt here. Any reduction
in the number of required qubits is significant in the era of noisy
intermediate-scale quantum (NISQ) computers, where at most few thousands of
qubits are available. Even more dramatic is the number of connections
necessary to solve the QUBO problems. In the worst case, every logical qubits
must be fully connected (see the coefficient $b_{rs}^{ij}$ in Eq. (7)). The
block Gauss-Seidel method can reduce the number of connections among qubits by
a factor $N^{2}R^{2}\left(1-1/D^{2}\right)$. Additionally, when $A$ is sparse,
connectivity is reduced considerably.

If $\epsilon$ is a bound for the accuracy of the floating-point numerical
representation of a physical variable $x_{i}$, for a solution represented with
$R$ qubits, then

$\displaystyle\frac{2c_{i}}{2^{R}}\leq\epsilon\ ,$

since $2c_{i}$ is the length of the open interval $[-d_{i},2c_{i}-d_{i})$
where the solution is represented and $2^{R}$ is the number of different
points that can be represented with $R$ bits inside the aforementioned
interval. From the expression above, the number of qubits required to
represent one variable with floating-point accuracy $\epsilon$ scales as

$\displaystyle R\geq\log_{2}\left(2c_{i}/\epsilon\right)\ .$

Since the total computing time $\mathcal{T}$ depends on the energy gap between
the ground and first excited states of the system, it is hard to obtain
precise estimates of $\mathcal{T}$. As noted in [26], the task of analytically
extracting the minimum gap scaling has been extremely difficult in practice,
except for a few special cases, such as for Deutsch-Jozsa [27] and Grover
algorithms [21, 28].

In Table 1 we summarize the consuming quantum resources used in the iterative
and non-iterative QUBO algorithms to solve the linear system of equations.

Method |  # qubits | time
---|---|---
Original QUBO | $N\cdot R$ | $\mathcal{T}_{N,R}$
Iterative QUBO | $(N/D)\cdot R$ | $\mathcal{T}_{N/D,R}\cdot D\cdot n_{\text{iter}}$
Table 1: Comparison between original QUBO and iterative QUBO methods.
$\mathcal{T}_{N,R}$ indicates the time for running one QUBO algorithm for a
linear system of size $N$ and $R$ qubits for representing the solution, while
$\mathcal{T}_{N/D,R}$ is the time for solving a system of size $N/D$.
$n_{\text{iter}}$ is the number of iterations in the case of the iterative
QUBO method. Both times $\mathcal{T}_{N,R}$ and $\mathcal{T}_{N/D,R}$ depend
on the annealing time and the number of runs to sampled solutions of the
problem.

##  III Results and Discussion: Two Dimensional Heat Equation

Figure 1:  Temperature distribution in the square plate. Figure 2:  The
relative error as a function of the number of iterations of the block Gauss-
Seidel method for different numbers of qubits used in the numerical
representation of the solution $\mathbf{x}$. The solutions obtained by the
2000Q (red) and Advantage (green) systems are compared with the classic Gauss-
Seidel implementation (blue). Figure 3:  Average chain break fraction as a
function of the number of iterations in the block Gauss-Seidel method for
different numbers of precision qubits. Figure 4:  The relative error as a
function of the number of iterations of the block Gauss-Seidel method, with a
shrink factor of $0.8$ of the search domain at each iteration, for different
numbers of qubits used in the numerical representation of the solution
$\mathbf{x}$. The solutions obtained by the 2000Q (red) and Advantage (green)
systems are compared with the classic Gauss-Seidel implementation (blue).

In order to assess the viability and performance of the proposed method, we
solve the heat equation for a square plate of length $L$ subject to fixed
temperature at its edges and multiple heat sources and sinks within the
domain. The steady state solution $T(x,y)$ is then determined with the
following PDE and boundary conditions,

$\displaystyle\dfrac{\partial^{2}T}{\partial
x^{2}}+\dfrac{\partial^{2}T}{\partial y^{2}}=0\,\,\,,$ (25) $\displaystyle
T(x,0)$ $\displaystyle=0^{\circ}C\ ,\ \ \quad\quad\ \ x\in[0,L]$
$\displaystyle T(x,L)$ $\displaystyle=\frac{x}{L}\cdot 100^{\circ}C\ ,\ \
x\in[0,L]$ $\displaystyle T(0,y)$ $\displaystyle=0^{\circ}C\ ,\ \ \quad\quad\
\ y\in[0,L]$ $\displaystyle T(L,y)$ $\displaystyle=\frac{y}{L}\cdot
100^{\circ}C\ ,\ \ y\in[0,L]\ .$ (26)

The problem is solved numerically by the finite difference method. The domain
is divided into $m$ segments of equal length, with nodes located at positions
$x_{i}=(L/m)\cdot i$ and $y_{j}=(L/m)\cdot j$, where $i,j=0\dots m$. The
points $(x_{i},y_{j})$ of the grid corresponding to $i=0$, $i=m$, $j=0$ or
$j=m$ form the boundary of the plate and are determined by the boundary
conditions (III). The points inside the grid are the unknowns, and the partial
derivatives can be written as

$\displaystyle\begin{split}\frac{\partial^{2}T}{\partial
x^{2}}(x_{i},y_{j})=&\frac{T(x_{i+1},y_{j})-2T(x_{i},y_{j})+T(x_{i-1},y_{j})}{(L/m)^{2}}\\\
+&\frac{(L/m)^{2}}{12}\dfrac{\partial^{4}T}{\partial
x^{4}}(x_{i}^{\ast},y_{j})\,\,,\\\ \frac{\partial^{2}T}{\partial
y^{2}}(x_{i},y_{j})=&\frac{T(x_{i},y_{j+1})-2T(x_{i},y_{j})+T(x_{i},y_{j-1})}{(L/m)^{2}}\\\
+&\frac{(L/m)^{2}}{12}\dfrac{\partial^{4}T}{\partial
y^{4}}(x_{i},y_{j}^{\ast})\,\,\,,\end{split}$ (27)

for some $x_{i}^{\ast}\in[x_{i-1},x_{i+1}]$ and
$y_{j}^{\ast}\in[y_{j-1},y_{j+1}]$. Ignoring terms proportional to $(L/m)^{2}$
or higher orders, the PDE becomes,

$\displaystyle\begin{split}T(x_{i+1},y_{j})&-4T(x_{i},y_{j})+T(x_{i-1},y_{j})\\\
&+T(x_{i},y_{j+1})+T(x_{i},y_{j-1})=0\,,\end{split}$ (28)

where $i,j=1,2,\ldots,m-1$. Equation (28) can be written in the form of a
linear system $A\textbf{x}=\mathbf{b}$, where $A$ is a sparse matrix with few
non-zero entries and $\mathbf{b}$ is defined by the boundary conditions and
eventual source terms.

The spatial discretization forms a $9\times 9$ matrix ($m=11$), resulting in a
system of coupled linear equations with $81$ unknown variables, which is
solved by the iterative QUBO algorithm on the Advantage and D-Wave 2000Q
systems [29]. A sample solution of the problem is shown in Fig. 1.

This approach enables us to solve larger linear systems of equations,
circumventing the limited number of qubits in the D-Wave quantum processing
units. For sake of comparison, in this work we use at most $7$ precision
qubits, which demands at most $56$ logical qubits to describe all unknown
variables in the block Gauss-Seidel approach, while $567$ would be required to
solve the full system of linear equations non-iteratively. We will see that
the number of logical qubits is much smaller than the number of physical
qubits in both D-Wave systems, because of the limited connectivity between
physical qubits.

In Fig. 2 the dependence of the relative error defined in Eq. (23) is shown as
a function of the number of iterations of the block Gauss-Seidel method, for
different numbers of qubits used in the numerical representation of the
unknown variables. In this example, the smallest error is found using three
qubits in the Advantage system. We also notice that the error does not
continually decrease, as would be expected in a classical algorithm. This
occurs in part because of the limited floating-point accuracy of the variables
$x_{i}$, such that after a certain number of iterations, an improvement of the
solution is not possible. Another observation is that the increase in qubits
does not necessarily translate into smaller errors. This result was not
expected. However, it can be related to the appearance of more noise in the
system, since there are more qubits involved in the numerical representation
of a given variable $x_{i}$, which in turn increases the the probability of
chain breaks. This subject will be discussed in the following. In general, the
Advantage system performs marginally better than 2000Q.

In order to create a logical qubit, generally more than one physical qubit is
necessary. This is a consequence of the limited number of couplers
(connections) between qubits. In the D-Wave quantum processing units, a set of
physical qubits that describe a vertex of a graph form a chain. For this chain
to behave as a basic quantum processing unit, all of its physical qubits must
be in the same state. When this does not occur, the chain is broken and the
state is determined by a majority vote. In Fig. (3) we plot the average chain
break fraction as function of the number of iterations and the number of
qubits used to represent each value of temperature. The average chain break
fraction for each number of precision qubits is greater in the 2000Q than in
the Advantage system and is on average independent on the iteration number.
This result is expected, since in the Chimera topology of the 2000Q system
each qubit is connected to $6$ other qubits, while in the Pegasus topology of
the Advantage system, each qubit connects to $15$ other qubits. Furthermore,
the Advantage system has $5000+$ qubits compared to $2000+$ qubits in the
2000Q system. These features allow for more compact chains that are less
susceptible to noise [30].

A significant improvement in the relative error of Fig. 2 can be obtained by
reducing the length of the interval $2c_{i}$ represented by a given number of
qubits at each iteration. The iterative QUBO algorithm is modified by
shrinking the interval by factor of $\gamma=0.8$ at each iteration, i.e., the
interval for estimating the solution at the $k$-th iteration is given by
$[x_{i}^{k-1}-c_{i}\gamma^{k-1},x_{i}^{k-1}+c_{i}\gamma^{k-1}]$, where
$x_{i}^{k-1}$ is the estimate for $x_{i}$ at the ($k-1$)-th iteration. This
procedure increases the capacity of the QUBO solver to represent the solution
of the system with a fixed number $R$ of qubits, as can be seen in Figure 4.
The implementation of this updated search interval length brings the quantum
solvers on par with the error obtained on classical hardware. It must be
noted, however, that the update of the rhs of equations (16) and (17) during
the Gauss-Seidel iterations is performed on a classical computer, with double
precision, even in the quantum algorithm.

A comparison of Figs. 2 and 4 shows that the 2000Q performs slightly better
than the Advantage system when the update of the search interval length is
incorporated into the algorithm. There is not a clear interpretation why this
occurs, since the Advantage system is a improved version of the 2000Q system,
where the chains are more compact and less susceptible to noise [30]. The
results obtained from both figures also suggest that a greater number of
precision qubits does not necessarily translate to a smaller error, as already
noticed earlier.

###### Acknowledgements.

The authors acknowledge the financial support from the Brazilian funding
agencies Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - CAPES,
Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq, Fundação
de Amparo à Pesquisa e Inovação do Estado de Santa Catarina - FAPESC, and the
Instituto Nacional de Ciência e Tecnologia de Informação Quântica - INCT-IQ.

##  IV  Conclusions and perspectives

The hybrid classical-quantum algorithm we have developed solves partial
differential equations by the finite difference method, yielding a system of
coupled linear equations. This technique was tested with success for the 2D
steady-state heat equation, but in principle can be applied to any partial
differential equation, encompassing a large array of applications for solution
by quantum annealers. The method with updated search interval length allows us
to obtain comparable accuracy to a classical algorithm with double precision
using 3 qubits to numerically represent each variable. The main limitation of
the quantum hardware is the reduced connectivity between qubits, albeit great
strides have been made in order to increase the number of couplers for each
qubit. The iterative algorithm allows for a significant reduction in the
number of logical qubits, from $N.R$ to $N.R/D$, at the expense of an increase
in the number of required iterations for a given accuracy.

It is our understanding that there are good prospects for the iterative QUBO
algorithm and that there is room for further improvements, such as robust
heuristics for update of the search interval length. Another improvement can
be obtained by using an application-oriented embedding, as many linear systems
typical of finite difference and finite volume methods are characterized by
sparseness of the coefficient matrix $A$. This has the advantage of requiring
less couplers than a fully connected problem. Possibly a specific embedding
can be found for this class of problems, allowing for less qubits and enabling
the solution of larger systems of coupled linear equations. Owing to our
limited access to the D-Wave QPUs, we believe that more tests are required to
fully benchmark the performance of 2000Q and Advantage systems for the
solution of systems of linear equations.

## References

  * Harrow et al. [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Physical review letters 103, 150502 (2009).
  * Leyton and Osborne [2008] S. K. Leyton and T. J. Osborne, arXiv preprint arXiv:0812.4423 (2008).
  * Clader et al. [2013] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, Physical review letters 110, 250504 (2013).
  * Cao et al. [2013] Y. Cao, A. Papageorgiou, I. Petras, J. Traub, and S. Kais, New Journal of Physics 15, 013021 (2013).
  * Berry [2014] D. W. Berry, Journal of Physics A: Mathematical and Theoretical 47, 105301 (2014).
  * Montanaro and Pallister [2016] A. Montanaro and S. Pallister, Physical Review A 93, 032324 (2016).
  * Berry et al. [2017] D. W. Berry, A. M. Childs, A. Ostrander, and G. Wang, Communications in Mathematical Physics 356, 1057 (2017).
  * Costa et al. [2019] P. C. Costa, S. Jordan, and A. Ostrander, Physical Review A 99, 012323 (2019).
  * Wang et al. [2020] S. Wang, Z. Wang, W. Li, L. Fan, Z. Wei, and Y. Gu, Quantum Information Processing 19, 1 (2020).
  * Scherer et al. [2017] A. Scherer, B. Valiron, S.-C. Mau, S. Alexander, E. Van den Berg, and T. E. Chapuran, Quantum Information Processing 16, 1 (2017).
  * Linden et al. [2020] N. Linden, A. Montanaro, and C. Shao, arXiv preprint arXiv:2004.06516 (2020).
  * Srivastava and Sundararaghavan [2019] S. Srivastava and V. Sundararaghavan, Physical Review A 99, 052355 (2019).
  * O’Malley and Vesselinov [2016] D. O’Malley and V. V. Vesselinov, in _2016 IEEE High Performance Extreme Computing Conference (HPEC)_ (IEEE, 2016), pp. 1–7.
  * Rogers and Singleton Jr [2020] M. L. Rogers and R. L. Singleton Jr, Frontiers in Physics 8, 265 (2020).
  * Borle and Lomonaco [2019] A. Borle and S. J. Lomonaco, in _International Workshop on Algorithms and Computation_ (Springer, 2019), pp. 289–301.
  * Glover and Kochenberger [2018] F. Glover and G. Kochenberger, arXiv preprint arXiv:1811.11538 (2018).
  * Kadowaki and Nishimori [1998] T. Kadowaki and H. Nishimori, Physical Review E 58, 5355 (1998).
  * Apolloni et al. [1989] B. Apolloni, C. Carvalho, and D. De Falco, Stochastic Processes and their Applications 33, 233 (1989).
  * Jansen et al. [2007] S. Jansen, M.-B. Ruskai, and R. Seiler, Journal of Mathematical Physics 48, 102111 (2007).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Reviews of Modern Physics 90, 015002 (2018).
  * Roland and Cerf [2002] J. Roland and N. J. Cerf, Physical Review A 65, 042308 (2002).
  * An and Lin [2019] D. An and L. Lin, arXiv preprint arXiv:1909.05500 (2019).
  * Chang et al. [2019] C. C. Chang, A. Gambhir, T. S. Humble, and S. Sota, Scientific reports 9, 10258 (2019).
  * Inc [2021] D.-W. S. Inc, _Ocean software documentation_ (2021), eprint https://docs.ocean.dwavesys.com/en/stable/index.html.
  * Golub and Loan [1996] G. H. Golub and C. F. V. Loan, _Matrix computations_ , Johns Hopkins Studies in the Mathematical Sciences (Johns Hopkins University Press, 1996), 3rd ed., ISBN 9780801854149,0801854148,080185413X.
  * Amin and Choi [2009] M. Amin and V. Choi, Physical Review A 80, 062326 (2009).
  * Das et al. [2002] S. Das, R. Kobes, and G. Kunstatter, Physical Review A 65, 062310 (2002).
  * Farhi et al. [2000] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, arXiv preprint quant-ph/0001106 (2000).
  * Inc [2020a] D.-W. S. Inc, _Wave quantum processing unit, tech. rep. (d-wave systems inc, burnaby, bc, canada, 2020) d-wave user manual 09-1109a-v._ (2020a), eprint https://docs.dwavesys.com/docs/latest/doc_qpu.html.
  * Inc [2020b] D.-W. S. Inc, _The d-wave advantage system: An overview, tech. rep. (d-wave systems inc, burnaby, bc, canada, 2020) d-wave technical report series 14-1049a-a._ (2020b), eprint https://dwavesys.com.


# Tractability of sampling recovery on unweighted function classes

David Krieg  Institut für Analysis, Johannes Kepler Universität Linz, Austria.
<EMAIL_ADDRESS>
###### Abstract.

It is well-known that the problem of sampling recovery in the $L_{2}$-norm on
unweighted Korobov spaces (Sobolev spaces with mixed smoothness) as well as
classical smoothness classes such as Hölder classes suffers from the curse of
dimensionality. We show that the problem is tractable for those classes if
they are intersected with the Wiener algebra of functions with summable
Fourier coefficients. In fact, this is a relatively simple implication of
powerful results from the theory of compressed sensing. Tractability is
achieved by the use of non-linear algorithms, while linear algorithms cannot
do the job.

We consider the problem of recovering a high-dimensional function
$f\colon[0,1]^{d}\to\mathbb{C}$ from a class $F_{d}$ using algorithms of the
form

(1) $A_{n}\colon F_{d}\to L_{2}([0,1]^{d}),\quad
A_{n}(f)=\phi(f(x_{1}),\ldots,f(x_{n}))$

with sampling points $x_{j}\in[0,1]^{d}$ and a recovery map
$\phi\colon\mathbb{C}^{n}\to L_{2}$. The error is measured in the $L_{2}$-norm
and in a worst-case setting, i.e.,

${\rm err}(A_{n},F_{d},L_{2})\,:=\,\sup_{f\in F_{d}}\|f-A_{n}(f)\|_{2}.$

It is known that this approximation problem suffers from the _curse of
dimensionality_ for most classical function classes $F_{d}$, including the
smoothness classes of $k$-times continuously differentiable functions [16, 17]
or Sobolev spaces of mixed smoothness [15]. That is, there exist
$C,\gamma,\varepsilon>0$ such that the (information) complexity

(2) $n(\varepsilon,F_{d},L_{2})\,:=\,\min\big{\\{}n\in\mathbb{N}\ |\ \exists
A_{n}\colon{\rm err}(A_{n},F_{d},L_{2})\leq\varepsilon\big{\\}}$

is bounded below by $C(1+\gamma)^{d}$ for all $d\in\mathbb{N}$.

A classical approach to make this problem tractable in high dimensions is to
consider weighted function classes $F_{d}$, assuming a decreasing importance
of the input variables. This approach goes back to [27] and has gained
significant popularity, see also Remark 14 about weighted Korobov spaces.
Unfortunately, it is much harder to provide reasonable unweighted function
classes where the approximation problem is tractable. By _unweighted_ , we
mean that all variables are of equal importance. Formally, one may call a
function class $F_{d}$ unweighted if $f\circ\pi\in F_{d}$ for any $f\in F_{d}$
and any permutation matrix $\pi$.

It was recently discovered by Goda [13] that _numerical integration_ is
polynomially tractable on the unweighted class

(3) $F_{d}^{\log}\,:=\,\bigg{\\{}f\in C([0,1]^{d})\ \Big{|}\
\sum_{k\in\mathbb{Z}^{d}}|\hat{f}(k)|\max(1,\log\|k\|_{\infty})\leq
1\bigg{\\}}.$

Namely, the number of sampling points that are needed to approximate the
integral of any function from $F_{d}^{\log}$ up to an error $\varepsilon$ is
bounded above by a polynomial in $\varepsilon^{-1}$ and $d$. Goda refines a
method of Dick who already showed in [6] that integration is tractable on the
class of $\alpha$-Hölder continuous functions with absolutely convergent
Fourier series. The class $F_{d}^{\log}$ avoids the Hölder condition by
strengthening the condition of absolute convergence in a very slight manner.

Those results raise the question whether the problem of $L_{2}$-approximation
on $F_{d}^{\log}$, whose complexity can only be higher than the complexity of
the integration problem, is polynomially tractable as well. We give a positive
answer.

######  Theorem 1.

There is a constant $c>0$ such that for all $d\in\mathbb{N}$ and
$\varepsilon\in(0,1/2)$, we have

$n(\varepsilon,F_{d}^{\log},L_{2})\,\leq\,c\,d\,\varepsilon^{-3}\log^{3}(\varepsilon^{-1}).$

More generally, the problem of sampling recovery on $F_{d}^{\log}$ in $L_{p}$
is tractable for all $1\leq p<\infty$, see Corollary 13. In fact, we will show
that $L_{p}$-approximation is polynomially tractable for many classical
unweighted classes of functions, like Sobolev spaces of mixed smoothness or
Hölder continuous functions, if it is additionally assumed that the Fourier
series converges absolutely. See Theorem 11 for a general result and
Corollaries 15 and 17 for the corresponding examples. Those findings are
consequences of powerful results from the theory of compressive sensing, as
will be discussed in Section 1.

Interestingly, the curse of dimensionality is only mitigated due to the use of
non-linear sampling algorithms. If one only allows linear algorithms, using
linear recovery maps $\phi$ in (1), the $L_{2}$-approximation problem on
$F_{d}^{\log}$ suffers from the curse of dimensionality, see Lemma 18. This
fact also implies that the sampling recovery on $F_{d}^{\log}$ in $L_{\infty}$
is intractable, see Corollary 19.

######  Remark 2.

For numerical integration on the class $F_{d}^{\log}$, Goda [13] obtained an
upper bound of order $d^{3}\varepsilon^{-3}$. Since the upper bound from
Theorem 1 also holds for numerical integration, we improve upon this bound if
the dimension $d$ is larger than $\log^{3/2}(\varepsilon^{-1})$. While the
present paper was under review, Goda’s upper bound on the complexity of
integration has been further improved to $d\varepsilon^{-3}$, see [3].

######  Remark 3.

The function class $F_{d}^{\log}$ defined in (3) is from the first version of
the paper [13] that appeared on arXiv in October 2022. In a later version,
Goda replaced the term $\|k\|_{\infty}$ by $\min_{i\in{\rm supp}(k)}|k_{i}|$
and thus obtained tractability on an even larger class. I believe that this
replacement is only possible for the integration problem and it will lead to
intractability in the case of the approximation problem.

######  Remark 4.

Other unweighted classes of functions where the approximation problem is
tractable are high-dimensional functions that show some “low-dimensional
structure” like tensor products of univariate functions [1, 21, 24], sums of
low-dimensional functions [30], or compositions of univariate functions with
linear functionals [23]. Moreover, it was proven in [31] that the
approximation problem is weakly tractable for the class of analytic functions
defined on the cube with directional derivatives of all orders bounded by 1.
For the integration problem, we mention the star-discrepancy as a further
example [14].

######  Remark 5.

Recently, the papers [5, 19] also derived new bounds on the complexity of
sampling recovery in $L_{2}$ using non-linear algorithms. It is shown in [19,
Thm. 3.2] that the $n$-th minimal error is essentially bounded by the $s$-term
widths of the class $F_{d}$ in $L_{\infty}$ where it suffices to take $n$
quasi-linear in both $s$ and $d$. The paper [19] provides several examples of
classical smoothness classes $F_{d}$ where non-linear sampling algorithms have
a better rate of convergence than any linear algorithm. Here, we use instead a
bound by $s^{-1/2}$ times the $s$-term widths of the class $F_{d}$ in the
Wiener algebra, see Lemma 9 and Remark 12. The extra factor $s^{-1/2}$ eases
the tractability analysis significantly. It would be interesting to know
whether [19, Thm. 3.2] also leads to tractability results for the class
$F_{d}^{\log}$ or related ones.

######  Remark 6.

It is remarkable that the complexity bounds in this paper are not obtained
with a sophisticated choice of interpolation nodes but rather with i.i.d.
uniformly distributed sampling points, see Lemma 9. Recently, there has been
much interest in the surprising power of i.i.d. random information in
comparison to optimal information. I refer to the survey paper [28].

##  1\. A general complexity bound

Let $\mu$ be a probability measure on a set $D$ and let
$\mathcal{B}=\\{b_{k}\mid k\in I\\}$ be a countable orthonormal system in
$L_{2}(D,\mu)$ such that

(4) $C_{\mathcal{B}}\,:=\,\sup_{k\in I}\|b_{k}\|_{\infty}\,<\,\infty.$

We say that $\mathcal{B}$ is a bounded orthonormal system. Whenever
$(c_{k})_{k\in I}$ is an absolutely summable sequence of complex numbers, the
series $f=\sum_{k\in I}c_{k}b_{k}$ converges uniformly on $D$ and we have

$\hat{f}(k)\,:=\,\int_{D}f(x)\overline{b_{k}(x)}\,{\rm
d}\mu(x)\,=\,c_{k},\quad k\in I.$

We consider the Wiener algebra

(5) $\mathcal{A}\,:=\,\mathcal{A}(\mathcal{B})\,:=\,\left\\{f=\sum_{k\in
I}c_{k}b_{k}\,\,\Big{|}\,\,\|f\|_{\mathcal{A}}:=\sum_{k\in
I}|c_{k}|<\infty\right\\}$

and denote its unit ball by $\mathcal{A}_{0}(\mathcal{B})$. For a finite index
set $\Lambda\subset I$ and $f\in L_{2}$, we write

$P_{\Lambda}(f)\,:=\,\sum_{k\in\Lambda}\hat{f}(k)b_{k}$

and denote by $\mathcal{T}(\Lambda)$ the span of the functions $b_{k}$ with
$k\in\Lambda$. Given a class $F\subset\mathcal{A}$, we define the projection
error

$E_{\Lambda}^{\infty}(F)\,:=\,\sup_{f\in F}\|f-P_{\Lambda}f\|_{\infty}.$

We consider the following algorithm.

######  Algorithm 7.

Given a class $F\subset\mathcal{A}$, a finite index set $\Lambda\subset I$,
and a finite and non-empty point set $X\subset D$, we consider the algorithm
$A_{F,\Lambda,X}$ which maps a function $f\in F$ to any solution of

(6) $\displaystyle\min\|g\|_{\mathcal{A}}$ subject to $\displaystyle
g\in\mathcal{T}(\Lambda)$ and $\displaystyle\frac{1}{\\#X}\sum_{x\in
X}|f(x)-g(x)|^{2}\,\leq\,E_{\Lambda}^{\infty}(F)^{2}.$

The algorithm is well-defined since the function $g=P_{\Lambda}(f)$ satisfies
the constraints of the minimization problem. Moreover, the output of the
algorithm merely depends on $f$ via the data $(f(x))_{x\in X}$ and thus
$A_{F,\Lambda,X}$ is a sampling algorithm of the form (1).

######  Remark 8.

Algorithm 7 is a direct translation of a classical algorithm from compressed
sensing, called basis pursuit denoising. Namely, if we denote
$x=(\hat{g}(k))_{k\in\Lambda}$, $y=(f(x))_{x\in X}$ and
$G=(b_{k}(x))_{k\in\Lambda,x\in X}$, solving (6) is equivalent to the
$\ell^{1}$-minimization problem

(7) $\min\|x\|_{1}\quad\text{subject to}\quad\|Gx-y\|_{2}\leq
E_{\Lambda}^{\infty}(F)\sqrt{\\#X}.$

The insight that $\ell^{1}$-minimization can be of great advantage over linear
algorithms goes back at least to [2, 8, 11] and the problem of
$\ell^{2}$-recovery on $\ell^{1}$-balls in $\mathbb{R}^{m}$. Different methods
of computing a solution to (7) can be found in [10, Ch. 15].

We now give an error bound for Algorithm 7. The error is expressed in terms of
the best $s$-term widths of $F$ in $\mathcal{A}$, defined by

$\sigma_{s}(F,\mathcal{A})\,:=\,\sup_{f\in F}\,\inf_{\\#\Lambda\leq
s\vphantom{|^{|}}}\,\sum_{k\in I\setminus\Lambda}|\hat{f}(k)|.$

We derive our error bound as a simple implication of [26, Thm. 6.1], where we
reinterpret our data on $f$ as noisy data on $P_{\Lambda}f$. The key
ingredient to error bounds such as [26, Thm. 6.1] is the restricted isometry
property of the matrix $G=(b_{k}(x))_{k\in\Lambda,x\in X}$; For $m$ random
points, the matrix $m^{-1/2}G$ acts as a quasi-isometry on the set of vectors
with at most $s$ non-zero entries, where it suffices that $m$ is quasi-linear
in $s$ and logarithmic in $\\#\Lambda$, see [10, Thm. 12.32] or [26, Thm.
5.2]. This property implies that vectors $c$ from the unit ball in
$\ell^{1}(\Lambda)$ can be recovered well from data of the form $y=Gc+e$ if
the noise $e$ is sufficiently small, see [10, Thm. 6.12]. In our case, we want
to recover the vector $c=(\hat{f}(k))_{k\in\Lambda}$. The data is given by
$y=(f(x))_{x\in X}$ while $Gc=(P_{\Lambda}f(x))_{x\in X}$ and so the entries
of $e$ are bounded by $E_{\Lambda}^{\infty}(F)$.

######  Lemma 9 (Compare [26, Thm. 6.1]).

There is a universal constant $c\geq 1$ such that the following holds. Let
$\mathcal{B}=\\{b_{k}\mid k\in I\\}$ be a bounded orthonormal system with
respect to a probability measure $\mu$ and let
$F\subset\mathcal{A}(\mathcal{B})$. For $\gamma\in(0,1)$, $s\geq
2C_{\mathcal{B}}^{2}$, and $\Lambda\subset I$, let $X_{m}$ be a set of

$m\,\geq\,c\,\log(\gamma^{-1})\,s\log^{3}(s)\log(\\#\Lambda)$

independent random points with distribution $\mu$. Then, with probability at
least $1-\gamma$, we have for all $2\leq p\leq\infty$ and $f\in F$ that

$\|f-A_{F,\Lambda,X_{m}}(f)\|_{p}\,\leq\,c\,C_{\mathcal{B}}\,s^{-1/p}\sigma_{s}(F,\mathcal{A})+c\,s^{1/2-1/p}E_{\Lambda}^{\infty}(F).$

###### Proof.

We use [26, Thm. 6.1] with $b_{j}$ instead of $\phi_{j}$ and
$\omega_{j}=C_{\mathcal{B}}$. We choose $c_{0}$ as the maximum of all
constants appearing in [26, Thm. 6.1] and $c=c_{0}+1$. Then
$X_{m}=\\{t_{1},\ldots,t_{m}\\}$ satisfies the conclusion of [26, Thm. 6.1]
with probability $1-\gamma$. Given $g\in F$, we apply this conclusion to the
function $f=P_{\Lambda}g$ and samples $y_{\ell}=f(t_{\ell})+\xi_{\ell}$ with
$\xi_{\ell}=g(t_{\ell})-f(t_{\ell})$ and hence $\|\xi\|_{2}\leq
E_{\Lambda}^{\infty}(F)\sqrt{m}=:\eta$. Then the function $f^{\\#}$ defined in
[26, Thm. 6.1] equals the function $A_{F,\Lambda,X_{m}}(g)$. By the triangle
inequality and Jensen’s inequality, we have

$\|g-A_{F,\Lambda,X_{m}}(g)\|_{p}\leq
E_{\Lambda}^{\infty}(F)+\|f-f^{\\#}\|_{p}.$

As $\sigma_{s}(f)_{\omega,1}$ as defined in [26] is bounded by
$C_{\mathcal{B}}\sigma_{s}(F,\mathcal{A})$, we get

$\displaystyle\|f-f^{\\#}\|_{2}\,\leq\,c_{0}C_{\mathcal{B}}\,s^{-1/2}\sigma_{s}(F,\mathcal{A})+c_{0}E_{\Lambda}^{\infty}(F)\,=:\,R,$
$\displaystyle\|f-f^{\\#}\|_{\infty}\,\leq\,c_{0}C_{\mathcal{B}}\,\sigma_{s}(F,\mathcal{A})+c_{0}s^{1/2}E_{\Lambda}^{\infty}(F)\,=\,Rs^{1/2}.$

It only remains to apply
$\|h\|_{p}^{p}\leq\|h\|_{2}^{2}\cdot\|h\|_{\infty}^{p-2}$ for $h=f-f^{\\#}$. ∎

This leads to the following complexity bound. Recall the definition of the
(information) complexity of sampling recovery in $L_{2}$ from (2) which is
defined analogously for sampling recovery in $L_{p}$. We also define the
minimal size of an index set needed for a projection error $\varepsilon>0$,
i.e.,

(8) $N_{\mathcal{B}}^{\infty}(\varepsilon,F)\,:=\,\min\\{\\#\Lambda\mid
E_{\Lambda}^{\infty}(F)\leq\varepsilon\\}.$

######  Proposition 10.

There is a universal constant $C\geq 2$ such that, for any bounded orthonormal
system $\mathcal{B}$, any $F\subset\mathcal{A}_{0}(\mathcal{B})$, all $1\leq
p<\infty$ and $\varepsilon\in(0,1)$, we have

$n(\varepsilon,F,L_{p})\,\leq\,C\,\tilde{\varepsilon}^{-r}\log^{3}(\tilde{\varepsilon}^{-r})\,\log
N_{\mathcal{B}}^{\infty}\big{(}\tilde{\varepsilon}^{r/2},F\big{)}$

with $\tilde{\varepsilon}=\varepsilon/(CC_{\mathcal{B}})$ and
$r=\max\\{p,2\\}$.

###### Proof.

The case $p<2$ follows from the case $p=2$, so let $p\geq 2$. Put $C=2c$ with
$c$ from Lemma 9. We fix $\gamma={\rm e}^{-1}$ and choose $\Lambda\subset I$
with $\\#\Lambda\leq N_{\mathcal{B}}^{\infty}(\tilde{\varepsilon}^{p/2},F)$
and $E_{\Lambda}^{\infty}(F)\leq\tilde{\varepsilon}^{p/2}$. Moreover, we put
$s=\lceil\tilde{\varepsilon}^{-p}\rceil$ so that $s\geq 2C_{\mathcal{B}}^{2}$.
Observing that

$\sigma_{s}(F,\mathcal{A})\,\leq\,1,$

Lemma 9 yields the statement. ∎

Proposition 10 leads to our main result on the tractability of
$L_{p}$-approximation. Here, given classes $F_{d}\subset\mathcal{A}_{d}$ for
each $d\in\mathbb{N}$, we say that sampling recovery on $F_{d}$ in $L_{p}$ is
polynomially tractable if

$\exists C,q,r\geq 0\colon\forall d\in\mathbb{N}\
\forall\varepsilon>0\colon\quad
n(\varepsilon,F_{d},L_{p})\,\leq\,Cd^{q}\varepsilon^{-r}.$

Recall the definition of $C_{\mathcal{B}}$, $\mathcal{A}_{0}(\mathcal{B})$ and
$N_{\mathcal{B}}^{\infty}(\varepsilon,F)$ from (4), (5) and (8).

######  Theorem 11.

For every $d\in\mathbb{N}$, let $\mathcal{B}_{d}$ be a bounded orthonormal
system and $F_{d}\subset\mathcal{A}_{0}(\mathcal{B}_{d})$. Assume that there
are positive constants $c_{1},c_{2},\alpha,\beta,$ and $\gamma$ such that

$C_{\mathcal{B}_{d}}\,\leq\,c_{1}d^{\alpha}\qquad\text{and}\qquad
N_{\mathcal{B}}^{\infty}(\varepsilon,F_{d})\,\leq\,\exp(c_{2}d^{\beta}\varepsilon^{-\gamma})$

for all $\varepsilon>0$ and $d\in\mathbb{N}$. Then the problem of sampling
recovery on $F_{d}$ in $L_{p}$ is polynomially tractable for all $1\leq
p<\infty$.

That is, sampling recovery in $L_{p}$ on classes with absolutely convergent
(generalized) Fourier series is tractable even if the number of Fourier
coefficients needed for an $\varepsilon$-approximation of functions $f\in
F_{d}$ in the uniform norm grows super-exponentially in $d$ and
$\varepsilon^{-1}$, as long as the growth is not double-exponential.

######  Remark 12 (Rate of convergence).

For our tractability analysis, it was enough to observe that
$\sigma_{n}(F,\mathcal{A})\leq 1$ for all $F\subset\mathcal{A}_{0}$. In fact,
it would not help for the classes considered in this paper to take the decay
of the widths $\sigma_{n}(F,\mathcal{A})$ into account as $n$ would have to be
exponentially large in $d$ if we want widths significantly smaller than one.
This is different if one is interested in studying the rate of convergence of
the $n$th minimal error

${\rm err}(n,F_{d},L_{p})\,:=\,\inf_{\begin{subarray}{c}x_{1},\dots,x_{n}\in
D\\\ \phi\colon\mathbb{C}^{n}\to L_{p}\end{subarray}}\,\sup_{f\in
F}\,\|f-A_{n}(f)\|_{p}.$

Often (e.g., for Sobolev classes of mixed smoothness) one can choose the size
of $\Lambda$ polynomial in $n$ in order to obtain from Lemma 9 that

${\rm
err}(Cn\log^{4}n,F,L_{p})\,\lesssim\,n^{-1/p}\,\sigma_{n}(F,\mathcal{A}).$

In comparison, the paper [19] recently revealed estimates of the form

${\rm
err}(Cn\log^{4}n,F,L_{2})\,\lesssim\,\sigma_{n}(F,\mathcal{B})_{L_{\infty}}$

by means of the best $n$-term widths in $L_{\infty}$ instead of $\mathcal{A}$,
where the power of the logarithm can be reduced to 3 for the trigonometric
system.

##  2\. Results for specific classes

We now consider the Fourier system

$\mathcal{B}_{d}\,:=\,\\{b_{k}=e^{2\pi i\langle k,\cdot\rangle}\mid
k\in\mathbb{Z}^{d}\\}$

that is orthonormal in $L_{2}([0,1]^{d})$ with the Lebesgue measure and
satisfies $C_{\mathcal{B}_{d}}=1$. We obtain the classical Wiener algebra
$\mathcal{A}_{d}:=\mathcal{A}(\mathcal{B}_{d})$ of periodic functions with
absolutely convergent Fourier series. We consider three examples, starting
with the class $F_{d}^{\log}$ from the introduction. The following corollary
also proves Theorem 1.

######  Corollary 13.

Let $2\leq p<\infty$. There is a constant $c>0$ such that for all $d\geq 2$
and $\varepsilon\in(0,1/2)$, we have

$n(\varepsilon,F_{d}^{\log},L_{p})\,\leq\,c\,d\,\varepsilon^{-3p/2}\log^{3}(\varepsilon^{-p}).$

###### Proof.

For any $m\in\mathbb{N}$ and $f\in F_{d}^{\log}$, we have

$\|f-P_{[-m,m]^{d}}f\|_{\infty}\,\leq\,\log^{-1}(m+1)$

and thus

$N_{\mathcal{B}_{d}}^{\infty}(\varepsilon,F_{d}^{\log})\,\leq\,\exp(2d\varepsilon^{-1}).$

Now the bound is obtained from Proposition 10. ∎

As a second example, we consider Sobolev spaces of mixed smoothness $s>1/2$
(also called unweighted Korobov spaces), namely

$H_{\rm mix}^{s,d}\,:=\,\bigg{\\{}f\in C([0,1]^{d})\ \Big{|}\
\sum_{k\in\mathbb{Z}^{d}}|\hat{f}(k)|^{2}\prod_{i=1}^{d}\max\\{1,|k_{i}|^{2s}\\}\leq
1\bigg{\\}}.$

It is well known that $L_{2}$-approximation on $H_{\rm mix}^{s,d}$ suffers
from the curse of dimensionality, see, e.g., [15]. The curse is relinquished
in the presence of absolutely summable Fourier coefficients, namely, for the
class

$F^{s,d}_{\rm mix}\,:=\,\bigg{\\{}f\in H_{\rm mix}^{s,d}\ \Big{|}\
\sum_{k\in\mathbb{Z}^{d}}|\hat{f}(k)|\leq 1\bigg{\\}}.$

######  Remark 14 (Weighted Korobov spaces).

A classical approach to relinquish the curse of dimensionality on $H_{\rm
mix}^{s,d}$ is by introducing (product) weights
$1\geq\gamma_{1}\geq\gamma_{2}\geq\ldots>0$ and considering the class

$H_{\rm mix}^{s,d,\gamma}\,:=\,\bigg{\\{}f\in
C([0,1]^{d})\mid\sum_{k\in\mathbb{Z}^{d}}|\hat{f}(k)|^{2}\prod_{i=1}^{d}\max\\{1,\gamma_{i}^{-1}|k_{i}|^{2s}\\}\leq
1\bigg{\\}}.$

see, e.g., [9, 22, 25] and [7, Ch. 13] and the references therein. The
sequence $\gamma$ models a decreasing importance of the variables. It is known
that $L_{2}$-approximation on $H_{\rm mix}^{s,d,\gamma}$ is polynomially
tractable if the weights decay fast enough, namely, if there is a constant
$K>0$ with

$\sum_{i=1}^{d}\gamma_{i}\,\leq\,K\log(d),\quad d\geq 2,$

see [25, Thm. 1]. The present approach is more general: For any $\gamma$ as
above, a simple computation shows the inclusion

(9) $H_{\rm mix}^{s,d,\gamma}\subset C\,d^{q}\,F^{s,d}_{\rm mix}$

with constants $C,q>0$ that only depend on $K$ and $s$.

######  Corollary 15.

For any $s>1/2$ and $1\leq p<\infty$, the problem of $L_{p}$-approximation is
polynomially tractable on $F^{s,d}_{\rm mix}$.
More precisely, for $p\geq 2$, there is a constant $C>0$, depending only on
$s$ and $p$, such that for all $d\geq 2$ and $\varepsilon\in(0,1/2)$, we have

$n(\varepsilon,F^{s,d}_{\rm
mix},L_{p})\,\leq\,C\,d^{2}\log(d)\,\varepsilon^{-p}\log^{4}(\varepsilon^{-1}).$

###### Proof.

For $m\in\mathbb{N}$ and $f\in F^{s,d}_{\rm mix}$, we get from Hölder’s
inequality that

$\|f-P_{[-m,m]^{d}}f\|_{\infty}^{2}\,\leq\,\bigg{(}\sum_{\|k\|_{\infty}>m}|\hat{f}(k)|\bigg{)}^{2}\,\leq\,\sum_{\|k\|_{\infty}>m}\prod_{i=1}^{d}\max\\{1,|k_{i}|\\}^{-2s}\\\
\leq\,d\cdot\sum_{\begin{subarray}{c}k\in\mathbb{Z}^{d}\\\
|k_{1}|>m\end{subarray}}\prod_{i=1}^{d}\max\\{1,|k_{i}|\\}^{-2s}\,\leq\,d\cdot
c_{s}^{d-1}\sum_{\begin{subarray}{c}k_{1}\in\mathbb{Z}\\\
|k_{1}|>m\end{subarray}}|k_{1}|^{-2s}\,\leq\,d\cdot c_{s}^{d}\cdot m^{-2s+1}$

with $c_{s}>0$ depending only on $s$, and thus

$N_{\mathcal{B}_{d}}^{\infty}(\varepsilon,F^{s,d}_{\rm
mix})\,\leq\,\big{(}d\,c_{s}^{d}\,\varepsilon^{-2}\big{)}^{d/(2s-1)}.$

Now the bound is obtained from Proposition 10. ∎

######  Remark 16.

In regard of (9), we note that

(10) $n(\varepsilon,rF_{d},L_{p})\,\leq\,n(\varepsilon/(2r),F_{d},L_{p})$

holds for any convex and symmetric class $F_{d}$ and any $r\geq 1$, so that
tractability on $F^{s,d}_{\rm mix}$ indeed implies tractability on the classes
$H_{\rm mix}^{s,d,\gamma}$. Equation (10) follows from the optimality of
homogeneous algorithms for linear problems, see [20, Thm. 1].

As a third example, we consider the class of Hölder functions that has been
studied by Dick [6] for numerical integration, namely,

$F_{d}^{\alpha}\,:=\,\bigg{\\{}f\in C([0,1]^{d})\ \Big{|}\
\sum_{k\in\mathbb{Z}^{d}}|\hat{f}(k)|\leq 1\,\land\,\sup_{x\neq
y}\frac{|f(x)-f(y)|}{\|x-y\|_{2}^{\alpha}}\leq 1\bigg{\\}}$

with some $0<\alpha\leq 1$. Here, the supremum is taken over all
$x,y\in\mathbb{R}^{d}$, considering $f$ as a 1-periodic function on
$\mathbb{R}^{d}$.

######  Corollary 17.

For any $\alpha\in(0,1]$ and $1\leq p<\infty$, the problem of
$L_{p}$-approximation is polynomially tractable on $F^{s,d}_{\rm mix}$.
More precisely, for $p\geq 2$, there is a constant $c>0$, depending only on
$\alpha$, such that for all $d\geq 2$ and $\varepsilon\in(0,1/2)$, we have

$n(\varepsilon,F_{d}^{\alpha},L_{p})\,\leq\,c\,d^{2}\log^{2}(d)\,\varepsilon^{-p}\log^{4}(\varepsilon^{-1}).$

###### Proof.

Let $S_{m}=P_{[-m,m]}$ be the univariate Fourier sum operator and
$C^{\alpha}_{\rm per}$ be the class of univariate $\alpha$-Hölder continuous
functions with Hölder constant one. It is well known that

$\|S_{m}f\|_{\infty}\,\leq\,c_{1}\log(m)\|f\|_{\infty}$

for all $m\geq 2$ and all $f\in L_{\infty}([0,1])$ and

$\|f-S_{m}f\|_{\infty}\,\leq\,c_{2}\log(m)m^{-\alpha}$

for all $m\geq 2$ and all $f\in C^{\alpha}_{\rm per}$, where $c_{1},c_{2}>0$
depend on nothing but $\alpha$, see [18]. Let now $f\in F_{d}^{\alpha}$. For
fixed $x_{1},\dots,x_{i-1},x_{i+1},\dots,x_{d}\in[0,1]$, the function $f$ as a
function of $x_{i}$ is in $C^{\alpha}_{\rm per}$. By the triangle inequality,

$\displaystyle\|f-P_{[-m,m]^{d}}f\|_{\infty}\,\leq\,\sum_{i=1}^{d}\|S_{m}^{x_{1}}\cdots
S_{m}^{x_{i-1}}(f-S_{m}^{x_{i}}f)\|_{\infty}$
$\displaystyle\leq\,\sum_{i=1}^{d}(c_{1}\log(m))^{i-1}\cdot
c_{2}\log(m)m^{-\alpha}\,\leq\,\frac{(c_{3}\log(m))^{d}}{m^{\alpha}}.$

This term is smaller than $\varepsilon$ for
$m\geq\exp(c_{4}d\log^{2}(d)\log(\varepsilon^{-1}))$, thus

$\log
N_{\mathcal{B}_{d}}^{\infty}(\varepsilon,F_{d}^{\alpha})\,\leq\,c_{5}d^{2}\log^{2}(d)\log(\varepsilon^{-1}).$

Now the bound is obtained from Proposition 10. ∎

It is interesting to note that the polynomial order $d^{2}\varepsilon^{-2}$
that we obtained for $L_{2}$-approximation is the same as the one obtained in
[6] for numerical integration. It is surely an interesting open problem to
improve the complexity bounds in Corollaries 13, 15 and 17 and/or provide
matching lower bounds.

##  3\. Intractability with linear algorithms

We only have tractability for the class $F_{d}^{\log}$ thanks to non-linear
algorithms. If we restrict to linear algorithms, we still have the curse of
dimensionality. This is implied by a result of [12], see also [29, Lem. 3.2].
Namely, it is known that for any $m,n\in\mathbb{N}$, $m>n$, and any linear
mapping $T_{n}\colon\mathbb{C}^{m}\to\mathbb{C}^{m}$ of rank $n$, there is
some $x\in\mathbb{C}^{m}$ with $\|x\|_{1}\leq 1$ and

(11) $\|x-T_{n}x\|_{2}\,\geq\,\left(\frac{m-n}{m}\right)^{1/2}.$

######  Lemma 18.

The problem of $L_{2}$-approximation on the classes $F_{d}^{\log}$ with linear
algorithms suffers from the curse of dimensionality. For any linear mapping
$A_{n}\colon F_{d}^{\log}\to L_{2}$ with rank $n\leq 5^{d}/2$ we have

${\rm err}(A_{n},F_{d}^{\log},L_{2})\,\geq\,\frac{1}{\sqrt{2}}.$

Note that the curse of dimensionality even holds for the class of all linear
algorithms, using arbitrary linear information, instead of only sampling based
algorithms.

###### Proof.

Let $A_{n}\colon F_{d}^{\log}\to L_{2}$ be a linear operator of rank $n\leq
5^{d}/2$ and let $\Lambda=\\{-2,-1,0,1,2\\}^{d}$ so that $\\#\Lambda=5^{d}$.
Consider the linear mapping
$T\colon\mathbb{C}^{\Lambda}\to\mathcal{T}(\Lambda)$ that maps a vector of
coefficients to the corresponding trigonometric polynomial. Then
$T_{n}=T^{-1}P_{\Lambda}A_{n}T\colon\mathbb{C}^{\Lambda}\to\mathbb{C}^{\Lambda}$
is a linear mapping with rank at most $n$ and (11) yields that there is some
$x\in\mathbb{C}^{\Lambda}$ with $\|x\|_{1}\leq 1$ such that

$\|Tx-
P_{\Lambda}A_{n}Tx\|_{2}\,=\,\|x-T_{n}x\|_{2}\,\geq\,\left(\frac{5^{d}-5^{d}/2}{5^{d}}\right)^{1/2}\,=\,\frac{1}{\sqrt{2}}.$

The function $f=Tx$ is contained in $F_{d}^{\log}$. Moreover,
$f\in\mathcal{T}(\Lambda)$ and thus

(12)
$\|f-A_{n}f\|_{2}\,\geq\,\|f-P_{\Lambda}A_{n}f\|_{2}\,\geq\,\frac{1}{\sqrt{2}},$

as claimed. ∎

Clearly, the lower bound (12) also holds for the uniform instead of the
$L_{2}$-norm. Since linear algorithms are optimal for the recovery problem in
the uniform norm, see [4], Lemma 18 implies that the problem of uniform
recovery on the class $F_{d}^{\log}$ is intractable. Namely, the following
holds for all algorithms of the form

(13) $A_{n}\colon F_{d}^{\log}\to L_{\infty}([0,1]^{d}),\quad
A_{n}(f)=\phi(L_{1}(f),\ldots,L_{n}(f))$

where $L_{1},\ldots,L_{n}\colon F_{d}^{\log}\to\mathbb{C}$ are linear
functionals, possibly chosen adaptively, and $\phi\colon\mathbb{C}^{n}\to
L_{\infty}$. This includes the sampling algorithms of the form (1).

######  Corollary 19.

The problem of (sampling) recovery in the uniform norm on the classes
$F_{d}^{\log}$ suffers from the curse of dimensionality. For any mapping
$A_{n}\colon F_{d}^{\log}\to L_{\infty}$ of the form (13) with $n\leq 5^{d}/2$
we have

$\sup_{f\in F_{d}^{\log}}\|f-A_{n}(f)\|_{\infty}\,\geq\,\frac{1}{\sqrt{2}}.$

In a similar way, one can show that polynomial tractability cannot be achieved
with linear algorithms for the classes $F^{s,d}_{\rm mix}$ and
$F_{d}^{\alpha}$ from Section 2. We omit the details.

Acknowledgement. This research was funded by the Austrian Science Fund (FWF)
Project M 3212-N. For the purpose of open access, I have applied a CC BY
public copyright licence to any author accepted manuscript version arising
from this submission. I thank Erich Novak, Mario Ullrich, Tino Ullrich and a
referee for helpful comments.

## References

  * [1] M. Bachmayr, W. Dahmen, R. DeVore, and L. Grasedyck, Approximation of high-dimensional rank one tensors, Constr. Approx., 39:385–395, 2014.
  * [2] E. J. Candes and T. Tao, Near-optimal signal recovery from random projections: universal encoding strategies?, IEEE Trans. Inform. Theory, 52:5406–5425, 2006.
  * [3] L. Chen, H. Jiang, On the information complexity for integration in subspaces of the Wiener algebra, J. Complexity, 81:101819, 2024.
  * [4] J. Creutzig and P. Wojtaszczyk, Linear vs. nonlinear algorithms for linear problems, J. Complexity, 20:807–820, 2004.
  * [5] F. Dai, V. Temlyakov, Random points are good for universal discretization, J. Math. Anal. Appl., 529(1):127570, 2024.
  * [6] J. Dick, Numerical integration of Hölder continuous, absolutely convergent Fourier, Fourier cosine, and Walsh series, J. Approx. Theory, 183:14–30, 2014.
  * [7] J. Dick, P. Kritzer, F. Pillichshammer. Lattice Rules. Springer, Cham, 2022.
  * [8] D. L. Donoho, Compressed sensing, IEEE Trans. Inform. Theory, 52:1289–1306, 2006.
  * [9] A. Ebert, F. Pillichshammer, Tractability of approximation in the weighted Korobov space in the worst-case setting – a complete picture, J. Complexity, 67:101571, 2021.
  * [10] S. Foucart, H. Rauhut, A Mathematical Introduction to Compressive Sensing, Birkhäuser, New York, 2013.
  * [11] A. Yu. Garnaev and E. D. Gluskin, The widths of a Euclidean ball, Soviet Math. Dokl., 30:200–204, 1984.
  * [12] E. D. Gluskin, Norms of random matrices and widths of finite-dimensional sets, Math. USSR Sb., 48:173–182, 1984.
  * [13] T. Goda, Polynomial tractability for integration in an unweighted function space with absolutely convergent Fourier series, Proc. Amer. Math. Soc., 151:3925–3933, 2023.
  * [14] S. Heinrich, E. Novak, G. W. Wasilkowski, and H. Woźniakowski, The inverse of the star-discrepancy depends linearly on the dimension, Acta Arith., 96: 279-302, 2001.
  * [15] F. Hickernell and H. Woźniakowski, Tractability of Multivariate Integration for Periodic Functions. J. Complexity, 17(4):660–682, 2001.
  * [16] A. Hinrichs, E. Novak, M. Ullrich, H. Woźniakowski, The curse of dimensionality for numerical integration of smooth functions. Math. Comp., 83:2853–2863, 2014.
  * [17] A. Hinrichs, E. Novak, M. Ullrich, H. Woźniakowski. Product rules are optimal for numerical integration in classical smoothness spaces. J. Complexity, 38:39–49, 2017.
  * [18] D. Jackson, The theory of Approximation, AMS Colloquium Publication Volume XI, New York, 1930.
  * [19] T. Jahn, T. Ullrich, F. Voigtlaender. Sampling numbers of smoothness classes via $\ell^{1}$-minimization, J. Complexity, 79:101786, 2023.
  * [20] D. Krieg, P. Kritzer, Homogeneous algorithms and solvable problems on cones, preprint, arXiv:2311.15767.
  * [21] D. Krieg, D. Rudolf, Recovery algorithms for high-dimensional rank one tensors, J. Approx. Theory, 237:17–29, 2019.
  * [22] F. Y. Kuo, Ch. Schwab and I. H. Sloan, Quasi-Monte Carlo methods for high-dimensional integration: The standard (weighted Hilbert space) setting and beyond, ANZIAM J., 53(1):1–37, 2011.
  * [23] S. Mayer, T. Ullrich, J. Vybíral. Entropy and sampling numbers of classes of ridge functions. Constr. Approx., 42(2):231–264, 2015.
  * [24] E. Novak, D. Rudolf, Tractability of the approximation of high-dimensional rank one tensors, Constr. Approx., 43:1–13, 2016.
  * [25] E. Novak, I. H. Sloan, H. Woźniakowski, Tractability of approximation for weighted Korobov spaces on classical and quantum computers, Found. Comput. Math. 4(2):121–156, 2004.
  * [26] H. Rauhut and R. Ward, Interpolation via weighted $\ell^{1}$-minimization. Appl. Comput. Harmon. Anal. 40, 321–351, 2016.
  * [27] I. H. Sloan and H. Woźniakowski, When are quasi-Monte Carlo algorithms efficient for high dimensional integrals? J. Complexity 14:1–33, 1998.
  * [28] M. Sonnleitner, M. Ullrich, On the power of iid information for linear approximation, J. Appl. Num. Anal., 1:88–126, 2023.
  * [29] J. Vybíral, Widths of embeddings in function spaces. J. Complexity, 24(4):545–570, 2008.
  * [30] G.W. Wasilkowski, H. Woźniakowski, Finite-order weights imply tractability of linear multivariate problems, J. Approx. Theory, 130(1):57–77, 2004.
  * [31] G. Xu, On weak tractability of the Smolyak algorithm for approximation problems, J. Approx. Theory, 192:347–361, 2015.