# Tensor Quantum Programming

A. Termanova  These authors contributed equally: A. Termanova and Ar.
Melnikov.  Terra Quantum AG, Kornhausstrasse 25, 9000 St. Gallen, Switzerland
Ar. Melnikov  These authors contributed equally: A. Termanova and Ar.
Melnikov.  Terra Quantum AG, Kornhausstrasse 25, 9000 St. Gallen, Switzerland
E. Mamenchikov  Terra Quantum AG, Kornhausstrasse 25, 9000 St. Gallen,
Switzerland  N. Belokonev  Terra Quantum AG, Kornhausstrasse 25, 9000 St.
Gallen, Switzerland  S. Dolgov  Terra Quantum AG, Kornhausstrasse 25, 9000 St.
Gallen, Switzerland  A. Berezutskii  Terra Quantum AG, Kornhausstrasse 25,
9000 St. Gallen, Switzerland  R. Ellerbrock  Terra Quantum AG, Kornhausstrasse
25, 9000 St. Gallen, Switzerland  C. Mansell  Terra Quantum AG,
Kornhausstrasse 25, 9000 St. Gallen, Switzerland  M. Perelshtein  Terra
Quantum AG, Kornhausstrasse 25, 9000 St. Gallen, Switzerland

###### Abstract

Running quantum algorithms often involves implementing complex quantum
circuits with such a large number of multi-qubit gates that the challenge of
tackling practical applications appears daunting. To date, no experiments have
successfully demonstrated a quantum advantage due to the ease with which the
results can be adequately replicated on classical computers through the use of
tensor network algorithms. Additionally, it remains unclear even in theory
where exactly these advantages are rooted within quantum systems because the
logarithmic complexity commonly associated with quantum algorithms is also
present in algorithms based on tensor networks. In this article, we propose a
novel approach called Tensor Quantum Programming, which leverages tensor
networks for hybrid quantum computing. Our key insight is that the primary
challenge of algorithms based on tensor networks lies in their high ranks
(bond dimensions). Quantum computing offers a potential solution to this
challenge, as an ideal quantum computer can represent tensors with arbitrarily
high ranks in contrast to classical counterparts, which indicates the way
towards quantum advantage. While tensor-based vector-encoding and state-
readout are known procedures, the matrix-encoding required for performing
matrix-vector multiplications directly on quantum devices remained unsolved.
Here, we developed an algorithm that encodes Matrix Product Operators into
quantum circuits with a depth that depends linearly on the number of qubits.
It demonstrates effectiveness on up to 50 qubits for several matrices
frequently encountered in differential equations, optimization problems, and
quantum chemistry. We view this work as an initial stride towards the creation
of genuinely practical quantum algorithms.

##  I Introduction

Quantum computing has garnered significant attention in both academia and
industry. Of particular interest is hybrid quantum computing, where the main
idea is to efficiently combine classical and quantum computation [1, 2].
Considerable research efforts are being devoted to areas ranging from
variational algorithms [3] to quantum machine learning [4, 5]. On the one
hand, a clear understanding of how these methods provide advantages over
purely classical ones is still somewhat lacking [6, 7]. On the other hand,
important insights are being made [8]. Despite the ongoing theoretical
progress, it is still very difficult to devise algorithms that would allow
even an ideal quantum computer (QC) to efficiently solve practical tasks, such
as simulating physical or chemical systems [9, 10]. Hence, the purpose of this
work is to develop a methodology for creating effective quantum algorithms.

In the circuit model of quantum computation, quantum algorithms consist of
three main steps: data loading/encoding (quantum state preparation), data
processing (multi-qubit quantum gates), and measurement (quantum state
tomography). The primary challenge facing gate-based quantum computing is the
significant complexity involved in implementing these steps for arbitrary
states and multi-qubit gates. We elaborate on these complexities later in the
paper but in short, the complexity of decomposing arbitrary multi-qubit gates
into single- and two-qubit gates is exponential with respect to the number of
qubits [11, 12, 13, 14], and similar results hold for the preparation and
readout of arbitrary states [15, 16].

Consequently, in our present work, we suggest restricting the considered
states and matrices to those in the tensor train (TT) format [17], a simple
type of tensor network (TN) [18]. The main reason for this is that such states
and matrices possess small entanglement and therefore can be efficiently
realized via shallow circuits. For example, there are protocols for
implementing the encoding of a Matrix Product State (MPS) into shallow
circuits [19, 20, 21, 22, 23], which also work well with poor connectivity
between qubits. MPS and TT-vector are identical from the mathematical point of
view – we use them interchangeably hereinafter. Moreover, any quantum circuit
can be naturally regarded as a TN [24]. This insight leads us to shift our
focus towards the naturally arising inverse task: encoding a TN into a quantum
circuit. By focusing on encoding TTs, we leverage the structured nature of
TNs, bypassing the complexity associated with arbitrary encoding.

The rest of the paper is organized as follows. In Sec. II, we present a
concise exposition of the mathematical principles underlying TTs. In Sec. III,
we introduce the hybrid computing concept of Tensor Quantum Programming. In
Sec. IV, we briefly review the existing algorithms for encoding MPSs into a
quantum computer. In Sec. V.1, we propose a method for encoding TT-matrices
into quantum circuits, which underlies our scheme of Tensor Quantum
Programming. In Sec. V.2, we present numerical results investigating the
efficacy of our encoding algorithm. In Sec. VI, we explain the strategy for
achieving efficient data-readout through MPS tomography. In Sec. VII, we
provide a detailed overview of the application domains wherein the proposed
scheme of hybrid computations can offer advantages through the use of a QC.

##  II Preliminaries

In this section, we give a brief introduction to the mathematics behind TTs.
Any $n$-dimensional tensor can be represented or approximated in TT format
[17], which is a contracted chain of $n$ lower-order tensors shown in Fig. 1.

Figure 1: Graphical notation showing an $n$-dimensional tensor in its TT
format. Circles and blocks indicate tensors. The dangling vertical bonds
represent physical indices. Unless stated otherwise, we assume that the
dimension of physical indices equals 2. The connected lines denote tensor
contractions over virtual indices.

Such a decomposition is also known as the MPS format [18]. The term “state”
refers to the original use of MPS to describe quantum states. One can write a
quantum state $\ket{\psi}$ of $n$ qubits in the MPS format by expanding the
state vector in the computational basis into a TT:

$\ket{\psi}_{\text{\tiny{MPS}}}=\sum\limits_{\\{s\\}}B\text{\small{[1]}}^{s_{1}}_{1\alpha_{1}}B\text{\small{[2]}}^{s_{2}}_{\alpha_{1}\alpha_{2}}\dots
B\text{\small{[n]}}^{s_{n}}_{\alpha_{n-1}1}\ket{s_{1}s_{2}\dots s_{n}},$ (1)

where the terms $B\text{\small{[i]}}^{s_{i}}_{\alpha_{i-1},\alpha_{i}}$ are
$n$ different $3$-dimensional tensors (TT-cores), and we use the Einstein
summation convention that repeated indices are summed over. Each tensor
contains a physical index $s_{i}\in\\{1,2\\}$ and bond indices
$\alpha_{i}\in\\{1,2,\dots,r_{i}\\}$, as depicted in Fig. 2. We call the
maximum value of the bond index, $r$, the bond dimension, or the MPS (TT-)
rank.

Figure 2:  MPS representation of the state vector $\ket{\psi}$ defined in Eq.
(1).

A crucial connection exists between the entanglement of a quantum state and
its tensor train rank. Specifically, for an MPS $\ket{\psi}$ with a rank $r$,
the entropy of the reduced subsystem
$\rho=\text{Tr}_{2}{[\ket{\psi}\bra{\psi}]}$ follows an upper bound given by
$S(\rho)\leq\log{r}$ [18]. Therefore, MPS states are often called weakly
entangled.

Similar to the MPS format for vectors, the MPO (Matrix Product Operator) is a
TT representation for matrices or linear operators. The MPO representation of
some linear operator, $H$, on $n$-qubits reads as follows:

$H_{\text{\tiny{MPO}}}=\sum\limits_{\\{s\\},\\{l\\}}A\text{\small{[1]}}^{s_{1}l_{1}}_{1\alpha_{1}}A\text{\small{[2]}}^{s_{2}l_{2}}_{\alpha_{1}\alpha_{2}}\dots
A\text{\small{[n]}}^{s_{n}l_{n}}_{\alpha_{n-1}1}\times\\\ \ket{s_{1}s_{2}\dots
s_{n}}\bra{l_{1}l_{2}\dots l_{n}},$ (2)

where now the TT-cores
$A\text{\small{[i]}}^{s_{i}l_{i}}_{\alpha_{i-1}\alpha_{i}}$ are $4$-index
tensors with two physical indeces $s_{i}$ and $l_{i}\in\\{1,2\\}$, as depicted
in Fig. 3. The rank $r$ of an MPO is again the maximum value among all bond
indexes $\alpha_{i}$.

Figure 3:  MPO representation of the linear operator $H$ defined in Eq. (2).

The TT format supports algebraic operations on tensors, allowing norms, dot
products, element-wise sums, and products of tensors to be computed without
involving the whole tensor but working only with TT-cores [17]. This provides
a computationally efficient tool for working with vectors, matrices, quantum
states, and even data in machine learning [25]. For instance, discretized
smooth functions are usually well-represented in the MPS format [26].
Moreover, there is an exact MPS decomposition for trigonometric, exponential,
polynomial, and rational functions [27]. When discussing matrices, the MPO
format frequently offers accurate approximations for structured matrices.
Take, for instance, Toeplitz matrices, which consistently exhibit an MPO rank
of $3$ [28]. Transforming any MPS into a diagonal MPO matrix preserves its
rank [17].

##  III Tensor Quantum Programming Paradigm

Figure 4: Our hybrid paradigm of Tensor Quantum Programming
$TN\xrightarrow{}QC\xrightarrow{}TN$. Quantum state preparation is performed
using an MPS-encoding protocol [29, 23]. Matrix-vector multiplication (the
main operation on a QC) is implemented with the MPO-encoding scheme presented
in this work. Quantum state tomography is achieved using an MPS tomography
algorithm [30]. In the presented paradigm, the states at the input and output
of a quantum computer are effectively represented in the form of TNs, and all
of the steps have a logarithmic complexity in terms of the size of the
vectors.

We start with the key principles of tensor networks and their connection to
quantum computing. A key characteristic of MPS (TT-vector) and MPO (TT-matrix)
is the rank, also known as the bond dimension, which determines the degree of
entanglement or correlations in the vector or matrix. If the rank is small,
then such a representation is very efficient and allows the number of
parameters required for data-encoding to be compressed from exponential
scaling to polynomial scaling. In addition, as mentioned in the previous
section, it is possible to prepare (step 1) and read out (step 3) a MPS on a
quantum computer efficiently with linear complexity in the number of qubits
[23, 31]. Remarkably, as we demonstrate below, it is possible to efficiently
implement one of the main operations – the Matvec operation (matrix-vector
multiplication, step 2) – on a quantum computer.

Consequently, this allows us to introduce the hybrid computing concept of
Tensor Quantum Programming: $TN\xrightarrow{}QC\xrightarrow{}TN$, shown in
Fig. 4. Firstly, we assume that an algorithm for solving a particular
numerical problem is written solely in the language of tensor networks. That
is, only operations with tensor networks must be present. It is worth noting
that often the algorithm consists of multiple repeated iterations.
Subsequently, the classical implementation of this algorithm is initiated.
This approach is chosen because, in many instances, the initial ranks tend to
be small. Commencing the implementation classically proved to be considerably
more efficient in these cases compared to an immediate quantum implementation.
Moreover, should we observe an increase in the ranks of the solution, we
proceed to transfer the vector from the intermediate solution, obtained during
an intermediate iteration, to a quantum computer, employing the MPS
preparation algorithm (step 1 in Fig. 4) [23]. The following iterations of the
algorithm in the presence of large ranks of the solution are carried out on a
quantum computer via the method developed in Sec. V. As the ranks diminish
once more, we can transfer the intermediate solution back to a classical
computer and resume the process there, employing the MPS tomography algorithm
[30].

This scheme enables us to identify the advantage of a quantum computer over a
classical one. It is that much larger values of TT-ranks (the entanglement of
the state) are accessible. That is, problems with stronger correlations can be
solved since the complexity of quantum calculations is determined only by the
complexity of the quantum circuit and does not depend on the state vector’s
entanglement, which can grow unlimitedly during execution.

Further, in Sec. V, we show that the complexity of the circuit implementing
the multiplication of the matrix encoded into it is determined by the matrix
rank, which in many practical problems is often small [32, 33, 34, 35].
Classically, tensor computing is limited by the ranks of all intermediate
computation results, while quantum computing requires low ranks to prepare a
quantum state efficiently, but while computing, the ranks can increase without
bounds. Yet, for effective quantum state reading, we must reduce the ranks.
Formulated exclusively in terms of tensor networks, we opt for the most
suitable hardware for each phase, depending on the size of the vector’s rank:
classical when the rank is small and quantum when it is large. The benefits of
this overarching paradigm are discussed in Sec. VIII.

##  IV Vector-encoding via Matrix Product States

Loading classical data into a quantum computer poses a significant challenge
in quantum computing [23]. This step is crucial in many quantum algorithms and
its implementation for arbitrary states is very complicated and inefficient
due to its exponential scaling with the number of qubits [36].

Nevertheless, preparation of an MPS can be done efficiently, i.e. with linear
scaling in the number of qubits [19]. Various other algorithms are available
for MPS preparation, e.g. Refs. [20, 21, 37, 38], which are quite similar and
all work well due to the naturalness of encoding these states into quantum
circuits. For a thorough analysis, please see our previous article, Ref. [23].
In that paper, we also proposed a new vector-encoding scheme based on
variational hardware-efficient circuits. By simulating these circuits with
tensor networks, this scheme enables variational gates to be optimized
according to calculations of Riemannian gradients. This provides considerably
better accuracy and speed. For example, we numerically simulated how vectors
sampled from analytical functions (such as polynomial, trigonometric, and
Gaussian distributions) could be efficiently encoded into shallow quantum
circuits for up to 100 qubits.

While much of the research focuses on the approximate preparation of a desired
quantum state, a quantum scheme that exactly prepares an arbitrary MPS of rank
$r$ is presented in Ref. [19]. According to this work, any MPS of rank $r$ can
be encoded into a step quantum circuit (as shown in Fig. 5), where each gate
acts on $[\log(r)]+1$ qubits, where $[\cdot]$ is the ceiling function. To
create such a circuit, the MPS cores are first transformed using QR
decomposition into isometric cores [39]. These isometries are then extended to
unitaries projected on zero-state vectors. The resulting network of unitary
matrices and zero states form the quantum circuit for preparing the MPS.

Figure 5: The quantum circuit for encoding an MPS of rank $2$ involves the
sequential application of two-qubit operations. This scheme is extendable to
the case of arbitrary rank.

##  V Matrix-encoding via Matrix Product Operators

As depicted in Fig. 6, an MPO is a representation of a matrix in TT format.
Since an MPS can be prepared quite effectively on a quantum computer, the
question arises about the same procedure for an MPO. Only one publication has
recently emerged dedicated to this particular topic [40] (see Sec.VIII for a
discussion). In this section, we propose our algorithm for effectively
encoding MPOs into quantum circuits.

Figure 6: The representation of an $N\times N$ matrix, $M$, as an MPO of rank
$r$. This representation occupies $O(\log_{2}(N)r^{2})$ memory, which is
significantly less than the initially required $N^{2}$ memory in the case of
small rank $r$.

###  V.1 Method

The main step in the exact preparation of an MPS is its orthogonalization.
However, for MPOs, there is no such algorithm; it is even not clear whether an
arbitrary MPO can be decomposed into an MPO with unitary and isometric cores:

$A=A_{s_{1}s_{2}\dots s_{n}}^{l_{1}l_{2}\dots
l_{n}}=V\text{\small{[1]}}^{s_{1}l_{1}}_{\alpha_{0}\alpha_{1}}V\text{\small{[2]}}^{s_{2}l_{2}}_{\alpha_{1}\alpha_{2}}\dots
V\text{\small{[n]}}^{s_{n}l_{n}}_{\alpha_{n-1}\alpha_{n}},$ (3)

where $\alpha_{0}=\alpha_{n}=1$. Here we refer to $V[k]$ as a unitary
(isometric) core, signifying that it constitutes a unitary (isometric) matrix
upon reshaping, aligning with the index division shown in Fig. 7.

Figure 7: (a) A unitary MPO. The dashed lines indicate the index partition by
which the cores should be reshaped to meet the unitarity (isometric)
condition. After the reshaping, each tensor
$V[k]^{s_{k}l_{k}}_{\alpha_{k-1}\alpha_{k}}$ is a matrix
$V[k]_{\left(l_{k}\alpha_{k}\right),(\alpha_{k-1}s_{k})}$ with row multi-index
$(l_{k}\alpha_{k})$ and column multi-index $(\alpha_{k-1}s_{k})$. The
unitarity condition for inner cores then reads as
$V^{*}[k]_{(\alpha_{k-1}^{\prime}s_{k}^{\prime}),(l_{k}\alpha_{k})}V[k]_{\left(l_{k}\alpha_{k}\right),(\alpha_{k-1}s_{k})}=\delta_{\alpha_{k-1}^{\prime}\alpha_{k-1}}\delta_{s_{k}^{\prime}s_{k}}$,
where the asterisk denotes complex conjugation and $\delta_{jj^{\prime}}$ is
the Kronecker delta. For simplicity of notation, we denote this unitarity as
condition $V^{\dagger}[k]V[k]=I$, meaning the matrix form of $V[k]$. (b) The
graphical notation of the unitarity condition $V^{\dagger}[k]V[k]=I$ for one
of the inner cores. For inner cores, the matrix forms $V[k]$ ($k$ = 2, …, n-1)
are unitary matrices, whereas $V[1]$ and $V^{\dagger}[n]$ are isometries and
satisfy the conditions $V^{\dagger}[1]V[1]=I$ and $V[n]V^{\dagger}[n]=I$,
respectively. (c) Graphical notation for the isometric condition
$V[n]V^{\dagger}[n]=I$. The right-hand sides of the equalities in (b) and (c)
are the identity matrices.

For inner cores, this unitarity condition reads as

$V^{*}[k]_{(\alpha_{k-1}^{\prime}s_{k}^{\prime}),(l_{k}\alpha_{k})}V[k]_{\left(l_{k}\alpha_{k}\right),(\alpha_{k-1}s_{k})}=\delta_{\alpha_{k-1}^{\prime}\alpha_{k-1}}\delta_{s_{k}^{\prime}s_{k}},$
(4)

where the asterisk denotes complex conjugation and $\delta_{jj^{\prime}}$ is
the Kronecker delta. To simplify the notation, we denote the unitarity
condition Eq. (4) as $V^{\dagger}[k]V[k]=I$, meaning the matrix form of
$V[k]$. Two boundary cores after reshaping form non-square isometric matrices
$V[1]$ and $V^{\dagger}[n]$, which satisfy the isometric conditions
$V^{\dagger}[1]V[1]=I$ and $V[n]V^{\dagger}[n]=I$, respectively. We refer to
an MPO with unitary and isometric cores as a unitary MPO.

Our main idea is to combine two key MPS-preparation methods: orthogonalization
and a variational approach. This lets us find an MPO with unitary cores that
is close to a given MPO such that encoding this unitary MPO into a quantum
circuit will be straightforward. All the steps of our algorithm are
illustrated in Fig. 8.

Figure 8: The proposed matrix-encoding algorithm. (a) A target $2^{n}\times
2^{n}$ matrix to be encoded. (b) Its rank-$r$ MPO representation with $n$
cores. (c) The higher rank-$R$ MPO approximation with unitary cores (unitary
MPO) achieved via Riemannian optimization. (d) The quantum circuit obtained
straightforwardly from the unitary MPO using auxiliary qubits and
measurements. (e) Overall, the resulting quantum circuit encodes the target
matrix and can be implemented in any quantum algorithm. The matrices that are
responsible for operations in the algorithm are denoted
$M_{1},M_{2},\ldots,M_{p}$. The simplest example is when the algorithm is a
series of matrix-vector multiplications.

In the first step, we represent the target matrix $M$ in terms of an MPO with
rank $r$ before running our encoding algorithm. Next, we approximate the
resulting MPO $M_{\mathrm{mpo}}$ with a unitary MPO of a higher rank, $R$. To
achieve this approximation, we employ a manifold gradient descent method,
optimizing the error in the Riemannian manifold of unitary cores [41, 42].
It’s important to note that to obtain unitary cores at the boundaries we
artificially utilize ancillas that effectively transform isometric cores into
unitary ones, and for such extension, we only need $[\log_{2}(R)]$ ancillas.

Subsequently, the resulting unitary cores are utilized in a quantum circuit as
multi-qubit gates, which are applied sequentially to the system and ancilla
qubits. After measuring the ancillas and post-selecting the outcomes, the
resulting quantum circuit implements the desired matrix. It is worth
emphasizing that our method allows for the implementation of arbitrary, not
necessarily unitary, operations, thanks to the post-selection subroutine.
Furthermore, the resulting quantum circuit can serve as a building block for
any quantum algorithm that requires matrix-vector multiplication by an MPO.

Figure 9: The scheme of encoding an MPO into a quantum circuit. (a) The matrix
to be encoded is represented as an MPO of rank $r$. (b) This rank-$r$ MPO is
approximated by a higher rank-$R$ unitary MPO using Riemannian optimization.
The green dashed lines indicate the index partition by which the cores must be
reshaped to form unitary matrices $U_{k}$ for $k=2,\dots,n-1$, and isometric
matrices for $k=1,n$ (see Fig. 7). (c) Two boundary isometric cores are
extended to unitary ones, $U_{1}$ and $U_{n}$. (d) The quantum scheme encoding
the target MPO is then produced by sequentially applying unitary gates
$U_{1},...,U_{n}$ to all system qubits and ancillary qubits in zero state,
followed by measurement of the ancillas and post-selection.

Since any matrix can be approximated by an MPO with some error, our algorithm
can approximately encode any matrix into a quantum computer. In addition, in
many applications, such as partial differential equations, quantum physics, or
quantum chemistry, an analytical expression of a matrix in MPO format is
readily available [43, 44, 45].

####  V.1.1 Preprocessing of matrices using the unitary MPO

To begin with, it is crucial to address a significant aspect of the
normalization of the matrix we want to encode. When implementing the action of
a matrix, $M$, on a state represented by a density matrix, $\rho$, on a
quantum computer, the resulting transformed state becomes $M\rho M^{\dagger}$.
However, it is essential to consider the constraint $\mathrm{Tr}[M\rho
M^{\dagger}]\leq 1$, as quantum operations do not increase the trace [46].
This constraint implies that the absolute values of the elements of matrix $M$
must not exceed a certain upper bound.

To overcome this limitation, we introduce a real normalization constant, $c$,
and our objective shifts towards realizing the $c^{-1}M$ matrix on the quantum
computer instead of just $M$. From a practical standpoint, this implies the
necessity of remembering this normalization factor and incorporating it into
the final calculations. By incorporating the normalization constant, $c$, we
ensure that the quantum computer can efficiently handle the transformation
$c^{-2}M\rho M^{\dagger}$, satisfying the trace constraint while preserving
the desired operation. It’s important to emphasize that considering this
normalization factor is merely a technical requirement, which can be
straightforwardly managed in the computations.

As mentioned above, our first step is to find a unitary MPO, $A$, and a
normalization constant, $c$, such that $cA$ approximates $M_{\mathrm{mpo}}$ as
closely as possible. The TT-rank, $R$, of $A$ should be taken to be larger
than the TT-rank, $r$, of $M_{\mathrm{mpo}}$ since unitarity takes away some
of the degrees of freedom. We formulate this approximation task in terms of
finding the minimum of the cost function:

$\displaystyle\min_{c,A}\hskip 10.0pt$ $\displaystyle\hbox{C}=||cA-
M_{\mathrm{mpo}}||^{2}$ subject to $\displaystyle
V[k]:V^{\dagger}[k]V[k]=I\mbox{ }\forall k=1,\dots,n-1,$ $\displaystyle\hskip
28.0ptV[n]V^{\dagger}[n]=I,$

where $||\cdot||$ denotes the Frobenius norm. Here, the optimum is taken over
both $c$ and $A$. We solve this optimization problem by alternating descent,
and the optimization algorithm (steps (a)-(b) in Fig. 9) reads as follows:

  1. Step 0.

Fix an initial approximation $A$.

  2. Step 1.

Find $c$ that minimizes the cost function $||cA-M_{\mathrm{mpo}}||^{2}$ with
the fixed $A$ via formula

$c=\real\frac{\mathrm{Tr}[A^{\dagger}M_{\mathrm{mpo}}]}{||A||^{2}}.$ (5)

  3. Step 2.

Find a new approximation $A$ via one step of the Riemannian gradient descent
for the problem

$\displaystyle\min\limits_{A}\hskip 10.0pt$ $\displaystyle C=||cA-
M_{\mathrm{mpo}}||^{2}$ (6a) subject to $\displaystyle
V[k]:V^{\dagger}[k]V[k]=I\mbox{ }\forall k=1,\dots,n-1,$ (6b)
$\displaystyle\hskip 28.0ptV[n]V^{\dagger}[n]=I,$ (6c)

and go to Step 1.

Step 2 is the most challenging because it involves optimization with
constraints (Eq. (6b), (6c)) that ensure the isometry of the MPO cores. Any
$m\times p$ $(m\geq p)$ isometric matrix belongs to the Stiefel manifold
$\text{St}(m,p,\mathbb{C}):=\\{V\in\mathbb{C}^{m\times p}|V^{\dagger}V=I\\}$,
which is a Riemannian manifold. The cost function Eq. (6a) is thus defined on
the Cartesian product of $n$ Stiefel manifolds. To find a set of isometric
matrices $V[k]$ minimizing this cost function, one can use gradient-based
optimization methods, such as Riemannian gradient descent which works for
optimization problems with constraints [47].

In regular gradient-based methods, such as gradient descent, the variables
being optimized are updated by moving in the opposite direction to the cost
function gradient: $V[k]\rightarrow V[k]-\alpha\partial_{k}C$, where $\alpha$
is a step size and $\partial_{k}C=\partial C/\partial V[k]$. In the Riemannian
generalization of these methods, to determine the direction of this movement,
instead of $\partial_{k}C$, its projection
$\partial_{k}^{R}C\in\mathcal{T}_{V[k]}$ onto the space tangent to the Stiefel
manifold at point $V[k]$ is used. $\partial_{k}^{R}$ is called the Riemannian
gradient, and to update the values of matrices $V[k]$, a retraction onto the
manifold is used:
$V[k]\rightarrow\mathcal{R}_{V[k]}(-\alpha\partial_{k}^{R}C)$ [41, 42].
Besides conventional gradient descent, there are Riemannian generalizations of
other, more advanced optimization methods, such as the adaptive moment
estimation algorithm (ADAM) [48, 49] that we use in this work.

To sum up, Riemannian gradient descent is an iterative process where, at each
step, the algorithm calculates the Riemannian gradient of the current point’s
cost function. Subsequently, it employs the retraction in the opposite
direction to determine the next point in the optimization procedure. Since
retraction along the manifold is employed, the updated optimization point
inherently belongs to the manifold. In our case, this means that conditions
Eq. (6b) and (6c) are automatically satisfied.

Once we get a unitary MPO, $A$, that is the approximation of
$c^{-1}M_{\mathrm{mpo}}$, all the inner cores of $A$ are already unitary
matrices (after the corresponding reshape). However, the tensors at the
boundaries ($V[1]$ and $V[n]$) are only isometries and they need to be
expanded so that they become unitary. Therefore, we extend them to unitary
matrices and project them onto $\ket{0}$ states to maintain the correspondence
with the matrix, see Fig. 9(c). Now all the cores are unitary and the
transition to the quantum circuit is straightforward.

####  V.1.2 Encoding using quantum circuit

In the context of quantum circuits, the $\ket{0}_{a_{1}}$ state (see Fig. 9
(c)) corresponds to preparing additional ancilla qubits in the zero state. At
the same time, the $\ket{0}_{a_{2}}$ state corresponds to the projection onto
the zero states, i.e. measurements with post-selection. Thus, we obtain the
following steps for the quantum part:

  1. 1.

Prepare ancillary qubits in the zero state $\ket{0}^{\otimes[\log_{2}(R)]}$
and leave the other qubits alone since our goal is only to encode the matrix.
We will receive an initial state as input to which we will then apply this
matrix.

  2. 2.

Sequentially apply gates $U_{1},...,U_{n}$, which should be decomposed into
one-qubit and two-qubit operations (see below).

  3. 3.

Measure the ancillary qubits in the computational basis. Note that one should
not measure the qubits that were originally introduced as ancillary, but the
qubits from the opposite end; see Fig. 9 (d).

If all measurement outcomes are equal to $0$, stop the algorithm. Else go back
to step 1.

As stated, the multi-qubit gates $U_{1},...,U_{n}$ must be decomposed into
one-qubit and two-qubit operations. Each $U_{k}$ acts on $([\log_{2}(R)]+1)$
qubits, where $R$ is the TT-rank of the unitary MPO, $A$. Ref. [14] proposes a
numerical strategy to approximately decompose a general multi-qubit unitary
with a CNOT gate count very close to the theoretical lower limit presented in
Ref. [11]. According to Ref. [11], decomposing an arbitrary $m$-qubit gate
requires $\simeq 4^{m-1}$ CNOT gates. In our case, we have to use $\simeq
n\cdot 4^{\log_{2}R}=nR^{2}$ CNOT gates in total. Therefore, the lower the
rank we take for the approximation, the simpler the quantum circuit turns out,
but the worse the accuracy of the approximation.

####  V.1.3 Post-selection success probability analysis

The introduction of ancillary qubits makes it possible to implement non-
unitary matrices on a quantum computer. The post-selection of ancilla
measurements makes the encoding algorithm probabilistic. In this section, we
examine the probability of successful encoding.

The quantum circuit $\mathbf{U}$ (see Fig. 9) implements the action of the
unitary MPO $A$ only if the ancillaries are all initialized and measured in
the state $\ket{0}$. We refer to the probability of achieving this successful
measurement as the success probability.

The success probability depends on the initial state
$\ket{\psi_{\mathrm{in}}}$ to which the quantum circuit is applied and can be
expressed as follows:

${\Pr}_{\text{suc}}=\mathrm{Tr}\big{[}\bra{0}_{a_{2}}\mathbf{U}\outerproduct{\psi_{\mathrm{in}}}{\psi_{\mathrm{in}}}\otimes\outerproduct{0}{0}_{a_{1}}\mathbf{U}^{\dagger}\ket{0}_{a_{2}}\big{]},$
(7)

Using $\bra{0}_{a_{2}}\mathbf{U}\ket{0}_{a_{1}}=A$, we can express the
probability as:

${\Pr}_{\text{suc}}=\mathrm{Tr}\big{[}A^{\dagger}A\outerproduct{\psi_{\mathrm{in}}}{\psi_{\mathrm{in}}}\big{]}=||A\ket{\psi_{\mathrm{in}}}||^{2}.$
(8)

We can represent the matrix $A^{\dagger}A$ in terms of its eigenvalues and
eigenstates:
$A^{\dagger}A=\sum_{k}\lambda_{k}^{2}(A)\ket{\lambda_{k}(A)}\bra{\lambda_{k}(A)}$.
In this case, the eigenvalues of $A^{\dagger}A$ are the squared singular
values $\lambda_{k}(A)$ of the matrix $A$. Thus the success probability reads

${\Pr}_{\text{suc}}=\sum\limits_{k}\lambda_{k}^{2}(A)\left|\innerproduct{\lambda_{k}(A)}{\psi_{\mathrm{in}}}\right|^{2}.$
(9)

Therefore, the success probability depends on the singular value spectrum of
$A$, as well as the initial state $\ket{\psi_{\mathrm{in}}}$ and its overlap
with the eigenvectors of $A^{\dagger}A$.

In certain scenarios, however, the initial state does not affect the success
probability. This happens when all the singular values are identical. One such
example is solving partial differential equations (PDEs) in closed quantum
systems, where encoded matrices are unitary with singular values all equal to
one [50].

According to Eq. (9), the maximum success probability for a fixed $A$ is
achieved when $\ket{\psi_{\mathrm{in}}}=\ket{\lambda_{\max}(A)}$:

$\max\limits_{\psi_{\mathrm{in}}}{\Pr}_{\text{suc}}=\lambda_{\max}^{2}(A)\leq
1,$ (10)

which also implies a limitation on the $A$’s singular values:
$\lambda_{k}(A)\leq 1$ for all $k$.

The above reasoning is general and also applies to other matrix encoding
methods that use auxiliary qubits. In our algorithm, $A$ is a unitary MPO
approximation of the normalized original matrix $M$. In the case of an exact
approximation, $A=c^{-1}M$ and
$\Pr_{\mathrm{suc}}=c^{-2}\sum_{k}\lambda_{k}^{2}(M)\left|\innerproduct{\lambda_{k}(M)}{\psi_{\mathrm{in}}}\right|^{2}$.
Then, taking into account Eq. (10), the probability of success for each fixed
$\ket{\psi_{\mathrm{in}}}$ reaches its maximum if $c=\lambda_{\max}(M)$.

Thus, a good strategy to increase the success probability is to normalize the
original matrix $M$ so that its maximum singular value is one, and then
approximate it with a unitary MPO for subsequent encoding into a quantum
circuit. However, finding the maximum singular value in the general case is a
difficult task for high-dimensional matrices. Therefore, we propose to
optimize the normalization constant when finding an approximation of the
target matrix by a unitary MPO.

####  V.1.4 The average success probability

As stated in the previous section, the initial state
$\ket{\psi_{\mathrm{in}}}$ to which the quantum circuit is applied can affect
the success probability. The initial state depends on the particular task. To
study the success probability without considering specific initial states, we
can uniformly average over all possible initial states of the system.

The average of uniformly distributed pure $n$-qubit quantum states with
respect to the Haar measure [51] equals $\frac{1}{2^{n}}I$, where $I$ is the
identity matrix. Thus, according to the Eq. (8) the average success
probability reads

$\langle{\Pr}_{\text{suc}}\rangle=\frac{||A||^{2}}{2^{n}},$ (11)

where $||\cdot||$ denotes the Frobenius norm. This equation shows that the
Frobenius norm of the unitary MPO, which depends on its singular value
spectrum, determines the average success probability. At first glance, it
might seem that $\langle{\Pr}_{\text{suc}}\rangle$ drops exponentially with
the number of qubits. However, this is not the general case, since $||A||^{2}$
can grow exponentially with the number of qubits. For instance, if $A$ is a
unitary matrix, then $||A||^{2}=2^{n}$ and the average success probability is
$1$. In Sec. V.2, we numerically simulate several commonly encountered
matrices and find that there is no exponential decrease in the average success
probability.

Figure 10: The minimum relative error, $\varepsilon$, achieved in the
approximation of (a) the diagonal and (b) the Laplace matrices by unitary MPOs
depending on its rank, $R$, for six different matrix sizes $2^{n}\times
2^{n}$, where $n$ is the number of qubits. The corresponding average success
probabilities for (c) the diagonal and (d) the Laplace matrices. For
convenience, we denote the abscissa axis as the number of used ancillas equal
to $\log_{2}R$. Examples of the matrix approximation for five qubits (i.e.,
the matrix sizes are $2^{5}\times 2^{5}$) with errors $3\cdot 10^{-5}$ for (e)
the diagonal matrix and $2\cdot 10^{-15}$ for (f) the Laplace matrix. The
relative error matrix is calculated using the formula $||M||^{-1}|cA-M|$.

###  V.2 Results

####  V.2.1 Diagonal and Laplace matrices

To demonstrate the capabilities of our method, we consider matrices that are
frequently used in computational mathematics and encode them into multi-qubit
quantum circuits. In this work, we consider a matrix, $M_{D}$, with a
discretized linear function on the diagonal and the Laplace operator, $M_{L}$,
corresponding to the second derivative matrix after finite difference
discretization [52], both of which are not unitary:

$M_{D}=\begin{pmatrix}0\\\ &1\\\ &&\ddots\\\ &&&2^{n}-1\end{pmatrix};\hskip
5.0ptM_{L}=\begin{pmatrix}2&-1\\\ -1&2&\ddots\\\ &\ddots&\ddots&-1\\\
&&-1&2\end{pmatrix}.$ (12)

The Laplace operator, $M_{L}$, is necessary for the solution of PDEs [53, 54]
while the diagonal matrix, $M_{D}$, as we describe in the Applications Sec.
VII, can be used in both optimization [55] and PDEs. For any $n>1$, the MPO
forms of the $M_{D}$ and $M_{L}$ matrices have TT-ranks $r$ equal to 2 and 3,
respectively. To investigate the performance of our method, we conducted the
following experiment for both matrices: fixing the number of qubits, $n$, and
the rank of the unitary MPO, $R$, we estimate the maximum achievable
approximation accuracy, $\varepsilon=||cA-M||^{2}/||M||^{2}$, and
corresponding average success probability defined in Eq (11).

The results of this experiment are presented in Fig. 10 – the minimum relative
error achieved during the approximation of these matrices by the unitary MPOs
is shown in the top plots and the corresponding success probability is shown
in the middle plots. For both matrices, our method can achieve reasonable
accuracy (i.e., errors of about $0.01\%$ and $0.5\%$ for $M_{D}$ and $M_{L}$,
respectively) and a sufficient success probability (i.e., more than $10-20\%$
up to 50 qubits).

Importantly, the error and the success probability almost do not depend on the
number of qubits. In the bottom part of Fig. 10, we show examples of the
resulting approximations and their comparison with exact matrices.

####  V.2.2 Approximation of multi-controlled Toffoli gate

Our scheme can be also applied to decompose multi-controlled quantum gates,
which are widely used in quantum computing. As an example, in this work we
consider the decomposition of multi-controlled Toffoli gate (MCT), a crucial
matrix in quantum computing that finds application in various quantum
algorithms, including Grover’s algorithm [56], Shor’s factorization algorithm
[57], and solving a linear system of equations [58, 59]. Therefore, the
efficient decomposition of this gate is one of the key tasks of quantum
computing [60, 61, 62].

The MCT gate for $n$ qubits can be represented in matrix form as follows:

$M_{\text{MCT}}=\begin{pmatrix}1\\\ &\ddots\\\ &&1\\\ &&&0&1\\\ &&&1&0\\\
\end{pmatrix}.$ (13)

The important feature of this matrix is that its TT-rank in the MPO
representation equals $2$, which allows us to successfully employ our encoding
method. It is possible to approximate the MCT gate with high accuracy by a
unitary MPO with sufficient rank $R$, achieving an error $\varepsilon$ less
than $10^{-9}$. An illustrative example of such an approximation for a
$5$-qubit system is presented in Fig. 11.

Figure 11: Example of the multi-controlled Toffoli gate approximation by a
unitary MPO for $5$ qubits with an error $\varepsilon=6\cdot 10^{-10}$. The
relative error matrix is calculated using the formula $||M||^{-1}|cA-M|$,
where $M$ is a MCT gate matrix (13). Figure 12: (a) The dependence of the MPO
representation rank $r$ of the matrix $e^{-iHdt}$ with respect to $dt$. At
small $dt$, the rank $r$ weakly depends on the number of qubits $n$. However,
as $dt$ increases, the rank $r$ plateaus and reaches a maximum value that is
solely based on the number of qubits and equals $4^{[\frac{n}{2}]}$. The
rounding tolerance for the MPO approximation is $10^{-6}$. (b) The
relationship between $dt$ and the minimum achieved error $\varepsilon$ in the
approximation of $e^{-iHdt}$ by a unitary MPO of rank $R$ for 3 different
values of $R$ (solid) and by Trotter product formulae of 1st and 2nd order
(dashed). The number of qubits equals $8$ in all cases. (c) An example of the
evolution operator approximation for $5$ qubits and $t=0.05$ with an error
$2\cdot 10^{-8}$. The relative error matrix is calculated using the formula
$||M||^{-1}|cA-M|$.

####  V.2.3 Hamiltonian simulation

One of the primary areas of use for quantum computation is the simulation of
Hamiltonian dynamics. The challenge of simulating the behavior of quantum
systems was the original inspiration for quantum computers, and it still
serves as one of their major potential applications. Hamiltonian simulation
acts as the fundamental approach for simulating quantum systems and is also
applicable to quantum algorithms and quantum chemistry.

To demonstrate how well our algorithm can perform in approximating an
evolution operator, we pick an open quantum Ising chain with both transverse
and longitudinal fields for the Hamiltonian:

$H=J\sum\limits_{k=1}^{n-1}Z_{k}Z_{k+1}+g\sum\limits_{k=1}^{n}X_{k}+h\sum\limits_{k=1}^{n}Z_{k}.$
(14)

Here, $X$ and $Z$ are Pauli matrices and a non-integrable point is chosen for
the Hamiltonian parameters: $J=2$, $g=1$, and $h=1$.

Our task is to encode the evolution operator $e^{-iHdt}$ into a quantum
circuit with a small value of $dt$, for the evolution operator to have small
ranks; see Fig. 12(a). Then, by repeating this circuit several times, we can
perform the evolution at any desired time $T$. In Fig. 12(b), we present the
results of the comparison in terms of the accuracy of our algorithm with the
Trotter decomposition that was used in a recent quantum utility work [63]. It
can be seen that our algorithm demonstrates better accuracy when choosing a
sufficient rank $R$ at reasonable $dt$ values.

##  VI Readout via MPS tomography

The last important step of any quantum algorithm is a readout of the output
state. In general, to implement efficient tomography of $n$-qubit quantum
states, it is necessary to perform an exponential (in $n$) number of
measurements, and an exponential amount of classical memory and computing
power is required. To overcome this, various approaches have been introduced:
classical shadow [64], neural network methods [65], and TNs [31, 66, 30].
Here, we propose to use TN-based tomography [31] for several reasons:

  * •

It provides a more comprehensive representation of the vector in the MPS
format, extending beyond the mere capability to measure observables [64].

  * •

Having the output in the TN format serves as a valuable input for subsequent
stages of a hybrid algorithm. Our scheme, for instance, involves the
progression from MPS to quantum computing and back to MPS.

  * •

TNs offer a more transparent and explainable approach compared to alternatives
such as neural networks [67, 68]

###  VI.1 Method

In our experiments, we employed MPS tomography as outlined in Ref. [30]. The
algorithm can be summarized as follows:

  * •

Firstly, choose the initial MPS approximation with fixed rank $r$ for the
output state. Note that with this method, complex-valued vectors can also be
measured.

  * •

Secondly, measure the output state in different randomly chosen bases (for
example, XXYZ…ZXY), where the outcomes, $x_{i}$, in bases, $b$, can be denoted
as $x_{i}^{(b)}$.

  * •

Next, introduce a loss function to quantify the relative entropy between the
predictions from the MPS approximation and the actual measured probabilities.

$\text{Loss}=-\frac{1}{N}\sum_{i=1}^{N}\log P(x_{i}^{(b)}).$ (15)

Here, $P(x_{i}^{(b)})$ represents the probability of obtaining the measured
outcome $x_{i}$ in basis $b$, with $N$ denoting the number of measurements. We
obtain this probability by projecting the MPS at the current iteration onto
the corresponding basis vector, which is done effectively in the format of
tensor networks.

  * •

Subsequently, employ stochastic gradient descent on the parameters of the MPS
cores to minimize the loss as much as possible [69].

###  VI.2 Results

The key aspect of the above method lies in its ability to maintain a linear
scaling of the number of measurements (or samples) with the number of qubits
for tomographing a state while ensuring a specified level of accuracy. To
explore this, we examine the MPS tomography of a variational quantum circuit
(VQC) output state, typically well-described in the form of MPS when the
number of layers is limited [70].

Fig. 13 depicts an example of a variational quantum circuit and the results of
the MPS tomography for such a circuit. Notably, we observe a linear
correlation between the number of measurements and the number of qubits when
the number of layers is fixed.

Figure 13: (a) A VQC example. The circuit is composed of multiple layers,
while each layer is a sequence of non-parametric CNOT gates and single-qubit
variational gates $u$. The parameters of the variational gates $u$ are
determined at random for our experiment. (b) MPS tomography results for the
output states of random VQCs. The plot shows the number of measurements
required to obtain 99% quantum fidelity, depending on the number of qubits $n$
in the circuit with a fixed numbers of layers, $L$.

##  VII Applications

At first glance, it may seem that this framework does not find a lot of
applications. However, in this section, we consider several examples. In
essence, what a quantum computer can do is sequentially perform Matvec
operations:

$X=U_{n}U_{n-1}\dots U_{2}U_{1}X_{0},$

where $X_{0}$ is the initial state-vector and $U_{i}$ are matrices (not
necessarily unitary). Therefore, regarding quantum algorithms, encoding
matrices into quantum circuits is one of the most important tasks. This is
precisely the task that we provide a solution to in this article - the main
feature is that the number of gates in the proposed circuit has a linear
scaling with the number of qubits and it is possible to control the size of
the circuit, varying the accuracy of the approximation. In contrast, the
encoding of an arbitrary $n$-qubit matrix scales exponentially [11].

To see how a quantum implementation can possess advantages over a classical
one, consider the following. If the task is to implement a matrix-vector
multiplication, $Ax$, then the classical realization via tensor networks has
complexity $O(ndr^{3}(A)r^{3}(x))$, where $n^{d}$ is the vector size, and
$r(A)$ and $r(x)$ are TT-ranks of $A$ and $x$, respectively. The complexity of
the quantum implementation does not depend on the rank $r(x)$, only on the
rank $r(A)$. Therefore, when there are large solution ranks and small
operation ranks, the quantum computer provides advantages. As will be shown
below, this approach has numerous applications, including the solution of
PDEs, optimization, machine learning, and quantum chemistry.

###  VII.1 Simulation

When almost any differential equation is discretized in time and space, it is
reduced to an iterative scheme, which in its simplest form is written as
follows:

$x_{i+1}=Mx_{i},$

where $M$ represents the derivative operator. This matrix $M$ is usually
diagonal or contains first and second derivatives, and is therefore well-
represented by an MPO with small TT-ranks [17, 28, 32]. Often matrix $M$ does
not depend on the iteration number, $i$, and is fixed. In this case, the
implementation of our numerical scheme on a quantum computer is advantageous
because it can be realized in a way that is independent of the ranks of the
solution, $x_{i}$, during the iterations. Consequently, increasing the ranks
during the iterations, which is often the case in such problems [32], does not
affect the computational complexity on a quantum computer - it will only
depend on the efficiency of encoding matrix $M$ into a quantum circuit.

For example, if we want to solve the heat-diffusion equation:

$\displaystyle\frac{\partial u}{\partial t}$ $\displaystyle=k^{2}\Delta
u,\quad x\in[0,1],t\geq 0;$ $\displaystyle u(x,0)$
$\displaystyle=u_{0}(x),\quad x\in[0,1];$ $\displaystyle u(0,t)$
$\displaystyle=u(1,t)=0,\quad t\geq 0;$

after discretization (explicit in time), it is reduced to

$\displaystyle u_{i+1}$ $\displaystyle=(I+\frac{\Delta t}{\Delta
x^{2}}k^{2}M_{L})u_{i},$ $\displaystyle u_{0}$ $\displaystyle=\hat{u}_{0},$

where $M_{L}$ is the Laplace operator from Section V.2. Here, $\Delta t$ and
$\Delta x$ are discretization steps in time and space, respectively, $u_{i}$
is the discretized solution at time $i\Delta t$. If the initial ranks of
$u_{0}$ are large then this scheme is better suited to a quantum
implementation than a classical one.

###  VII.2 Optimization

Tensor trains are used to tackle black-box optimization problems [71, 72, 73].
Remarkably, it is often possible to reduce an optimization problem to the task
of finding the maximum element in a tensor $y$, regardless of whether it is a
discrete, continuous, or black box optimization (including minimization) task
[74]. In this formulation, one can always apply the so-called Power method
[75]. That is, finding the maximal element in the normalized tensor $y^{n}$
instead of $y$, which is much easier since the largest element in $y^{n}$ is
much larger, in relative terms, compared to the largest element of $y$.
Ideally, with a sufficiently large degree, $n$, the normalized $y^{n}$ turns
into a delta function in a peak with a maximum of $y$.

We can approximate $y$ as a TT, $y_{TT}$. The overall accuracy of this
approximation can be poor as long as the maximal element of $y_{TT}$ is close
to the maximal element of $y$. We can now create a quantum circuit that
realizes a diagonal MPO matrix, $M_{y}$, that corresponds to $y_{TT}$ using
our algorithm. Thus, we can implement the Power method on a quantum computer,
preparing the initial state as a vector of units:

$y^{n}\simeq M_{y}M_{y}...M_{y}M_{y}\frac{1}{\sqrt{N}}\begin{pmatrix}1\\\ 1\\\
\vdots\\\ 1\end{pmatrix}.$

It is important to note that when implemented on a quantum computer, the
normalization will be automatically observed, due to post-selection (see Sec.
V.1.3). This approach yields similar quantum advantages as in the case of
solving PDEs [76]. That is when implementing the power method on tensor
networks, the ranks can grow significantly, while for a quantum computer, the
complexity will not depend on the ranks $y$ during the multiplications.

###  VII.3 Machine Learning

It is possible to consider arbitrary neural networks (convolutional neural
networks, transformers, etc.) in the form of tensor networks [77, 78, 79]. The
contraction of these tensor networks would be computationally challenging on a
classical computer due to the presence of large bond dimensions [24]. However,
by encoding these tensor networks on a quantum computer, the calculations
could be performed efficiently regardless of the architecture’s complexity,
which is a promising avenue for quantum machine learning [23].

###  VII.4 Quantum Chemistry

The description of molecules involves an accurate simulation of electrons and
nuclei. The most common representation is the nonrelativistic molecular
Hamiltonian

$\hat{H}=\sum_{N}\hat{T}_{N}+\sum_{e}\hat{T}_{e}+\sum_{N,M}\frac{Z_{N}\cdot
Z_{M}}{\hat{r}_{NM}}-\sum_{N,e}\frac{Z_{N}}{\hat{r}_{Ne}}+\sum_{e,f}\frac{1}{\hat{r}_{ef}}.$
(16)

Here, $N,M$ denote nuclei, $e,f$ denote electrons, $\hat{T}$ denotes the
kinetic energy, $\hat{r}$ is the distance of two particles, and $Z$ is the
nuclear charge. In this formulation, both the electrons and the nuclei are
described as quantum particles and they interact via the Coulomb force. The
adiabatic or Born-Oppenheimer approximation splits this problem into weakly
coupled electron and nuclei sub-problems. While certain problems require
including relativistic effects or going beyond the Born-Oppenheimer
approximation, a large class of quantum chemical problems are solved this way.

After invoking the Born-Oppenheimer approximation, the problem is reduced to
the so-called electronic structure problem. In the second quantization
formulation, the corresponding Hamiltonian reads

$\hat{H}_{el}=\sum_{p,q}h_{pq}\cdot\hat{a}_{p}\hat{a}_{q}^{\dagger}+\sum_{p,q,r,s}h_{pqrs}\cdot\hat{a}_{p}\hat{a}_{q}\hat{a}_{r}^{\dagger}\hat{a}_{s}^{\dagger}.$
(17)

Here, $p,q,r,s$ are atomic orbital indices, $\hat{a},\hat{a}^{\dagger}$ are
fermionic annihilation and creation operators, respectively, $h_{pq}$ entails
the kinetic energy of the electrons plus the nuclear-electron attraction for
fixed nuclear positions, and $h_{pqrs}$ encodes electron-electron repulsion.
The Hamiltonian elements stem from the corresponding integrals over atomic
orbitals and are evaluated using recursion formulae [80]. The brute force
approach to solving the problem scales exponentially with the number of
electrons and orbitals. DMRG has been a powerful tool in treating the
electronic structure problem [81, 82]. In this approach, the electronic state
is described as a TT, and self-consistent field schemes are used to solve for
ground- and excited states.

The quantum tensor network formulation of this problem facilitates simulating
strong entanglement cases on a quantum computer: first, a DMRG computation is
performed to obtain a low-rank approximation to the system of interest. This
serves as the initial state. Subsequently, the MPO is represented using a
unitary MPO and mapped to a quantum circuit using the scheme described in the
present work.

The electron-electron correlation determines the required rank in DMRG methods
and varies widely between molecules, their geometries, and considered states.
For certain instances of the problem, the tensor train description exceeds
classical resources and sufficiently reliable quantum computers can provide an
advantage in solving these systems. While ground state chemistry of standard
organic molecules [83] is often described sufficiently accurately using
polynomially scaling methods such as coupled cluster and/or perturbation
theory, certain systems are beyond the capabilities of today’s classical
computers. Potential candidate systems include metal complexes, e.g. chromium
dimer, active sites of proteins with metal atoms in their center, and systems
showing nonadiabatic interactions in which multiple electronic states are
interacting. The effectiveness of quantum computers still depends on the
robustness of the quantum computers and certain TN structures might provide
more accurate solutions at lower ranks.

##  VIII Conclusion

This article introduces a novel numerical method for encoding MPOs into
quantum circuits with a depth that depends linearly on the number of qubits.
The devised algorithm offers a compelling solution to the crucial challenge of
efficiently encoding matrices into quantum systems. It enables the effective
encoding of a wide variety of matrices, not strictly confined to unitary
matrices, that can be represented in MPO form. Such matrices are prevalent in
diverse practical applications, ranging from solving partial differential
equations to optimization.

To test the algorithm’s applicability, we addressed the encoding of matrices
like the Laplace operator, a diagonal matrix, the multi-controlled Toffoli
gate, and the evolution operator of the Ising model Hamiltonian. Our findings
illustrate the high-precision encoding of these matrices onto quantum
circuits, demonstrating the algorithm’s efficiency for matrices that can act
on up to 50 qubits. This demonstrates the algorithm’s versatility and
potential to enhance quantum computational capabilities in diverse domains.
Notably, the proposed circuits exhibit a linear scaling of gate count with
qubit number, establishing the efficiency of our matrix-encoding algorithm in
comparison to existing multi-qubit gate decomposition methods, which have
exponential scaling [11, 12, 13].

An effective solution to the problem of encoding matrices into a quantum
circuit using TTs allowed us to introduce the concept of hybrid Tensor Quantum
Programming. This strategy integrates classical and quantum computing to solve
numerical problems framed within tensor networks. Starting with a classical
phase, computations are performed classically, while the TT-ranks are small,
benefiting from the efficiency of classical implementation. When the ranks
increase, indicating potential quantum benefits, the quantum phase is
initiated and the intermediate solution is transferred to a quantum computer.
During this phase, the rank growth is unlimited, and once the ranks decrease
again, the process reverts to the classical phase. This hybrid approach
outperforms purely classical methods by efficiently handling large-rank
problems, utilizing classical computation for small ranks and quantum
computation for large ranks. Consequently, it optimizes computational
resources for addressing problems with varying degrees of entanglement and
correlation.

It is important to note that while our MPO-encoding algorithm demonstrates
good applicability, particularly in scenarios featuring low-rank MPOs, it does
not completely solve the problem of accurately encoding arbitrary matrices
into quantum circuits. In the general case, an arbitrary matrix has full rank
in the MPO format, which can significantly complicate the encoding scheme.
However, our algorithm can be successfully used for encoding a low-rank MPO
approximation of the original matrix obtained by cutting the rank with the
introduction of an additional error.

The utilization of Riemannian optimization allows us to achieve high levels of
accuracy in our MPO-encoding algorithm. However, such optimization techniques
do not inherently guarantee the absolute accuracy of the obtained solution.
Using block-encoding for the MPO [40] provides completely accurate encoding
but requires the incorporation of an additional auxiliary qubit for each MPO
core. Consequently, the number of auxiliary qubits increases in proportion to
the system qubits, which limits the applicability of this method, particularly
for encoding large matrices. Riemannian optimization allows for more flexible
MPO-encoding, which is imprecise but requires fewer auxiliary qubits. Thus, a
promising path for future development may be exploring the trade-off between
the encoding accuracy, the auxiliary qubit number, and the probability of
successful implementation. This is achievable by improving the optimization
strategy in our algorithm, combining it with block-encoding, or developing
completely new encoding techniques.

The success probability in the MPO encoding algorithm is determined not only
by the intrinsic properties of the encoded matrix itself and the encoding
accuracy but also by the initial state vector on which the matrix acts. For
specific initial states, the success probability can significantly exceed the
average success probability. The highest success probability is achieved when
the initial state aligns with the encoded matrix’s right singular vector with
the maximum singular value. Such initial states prove important not only for
maximizing success probability but also from a practical standpoint, since
finding the largest singular value (or eigenvalue) states is a prevalent
challenge. In this case, following the Tensor Quantum Programming scheme, the
initial state is a classically precomputed approximate solution, which is
close to the largest singular-value state. This approach can greatly increase
the likelihood of successfully implementing the quantum part, finding utility,
for instance, in approaches like the optimization power method [75].

This work opens up directions for further optimal development of quantum
algorithms and computations. The driving idea is to find the drawbacks and
difficulties in tensor network algorithms and try to eliminate them with the
help of quantum computers, drawing on the substantial connection between both
approaches. The primary example is a significant increase in the ranks
inherent in numerical schemes, which can be minimized by the ability of a
quantum computer to operate on states with arbitrary entanglement (ranks).

Our concept of Tensor Quantum Programming is universal and can be generalized
to different tensor network architectures [67, 84]. It can be tailored to
particular problems and available quantum device architectures with specified
qubit connectivity. This adaptability allows the Tensor Quantum Programming
scheme to find potential applications in a range of fields coupled to tensor
networks, for instance, in machine learning. Further research will be required
to explore these avenues.

###### Acknowledgements.

We thank Alexey Melnikov, Dmitry Morozov, Akshay Vishwanathan, Vyacheslav
Kuzmin, and Valerii Vinokur for useful discussions.

## References

  * [1] Suguru Endo, Zhenyu Cai, Simon C Benjamin, and Xiao Yuan.  Hybrid quantum-classical algorithms and quantum error mitigation.  Journal of the Physical Society of Japan, 90(3):032001, 2021.
  * [2] Michael Perelshtein, Asel Sagingalieva, Karan Pinto, Vishal Shete, Alexey Pakhomchik, Artem Melnikov, Florian Neukart, Georg Gesek, Alexey Melnikov, and Valerii Vinokur.  Practical application-specific advantage through hybrid quantum computing.  arXiv preprint arXiv:2205.04858, 2022.
  * [3] M. Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C. Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R. McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, and Patrick J. Coles.  Variational quantum algorithms.  Nature Reviews Physics, 3(9):625–644, Aug 2021.
  * [4] Jacob Biamonte, Peter Wittek, Nicola Pancotti, Patrick Rebentrost, Nathan Wiebe, and Seth Lloyd.  Quantum machine learning.  Nature, 549(7671):195–202, 2017.
  * [5] William Huggins, Piyush Patil, Bradley Mitchell, K Birgitta Whaley, and E Miles Stoudenmire.  Towards quantum machine learning with tensor networks.  Quantum Science and Technology, 4(2):024001, 2019.
  * [6] Scott Aaronson.  Read the fine print.  Nature Physics, 11(4):291–293, 2015.
  * [7] Ewin Tang.  Dequantizing algorithms to understand quantum advantage in machine learning.  Nature Reviews Physics, 4(11):692–693, 2022.
  * [8] John M Martyn, Zane M Rossi, Andrew K Tan, and Isaac L Chuang.  Grand unification of quantum algorithms.  PRX Quantum, 2(4):040203, 2021.
  * [9] M. R. Perelshtein, A. I. Pakhomchik, A. A. Melnikov, A. A. Novikov, A. Glatz, G. S. Paraoanu, V. M. Vinokur, and G. B. Lesovik.  Solving large-scale linear systems of equations by a quantum hybrid algorithm.  Annalen der Physik, 534(7):2200082, May 2022.
  * [10] Leigh Lapworth.  A Hybrid Quantum-Classical CFD Methodology with Benchmark HHL Solutions.  arXiv preprint arXiv:2206.00419, 2022.
  * [11] Vivek V. Shende, Igor L. Markov, and Stephen S. Bullock.  Minimal universal two-qubit controlled-NOT-based circuits.  Phys. Rev. A, 69:062321, Jun 2004.
  * [12] V.V. Shende, S.S. Bullock, and I.L. Markov.  Synthesis of quantum-logic circuits.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 25(6):1000–1010, Jun 2006.
  * [13] Anna M. Krol, Aritra Sarkar, Imran Ashraf, Zaid Al-Ars, and Koen Bertels.  Efficient decomposition of unitary matrices in quantum circuit compilers.  Applied Sciences, 12(2), 2022.
  * [14] Péter Rakyta and Zoltán Zimborás.  Approaching the theoretical limit in quantum gate decomposition.  Quantum, 6:710, May 2022.
  * [15] Martin Plesch and Časlav Brukner.  Quantum-state preparation with universal gate decompositions.  Physical Review A, 83(3):032302, 2011.
  * [16] Matteo Paris and Jaroslav Rehacek.  Quantum state estimation, volume 649.  Springer Science & Business Media, 2004.
  * [17] Ivan V Oseledets.  Tensor-train decomposition.  SIAM Journal on Scientific Computing, 33(5):2295–2317, 2011.
  * [18] Román Orús.  A practical introduction to tensor networks: Matrix product states and projected entangled pair states.  Annals of Physics, 349:117–158, Oct 2014.
  * [19] C. Schön, E. Solano, F. Verstraete, J. I. Cirac, and M. M. Wolf.  Sequential Generation of Entangled Multiqubit States.  Phys. Rev. Lett., 95:110503, Sep 2005.
  * [20] Shi-Ju Ran.  Encoding of matrix product states into quantum circuits of one- and two-qubit gates.  Phys. Rev. A, 101:032310, Mar 2020.
  * [21] Peng-Fei Zhou, Rui Hong, and Shi-Ju Ran.  Automatically differentiable quantum circuit for many-qubit state preparation.  Phys. Rev. A, 104:042601, Oct 2021.
  * [22] Manuel S Rudolph, Jing Chen, Jacob Miller, Atithi Acharya, and Alejandro Perdomo-Ortiz.  Decomposition of matrix product states into shallow quantum circuits.  arXiv preprint arXiv:2209.00595, 2022.
  * [23] Alexey Melnikov, Mohammad Kordzanganeh, Alexander Alodjants, and Ray-Kuang Lee.  Quantum machine learning: from physics to software engineering.  Advances in Physics: X, 8(1):2165452, 2023.
  * [24] Igor L Markov and Yaoyun Shi.  Simulating quantum computation by contracting tensor networks.  SIAM Journal on Computing, 38(3):963–981, 2008.
  * [25] Salman Ahmadi-Asl, Maame Gyamfua Asante-Mensah, Andrzej Cichocki, Anh-Huy Phan, Ivan Oseledets, and Jun Wang.  Cross tensor approximation for image and video completion, 2022.
  * [26] Adam Holmes and AY Matsuura.  Efficient quantum circuits for accurate state preparation of smooth, differentiable functions.  In 2020 IEEE International Conference on Quantum Computing and Engineering (QCE), pages 169–179. IEEE, 2020.
  * [27] Ivan Oseledets.  Constructive representation of functions in low-rank tensor formats.  Constructive Approximation, 37, 09 2010.
  * [28] Vladimir A Kazeev and Boris N Khoromskij.  Low-rank explicit QTT representation of the Laplace operator and its inverse.  SIAM journal on matrix analysis and applications, 33(3):742–758, 2012.
  * [29] Artem A. Melnikov, Alena A. Termanova, Sergey V. Dolgov, Florian Neukart, and Michael R. Perelshtein.  QPrep: Library for quantum state preparation using tensor networks.
  * [30] Giacomo Torlai and Matthew Fishman.  PastaQ: A package for simulation, tomography and analysis of quantum computers, 2020.
  * [31] Marcus Cramer, Martin B Plenio, Steven T Flammia, Rolando Somma, David Gross, Stephen D Bartlett, Olivier Landon-Cardinal, David Poulin, and Yi-Kai Liu.  Efficient quantum state tomography.  Nature communications, 1(1):1–7, 2010.
  * [32] Egor Kornev, Sergey Dolgov, Karan Pinto, Markus Pflitsch, Michael Perelshtein, and Artem Melnikov.  Numerical solution of the incompressible navier-stokes equations for chemical mixers via quantum-inspired tensor train finite element method.  arXiv preprint arXiv:2305.10784, 2023.
  * [33] Madeleine Udell and Alex Townsend.  Why are big data matrices approximately low rank?  SIAM Journal on Mathematics of Data Science, 1(1):144–160, 2019\.
  * [34] Bernhard Beckermann and Alex Townsend.  Bounds on the singular values of matrices with displacement structure.  SIAM Review, 61(2):319–344, 2019.
  * [35] Jielun Chen, EM Stoudenmire, Steven R White, et al.  Quantum Fourier transform has small entanglement.  PRX Quantum, 4(4):040318, 2023.
  * [36] Martin Plesch and Časlav Brukner.  Quantum-state preparation with universal gate decompositions.  Phys. Rev. A, 83:032302, Mar 2011.
  * [37] Javier Gonzalez-Conde, Thomas W Watts, Pablo Rodriguez-Grasa, and Mikel Sanz.  Efficient amplitude encoding of polynomial functions into quantum computers.  arXiv preprint arXiv:2307.10917, 2023.
  * [38] Raghav Jumade and Nicolas PD Sawaya.  Data is often loadable in short depth: Quantum circuits from tensor networks for finance, images, fluids, and proteins.  arXiv preprint arXiv:2309.13108, 2023.
  * [39] Daniel Malz and Rahul Trivedi.  Computational complexity of isometric tensor network states, 2024.
  * [40] Martina Nibbi and Christian B. Mendl.  Block encoding of matrix product operators, 2023.
  * [41] Markus Hauru, Maarten Van Damme, and Jutho Haegeman.  Riemannian optimization of isometric tensor networks.  SciPost Phys., 10:040, 2021.
  * [42] Ilia A Luchnikov, Mikhail E Krechetov, and Sergey N Filippov.  Riemannian geometry and automatic differentiation for optimization problems of quantum physics and quantum technologies.  New Journal of Physics, 23(7):073006, Jul 2021.
  * [43] W. Hackbusch and B. N. Khoromskij.  Low-rank Kronecker-product approximation to multi-dimensional nonlocal operators. II. HKT representation of certain operators.  Computing, 76(3-4):203–225, 2006.
  * [44] L. Grasedyck.  Existence and computation of low Kronecker-rank approximations for large systems in tensor product structure.  Computing, 72:247–265, 2004.
  * [45] D. V. Savostyanov, S. V. Dolgov, J. M. Werner, and I. Kuprov.  Exact NMR simulation of protein-size spin systems using tensor train formalism.  Phys. Rev. B, 90:085139, 2014.
  * [46] Michael A Nielsen and Isaac Chuang.  Quantum computation and quantum information, 2002.
  * [47] Ilia A. Luchnikov, Alexander Ryzhov, Sergey N. Filippov, and Henni Ouerdane.  QGOpt: Riemannian optimization for quantum technologies.  SciPost Phys., 10:079, 2021.
  * [48] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [49] Gary Bécigneul and Octavian-Eugen Ganea.  Riemannian adaptive optimization methods.  arXiv preprint arXiv:1810.00760, 2018.
  * [50] Angel Rivas and Susana F Huelga.  Open quantum systems, volume 10.  Springer, 2012.
  * [51] Antonio Anna Mele.  Introduction to Haar Measure Tools in Quantum Information: A Beginner’s Tutorial, 2024.
  * [52] Christian Grossmann, Hans-Görg Roos, and Martin Stynes.  Numerical treatment of partial differential equations, volume 154\.  Springer, 2007.
  * [53] Hai-Ling Liu, Yu-Sen Wu, Lin-Chun Wan, Shi-Jie Pan, Su-Juan Qin, Fei Gao, and Qiao-Yan Wen.  Variational quantum algorithm for the Poisson equation.  Phys. Rev. A, 104:022418, Aug 2021.
  * [54] Shengbin Wang, Zhimin Wang, Wendong Li, Lixin Fan, Zhiqiang Wei, and Yongjian Gu.  Quantum fast Poisson solver: the algorithm and complete and modular circuit design.  Quantum Information Processing, 19(6):1–25, 2020.
  * [55] Lars Grasedyck, Lukas Juschka, and Christian Löbbert.  Finding entries of maximum absolute value in low-rank tensors.  arXiv preprint arXiv:1912.02072, 2019.
  * [56] Lov K Grover.  A fast quantum mechanical algorithm for database search.  In Proceedings of the twenty-eighth annual ACM symposium on Theory of computing, pages 212–219, 1996.
  * [57] Peter W Shor.  Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer.  SIAM review, 41(2):303–332, 1999.
  * [58] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  Physical Review Letters, 103(15), Oct 2009.
  * [59] Andrew M Childs, Robin Kothari, and Rolando D Somma.  Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  SIAM Journal on Computing, 46(6):1920–1950, 2017.
  * [60] Vivek V Shende and Igor L Markov.  On the CNOT-cost of Toffoli gates.  arXiv preprint arXiv:0803.2316, 2008.
  * [61] Amit Saha, Ritajit Majumdar, Debasri Saha, Amlan Chakrabarti, and Susmita Sur-Kolay.  Asymptotically improved Grover’s algorithm in any dimensional quantum system with novel decomposed n-qudit Toffoli gate.  arXiv preprint arXiv:2012.04447, 2020.
  * [62] A. S. Nikolaeva, E. O. Kiktenko, and A. K. Fedorov.  Decomposing the generalized Toffoli gate with qutrits.  Phys. Rev. A, 105:032621, Mar 2022.
  * [63] Youngseok Kim, Andrew Eddins, Sajant Anand, Ken Xuan Wei, Ewout Van Den Berg, Sami Rosenblatt, Hasan Nayfeh, Yantao Wu, Michael Zaletel, Kristan Temme, et al.  Evidence for the utility of quantum computing before fault tolerance.  Nature, 618(7965):500–505, 2023.
  * [64] Hsin-Yuan Huang, Richard Kueng, and John Preskill.  Predicting many properties of a quantum system from very few measurements.  Nature Physics, 16(10):1050–1057, 2020.
  * [65] Juan Carrasquilla, Giacomo Torlai, Roger G. Melko, and Leandro Aolita.  Reconstructing quantum states with generative models.  Nature Machine Intelligence, 1(3):155–161, Mar 2019.
  * [66] Jens Eisert.  Entanglement and tensor network states.  arXiv preprint arXiv:1308.3318, 2013.
  * [67] Valentin Khrulkov, Alexander Novikov, and Ivan Oseledets.  Expressive power of recurrent neural networks.  arXiv preprint arXiv:1711.00811, 2017.
  * [68] Ivan Oseledets and Eugene Tyrtyshnikov.  TT-cross approximation for multidimensional arrays.  Linear Algebra and its Applications, 432(1):70–88, 2010.
  * [69] Sebastian Ruder.  An overview of gradient descent optimization algorithms.  arXiv preprint arXiv:1609.04747, 2017.
  * [70] M. R. Perelshtein, A. I. Pakhomchik, Ar. A. Melnikov, M. Podobrii, A. Termanova, I. Kreidich, B. Nuriev, S. Iudin, C. W. Mansell, and V. M. Vinokur.  NISQ-compatible approximate quantum algorithm for unconstrained and constrained discrete optimization.  Quantum, 7:1186, 2023.
  * [71] Asel Sagingalieva, Mo Kordzanganeh, Andrii Kurkin, Artem Melnikov, Daniil Kuhmistrov, Michael Perelshtein, Alexey Melnikov, Andrea Skolik, and David Von Dollen.  Hybrid quantum ResNet for car classification and its hyperparameter optimization.  Quantum Machine Intelligence, 5:38, 2023.
  * [72] Nikita Belokonev, Artem Melnikov, Maninadh Podapaka, Karan Pinto, Markus Pflitsch, and Michael Perelshtein.  Optimization of chemical mixers design via tensor trains and quantum computing.  arXiv preprint arXiv:2304.12307, 2023.
  * [73] Dmitry Morozov, Artem Melnikov, Vishal Shete, and Michael Perelshtein.  Protein-protein docking using a tensor train black-box optimization method.  arXiv preprint arXiv:2302.03410, 2023.
  * [74] Anastasia Batsheva, Andrei Chertkov, Gleb Ryzhakov, and Ivan Oseledets.  PROTES: Probabilistic Optimization with Tensor Sampling.  arXiv preprint arXiv:2301.12162, 2023.
  * [75] Andrei Chertkov, Gleb Ryzhakov, Georgii Novikov, and Ivan Oseledets.  Optimization of functions given in the tensor train format.  arXiv preprint arXiv:2209.14808, 2022.
  * [76] Akshay Vishwanathan, Artem Melnikov, Alena Termanova, and Michael Perelshtein.  Quantum power optimization algorithm via tensor network representations.  Manuscript in Preparation, 2024.
  * [77] Alexey Naumov, Ar Melnikov, Vadim Abronin, Fedor Oxanichenko, Kamil Izmailov, Markus Pflitsch, Alexey Melnikov, and Michael Perelshtein.  Tetra-AML: automatic machine learning via tensor networks.  arXiv preprint arXiv:2303.16214, 2023.
  * [78] Georgios Laskaris, Artem A Melnikov, Michael R Perelshtein, Reuben Brasher, Thomas Baeck, and Florian Neukart.  Comparison between tensor networks and variational quantum classifier.  arXiv preprint arXiv:2311.15663, 2023.
  * [79] V Abronin, A Naumov, D Mazur, D Bystrov, K Tsarova, Ar Melnikov, I Oseledets, S Dolgov, R Brasher, and M Perelshtein.  TQCompressor: improving tensor decomposition methods in neural networks via permutations.  arXiv preprint arXiv:2401.16367, 2024.
  * [80] Trygve Helgaker, Poul Jorgensen, and Jeppe Olsen.  Molecular electronic-structure theory.  John Wiley & Sons, 2013.
  * [81] Garnet Kin-Lic Chan and Sandeep Sharma.  The density matrix renormalization group in quantum chemistry.  Annual review of physical chemistry, 62:465–481, 2011.
  * [82] Garnet Kin Chan, Anna Keselman, Naoki Nakatani, Zhendong Li, and Steven R White.  Matrix product operators, matrix product states, and ab initio density matrix renormalization group algorithms.  The Journal of chemical physics, 145(1), 2016.
  * [83] Seunghoon Lee, Joonho Lee, Huanchen Zhai, Yu Tong, Alexander M Dalzell, Ashutosh Kumar, Phillip Helms, Johnnie Gray, Zhi-Hao Cui, Wenyuan Liu, et al.  Evaluating the evidence for exponential quantum advantage in ground-state quantum chemistry.  Nature Communications, 14(1):1952, 2023.
  * [84] Justin A Reyes and E Miles Stoudenmire.  Multi-scale tensor network architecture for machine learning.  Machine Learning: Science and Technology, 2(3):035036, 2021.