Iterative tensor network transformations for element-wise evaluation of elementary and filtering functions

# Iterative tensor network transformations for element-wise evaluation of elementary and filtering functions
Source: [https://arxiv.org/html/2608.17135](https://arxiv.org/html/2608.17135)
Xiao Wang[https://orcid.org/0000-0002-3022-7260](https://orcid.org/0000-0002-3022-7260)Affiliation:Clarendon Laboratory, University of Oxford, Parks Road, Oxford OX1 3PU, United KingdomEmail:[xw970921@gmail\.com](mailto:[email protected])Thanks:These authors contributed equally to this work\.Tomohiro Hashizume[https://orcid.org/0000-0002-7154-5417](https://orcid.org/0000-0002-7154-5417)Affiliation:The Hamburg Centre for Ultrafast Imaging, Luruper Chaussee 149, Hamburg 22761, Germany\.Affiliation:Institute for Quantum Physics, University of Hamburg, Luruper Chaussee 149, Hamburg 22761, GermanyEmail:[tomohiro\.hashizume@uni\-hamburg\.de](mailto:[email protected])Thanks:These authors contributed equally to this work\.Pia Siegl[https://orcid.org/0000-0003-2249-8121](https://orcid.org/0000-0003-2249-8121)Affiliation:Institute for Quantum Physics, University of Hamburg, Luruper Chaussee 149, Hamburg 22761, GermanyAffiliation:Institute of Software Methods for Product Virtualization, German Aerospace Center \(DLR\), Nöthnitzer Straße 46b, 01187 Dresden, GermanyDieter Jaksch[https://orcid.org/0000-0002-9704-3941](https://orcid.org/0000-0002-9704-3941)Affiliation:Clarendon Laboratory, University of Oxford, Parks Road, Oxford OX1 3PU, United KingdomAffiliation:The Hamburg Centre for Ultrafast Imaging, Luruper Chaussee 149, Hamburg 22761, Germany\.Affiliation:Institute for Quantum Physics, University of Hamburg, Luruper Chaussee 149, Hamburg 22761, Germany

August 17, 2026

###### Abstract

Tensor networks are powerful formats for compressing large\-scale data\. However, their application to general data processing has been limited by the difficulty of performing nonlinear operations\. Here, we introduce iterative tensor network transformations \(ITNTs\), a general algorithmic framework for the element\-wise evaluation of elementary and nonlinear filtering functions on data encoded as tensor trains \(TTs\), a class of tensor networks\. Our approach operates entirely in the compressed domain, enabling efficient computation on exponentially large datasets while maintaining a controlled computational cost\. We demonstrate its power in two key areas: \(I\) evaluating highly nonlinear elementary and filtering functions on a 3D reactive flow field, enabling high\-fidelity reaction rate computation and region filtering, and \(II\) finding extrema in complex optimization problems, such as solving Max\-SAT instances on spaces up to2702^\{70\}configurations\. These results establish ITNT as a foundational tool that provides tensor network methods with the capability for general\-purpose data science and large\-scale optimization\.

###### Keywords:

quantum\-inspired algorithm, tensor trains, discrete nonlinear maps, data processing, max\-3SAT problem

## IIntroduction

The efficient manipulation of large\-scale datasets is a central challenge in modern science, particularly in the simulation of complex systems and the training of machine learning models\. In quantum many\-body physics, this challenge manifests itself as thecurse of dimensionality, which has been successfully addressed by tensor network algorithms\. Among these, the tensor train \(TT\) is a one\-dimensional tensor network that represents a data\-encoded array of sizeNNas a product of smaller tensors\([1](https://arxiv.org/html/2608.17135#bib.bib1);[2](https://arxiv.org/html/2608.17135#bib.bib2);[3](https://arxiv.org/html/2608.17135#bib.bib3);[4](https://arxiv.org/html/2608.17135#bib.bib4);[5](https://arxiv.org/html/2608.17135#bib.bib5)\)\. In many cases, the TT structure enables an exponential reduction of the computational and memory requirements to𝒪⁡\(log⁡N\)\\mathcal\{O\}\(\\log N\)\. This powerful data compression has enabled large\-scale simulations of quantum systems, far beyond the reach of methods that operate on the full state vector\. Owing to this compression efficiency, TT algorithms have been extended beyond their quantum origins to classical problems, including fluid dynamics\([6](https://arxiv.org/html/2608.17135#bib.bib6);[7](https://arxiv.org/html/2608.17135#bib.bib7);[8](https://arxiv.org/html/2608.17135#bib.bib8);[9](https://arxiv.org/html/2608.17135#bib.bib9);[10](https://arxiv.org/html/2608.17135#bib.bib10);[11](https://arxiv.org/html/2608.17135#bib.bib11)\), plasma physics\([12](https://arxiv.org/html/2608.17135#bib.bib12)\), and machine learning\([13](https://arxiv.org/html/2608.17135#bib.bib13)\)\. Driven by hardware advancements and algorithmic progress in performing operations, tensor network simulations are beginning to enable results across diverse physical systems far beyond the regime where simple decompression into an exact numerical vector format is possible\([14](https://arxiv.org/html/2608.17135#bib.bib14);[15](https://arxiv.org/html/2608.17135#bib.bib15);[16](https://arxiv.org/html/2608.17135#bib.bib16)\)\.

In this regime, the most critical physical quantities in frontier science and engineering are often inherently nonlinear, necessitating that the corresponding transformations be performed directly within the compressed TT representation\. This requirement, however, runs against a structural limitation of the format\. A TT supports only a small set of operations with explicit, rank\-controlled implementations: the linear operations, \(partial\) integrations that take conditional sums of the element, and, as its only nonlinear primitive, the element\-wise \(Hadamard\) product under which the bond dimensions of the factors multiply\. A generic nonlinear function applied element\-wise to a compressed array therefore admits no direct TT implementation\. One route around this limitation is to construct the transformed data by interpolation\([17](https://arxiv.org/html/2608.17135#bib.bib17);[18](https://arxiv.org/html/2608.17135#bib.bib18);[19](https://arxiv.org/html/2608.17135#bib.bib19);[20](https://arxiv.org/html/2608.17135#bib.bib20)\); for complex or steep transformations, however, interpolation fails: sharp discontinuities, emergent singularities, and rugged underlying landscapes render the data non\-interpolatable\. As Jensen’s inequality\([21](https://arxiv.org/html/2608.17135#bib.bib21);[22](https://arxiv.org/html/2608.17135#bib.bib22)\)illustrates, loss of fine structure due to compression or interpolation introduces systematic biases when highly nonlinear observables are evaluated\([23](https://arxiv.org/html/2608.17135#bib.bib23);[24](https://arxiv.org/html/2608.17135#bib.bib24)\)\. This serves as a primary bottleneck to high\-accuracy tensor network data processing in domains reliant on resolving sharp gradients, including potential energy mapping of chemical compounds\([25](https://arxiv.org/html/2608.17135#bib.bib25)\), turbulent flow simulations\([26](https://arxiv.org/html/2608.17135#bib.bib26);[9](https://arxiv.org/html/2608.17135#bib.bib9);[11](https://arxiv.org/html/2608.17135#bib.bib11)\), and statistical models encoding combinatorial optimization problems\([27](https://arxiv.org/html/2608.17135#bib.bib27)\)\. In principle, the native operations can be composed into convergent schemes, such as polynomial expansions built from repeated element\-wise products, or fixed\-point iterations for steeper targets, to realize arbitrary functions\. Whether such repeated, truncated compositions remain accurate for the nonlinear transformations relevant to industrial and scientific applications has, however, not yet been investigated\.

In this article, we address this question by introducing iterative tensor network transformations \(ITNTs\), a general framework for applying arbitrary functions to data encoded as a TT\. ITNTs leverage the iterative application of a small set of efficient TT operations: linear operations, \(partial\) integrations, and element\-wise multiplications \(Fig\.[1](https://arxiv.org/html/2608.17135#S1.F1)\(a\)\)\. This enables the computation of functions ranging from elementary functions to highly nonlinear operations, such as filtering, to be applied overNNpoints in parallel, performed entirely within the compressed domain in a highly controlled manner\. To demonstrate the power of ITNTs, we benchmark their application to large\-scale continuous fields and discrete energy landscapes\. For continuous fields, we compute the reaction rate of a 3D reactive flow field by applying ITNTs in TT form, giving rise to an order\-of\-magnitude higher fidelity in comparison to the interpolation method \(Fig\.[1](https://arxiv.org/html/2608.17135#S1.F1)\(b\)\)\. For discrete optimization, we find near\-optimal configurations for the NP\-hard Max\-SAT problem with a rugged landscape by filtering out the irrelevant bulk structure via nonlinear operations \(Fig\.[1](https://arxiv.org/html/2608.17135#S1.F1)\(c\)\)\. Notably, we identify a solution matching those found by state\-of\-the\-art heuristic algorithms\([28](https://arxiv.org/html/2608.17135#bib.bib28)\)from a configuration space ofN=270≈1021N=2^\{70\}\\approx 10^\{21\}states\. Furthermore, we show that controlled truncation within the ITNT framework enables explicit verification of solution validity, establishing a concrete link between tensor network computational theory and fundamental computational complexity theory\. By enabling the element\-wise evaluation of arbitrary functions in a compressed format, our framework opens an expansive class of problems in large\-scale scientific computing, data science, and combinatorial optimization to tensor network methods\.

![Refer to caption](https://arxiv.org/html/2608.17135v1/figures/overview.png)Figure 1:Overview of the results\.\(a\) Schematic of the Iterative tensor network transformations \(ITNTs\) methodology\. The framework enables the evaluation of highly complex nonlinear transformations,𝒯\\mathcal\{T\}, directly on TT\-encoded dataΨ\\Psidefined on array indices𝒋∈\{0,1,…,N−1=2n−1\}\\bm\{j\}\\in\\\{0,1,\\ldots,N\-1=2^\{n\}\-1\\\}\. This is achieved by leveraging2n2^\{n\}parallel evaluations of the iteratively converging Newton\-Raphson method \(bottom\-right panels\)\. In these panels, the red curve represents the nonlinear function being iteratively solved, while the blue and green lines illustrate the successive tangent\-linear approximations during the iterative process, starting from the blue dot\. This convergence allows the transformation to be applied while maintaining the data in a compressed format\. \(b\) Validation of the framework on a simulated methane/air jet flame\([26](https://arxiv.org/html/2608.17135#bib.bib26)\)\. The ITNT method is used to compute the highly nonlinear Arrhenius reaction rate \(exp\(−Ec/T\)\\exp\(\-E\_\{c\}/T\)\)\. The plot shows the iso\-rate surface at 10% of the maximum reaction rate, computed directly from a compressed temperature fieldTT\[K\] using a normalized reaction energy ofEc=17865\.2E\_\{c\}=17865\.2\[K\]\. The comparison uses a TT bond dimension ofχ=200\\chi=200, requiring only 0\.4The results demonstrate high\-fidelity reconstruction, with ITNT achieving enclosed volume relative error approximately 100 times lower than that of TCI, representing a significant improvement in capturing the flame’s topological features\. \(c\) Schematic of combinatorial optimization utilizing ITNT\. The framework uses iterative nonlinear operations, specifically self\-multiplication, to amplify and sharpen dominant configurations, such as the global maximum \(EoptE\_\{\\mathrm\{opt\}\}\) and sub\-optimal peaks \(EsubE\_\{\\mathrm\{sub\}\}\), within the problem\-encoded energy landscape\. Simultaneously, TT truncation \(depicted by the lightning bolt\) prunes away the rugged bulk spanned by the low\-energy configurations by discarding the subspace corresponding to smaller singular values\. This allows the algorithm to isolate optimal solutions while maintaining the data in a compressed TT format\.
## IIResults

LetΨ\\Psidenote an array of lengthNN, andΨ\[j\]\\Psi^\{\[j\]\}\(j∈\{0,⋯,N−1\}j\\in\\\{0,\\cdots,N\-1\\\}\) be thejjthelement of this array\. Here, square brackets are used for superscripts to distinguish them from exponents\. The TT representation compresses the element\-wise numerical representation ofΨ\\Psiinto a product ofnnlocal rank\-33tensorsΨ⁡\[q\]\\Psi\[q\], such thatΨ\[j\]≈∑𝜶Ψ\[0\]α0,α1\[j0\]Ψ\[1\]α1,α2\[j1\]⋯Ψ\[n−1\]αn−1,αn\[jn−1\]\\Psi^\{\[j\]\}\\approx\\sum\_\{\\bm\{\\alpha\}\}\\Psi\[0\]^\{\[j\_\{0\}\]\}\_\{\\alpha\_\{0\},\\alpha\_\{1\}\}\\Psi\[1\]^\{\[j\_\{1\}\]\}\_\{\\alpha\_\{1\},\\alpha\_\{2\}\}\\cdots\\Psi\[n\-1\]^\{\[j\_\{n\-1\}\]\}\_\{\\alpha\_\{n\-1\},\\alpha\_\{n\}\}\. In this formulation,jq∈\{0,1\}j\_\{q\}\\in\\\{0,1\\\}corresponds to theqqthbit of the binary representation of indexjj, denoted by the sequence𝒋∈\{j0,j1,⋯,jn−1\}\\bm\{j\}\\in\\\{j\_\{0\},j\_\{1\},\\cdots,j\_\{n\-1\}\\\}, andαr\\alpha\_\{r\}corresponds to the bond index of therrthbond\. The summation runs over all the bond indices𝜶=\{α0,α1,⋯,αn\}\\bm\{\\alpha\}=\\\{\\alpha\_\{0\},\\alpha\_\{1\},\\cdots,\\alpha\_\{n\}\\\}, each of which may have a size up toχ\\chi, the maximum bond dimension that controls the level of compression \(Methods[V\.1](https://arxiv.org/html/2608.17135#S5.SS1)\)\. While a TT can exactly represent any arrayΨ\\Psigiven an exponentially large maximum bond dimensionχ=2n/2\\chi=2^\{n/2\}, in practice,χ\\chiis kept at a much smaller constant, e\.g\.,χ=200\\chi=200for Fig\.[1](https://arxiv.org/html/2608.17135#S1.F1)\(b\) and provides an approximation ofΨ\[j\]\\Psi^\{\[j\]\}\. This value allows us to balance the minimization of compression loss against the reduced memory requirements𝒪⁡\(n​χ2\)\\mathcal\{O\}\(n\\chi^\{2\}\), and the computational scaling𝒪⁡\(poly⁡\(n\)​poly​\(χ\)\)\\mathcal\{O\}\(\\mathrm\{poly\}\(n\)\\mathrm\{poly\}\(\\chi\)\)\.

To proceed, we introduce the elementary TT operations that constitute the ITNT framework, as illustrated in Fig\.[1](https://arxiv.org/html/2608.17135#S1.F1)\(a\): scalar multiplication, variational addition, element\-wise products, and partial integration\. Specifically, scalar multiplication is defined asκ​Ψ\\kappa\\Psi, which is implemented by applying the scalarκ\\kappato a single local tensor of the TT\. Addition of two TTsΨ1\+Ψ2\\Psi\_\{1\}\+\\Psi\_\{2\}is a well\-established, stable routine and requires a computational cost scaling as𝒪⁡\(n​χ3\)\\mathcal\{O\}\(n\\chi^\{3\}\)\([4](https://arxiv.org/html/2608.17135#bib.bib4)\)\. The element\-wise \(Hadamard\) product,Ψ1⊙Ψ2\\Psi\_\{1\}\\odot\\Psi\_\{2\}, is realized with a zip\-up algorithm using copy tensors\([29](https://arxiv.org/html/2608.17135#bib.bib29)\), resulting in a scaling of𝒪⁡\(n​χ4\)\\mathcal\{O\}\(n\\chi^\{4\}\)\([30](https://arxiv.org/html/2608.17135#bib.bib30)\)or allowing inflation in the intermediate bond dimension,𝒪⁡\(n​χ3\)\\mathcal\{O\}\(n\\chi^\{3\}\)\([20](https://arxiv.org/html/2608.17135#bib.bib20)\)\. Finally, the partial integration operatorℐ\(Ψ,q\)=∑𝒋\\\{jq\}∑𝜶Ψ\[0\]α0,α1\[j0\]Ψ\[1\]α1,α2\[j1\]⋯Ψ\[n−1\]αn−1,αn\[jn−1\]\\mathcal\{I\}\(\\Psi,q\)=\\sum\_\{\\bm\{j\}\\backslash\\\{j\_\{q\}\\\}\}\\sum\_\{\\bm\{\\alpha\}\}\\Psi\[0\]^\{\[j\_\{0\}\]\}\_\{\\alpha\_\{0\},\\alpha\_\{1\}\}\\Psi\[1\]^\{\[j\_\{1\}\]\}\_\{\\alpha\_\{1\},\\alpha\_\{2\}\}\\cdots\\Psi\[n\-1\]^\{\[j\_\{n\-1\}\]\}\_\{\\alpha\_\{n\-1\},\\alpha\_\{n\}\}constructs a two\-element array by summing all elements for whichjq=0j\_\{q\}=0andjq=1j\_\{q\}=1, respectively \(Methods[V\.2](https://arxiv.org/html/2608.17135#S5.SS2)\)\.

These primitive operations immediately enable the element\-wise evaluation of integer polynomials and exponentiation directly from their basic definitions\. Consequently, elementary functions such as sine, cosine, hyperbolic sine, and hyperbolic cosine, which are linear combinations of \(complex\) exponentiation, can also be implemented within this framework\. To implement other elementary functions𝒯\\mathcal\{T\}\([31](https://arxiv.org/html/2608.17135#bib.bib31)\), particularly those whose inverses are readily computable from the aforementioned primitives, we employ the Newton\-Raphson method \(Fig\.[1](https://arxiv.org/html/2608.17135#S1.F1)\(a\), bottom right\)\. This involves computing the rootΨ′\\Psi^\{\\prime\}that satisfies𝒯−1​\(Ψ′\)−Ψ=0\\mathcal\{T\}^\{\-1\}\(\\Psi^\{\\prime\}\)\-\\Psi=0, where𝒯−1​\(Ψ\)\\mathcal\{T\}^\{\-1\}\(\\Psi\)is the inverse transformation satisfying𝒯−1​\(𝒯​\(Ψ\)\)=Ψ\\mathcal\{T\}^\{\-1\}\(\\mathcal\{T\}\(\\Psi\)\)=\\Psi\. Thus, the elementary functions defined by Bourchteinet al\.\([31](https://arxiv.org/html/2608.17135#bib.bib31)\)can be fully implemented within the TT framework \(for details, see Methods[V\.3](https://arxiv.org/html/2608.17135#S5.SS3)\)\. To demonstrate the robustness and accuracy of this approach, in addition to the methane/air jet flame in Fig\.[1](https://arxiv.org/html/2608.17135#S1.F1)\(b\), we apply this method to solve the Kidder equation\([32](https://arxiv.org/html/2608.17135#bib.bib32);[33](https://arxiv.org/html/2608.17135#bib.bib33);[34](https://arxiv.org/html/2608.17135#bib.bib34)\)in SM[I](https://arxiv.org/html/2608.17135#S1a)\. The Kidder equation is a strongly nonlinear partial differential equation that models flow through porous media with a square\-root dependency on the field\. Our numerical results show excellent agreement with state\-of\-the\-art analytical solutions on the initial slope of the density field, matching them to six decimal places\.

Elementary iterative methods fail for highly nonlinear, non\-invertible functions, such as the sign function or extremum finding, both of which are often required in data processing\. Here, tools originally developed for homomorphic cryptography\([35](https://arxiv.org/html/2608.17135#bib.bib35)\)are required, as we will discuss below\.

#### Computing the sign function\.

Our ITNT sign function algorithm is based on the homomorphic comparison of two encrypted numbers\([36](https://arxiv.org/html/2608.17135#bib.bib36)\), which implements comparisons using only addition, scalar multiplication, and element\-wise multiplication\. The TTΨ\\Psiis compared to a zero field∅\\varnothing\(∅\[j\]=0\\varnothing^\{\[j\]\}=0for alljj\), yielding an element\-wise evaluation of the sign function with the outcomeSign​\{Ψ\}\[j\]=Sign⁡\{Ψ\[j\]\}\\mathrm\{Sign\}\\\{\\Psi\\\}^\{\[j\]\}=\\mathrm\{Sign\}\\\{\\Psi^\{\[j\]\}\\\}, whereSign​\{x\}=1\\mathrm\{Sign\}\\\{x\\\}=1\(x\>0x\>0\),Sign​\{x\}=0\\mathrm\{Sign\}\\\{x\\\}=0forx=0x=0, andSign​\{x\}=−1\\mathrm\{Sign\}\\\{x\\\}=\-1\(x<0x<0\)\. The implementation is efficiently achieved via the iterationR\(k\+1\)=F⁡\(R\(k\)\)R^\{\(k\+1\)\}=F\\left\(R^\{\(k\)\}\\right\)withR\(0\)=ΨR^\{\(0\)\}=\\Psiwhere

F⁡\(R\(k\)\)=−12​\(\(R\(k\)\)⊙3−3​R\(k\)\)\\displaystyle F\(R^\{\(k\)\}\)=\-\\frac\{1\}\{2\}\\left\(\(R^\{\(k\)\}\)^\{\\odot 3\}\-3R^\{\(k\)\}\\right\)\(1\)andΨ⊙θ\\Psi^\{\\odot\\theta\}denotes theθ\\thetathelement\-wise power\. Convergence of this iteration toSign​\{Ψ\}\\mathrm\{Sign\}\\\{\\Psi\\\}requires\|Ψ\[j\]\|<3\|\\Psi^\{\[j\]\}\|<\\sqrt\{3\}\(∀j\\forall j\)\. Therefore,Ψ\\Psimight first need to be rescaled to ensure that the iterations remain within the radius of convergence\. PreconditioningΨ\\Psifurther may speed up convergence \(for details, see Methods[subsection V\.4](https://arxiv.org/html/2608.17135#S5.SS4)\)\. A rough estimate of the normalization constant may be obtained by performing a few iterations of the extremum\-finding ITNT discussed next\.

#### Extremum\-finding ITNT\.

We utilize partial integration to determine, one bit at a time, the configuration that contributes most to the integral\. For a fieldΨ\\PsiwithΨ\[j\]≥0\\Psi^\{\[j\]\}\\geq 0, the toal integral obeys∑jq∈\{0,1\}ℐ​\(Ψ,q\)\[jq\]=∑jΨ\[j\]\\sum\_\{j\_\{q\}\\in\\\{0,1\\\}\}\\mathcal\{I\}\(\\Psi,q\)^\{\[j\_\{q\}\]\}=\\sum\_\{j\}\\Psi^\{\[j\]\}; if this sum is dominated by the maximum elementmax⁡\{Ψ\[j\]\}\\max\\\{\\Psi^\{\[j\]\}\\\}, thenargmax​\{ℐ\[jq\]​\(Ψ,q\)\}\\mathrm\{argmax\}\\\{\\mathcal\{I\}^\{\[j\_\{q\}\]\}\(\\Psi,q\)\\\}coincides with theqqthcomponent of the global maximum indexJ=argmax​\{Ψ\[j\]\}J=\\mathrm\{argmax\}\\\{\\Psi^\{\[j\]\}\\\}\. We therefore define the binary fixation \(BF\) as sequentially fixing theqqthindex of a TT tomq=argmax⁡\{ℐ\[jq\]​\(Ψ,q\)\}m\_\{q\}=\\mathrm\{argmax\}\\\{\\mathcal\{I\}^\{\[j\_\{q\}\]\}\(\\Psi,q\)\\\}, starting from the00thtensor\. BF recovers the full index of the maximum𝒎=\(m0,m1,⋯,mn−1\)=𝑱\\bm\{m\}=\(m\_\{0\},m\_\{1\},\\cdots,m\_\{n\-1\}\)=\\bm\{J\}afternniterations\.

In practice, a raw field rarely holds such a dominance condition\. Therefore, as a preprocessing, we propose to raise the field components to a power2P2^\{P\}via self\-multiplications by performingPPsteps of self\-multiplications\. This concentrates the field onto its largest elements while preserving the location of the global maximum, and for sufficiently largePP, the dominance condition holds for the exactΨ⊙2P\\Psi^\{\\odot 2^\{P\}\}\. However, this self\-multiplication causes the bond dimension to grow asχ2P\\chi^\{2^\{P\}\}without bond dimension truncation\. To keep the bond dimension of the output TT manageable, we use singular value decomposition \(SVD\) for truncation\([37](https://arxiv.org/html/2608.17135#bib.bib37)\)\. SVD, by design, minimizes theℓ2\\ell\_\{2\}error, thus the truncation flattens and removes the rugged bulk that requires a large bond dimension for its exact representation, while the landscape near maximal values is preserved due to its large contributions to the globalℓ2\\ell\_\{2\}norm\.

Theℓ2\\ell\_\{2\}error, defined asε⁡\(Ψ,ΨSVD\)=‖Ψ−ΨSVD‖2\\varepsilon\(\\Psi,\\Psi\_\{\\mathrm\{SVD\}\}\)=\\\|\\Psi\-\\Psi\_\{\\mathrm\{SVD\}\}\\\|\_\{2\}, whereΨSVD\\Psi\_\{\\mathrm\{SVD\}\}is theΨ\\PsiTT obtained after SVD truncation, bounds theℓ∞\\ell\_\{\\infty\}error \(ϵ∞=‖Ψ−ΨSVD‖∞\\epsilon\_\{\\infty\}=\\\|\\Psi\-\\Psi\_\{\\mathrm\{SVD\}\}\\\|\_\{\\infty\}\), the maximum perturbation added to the elements, from the above\. Thus, the algorithm is guaranteed to converge to a bond dimension of, at worst,χf\\chi\_\{f\}corresponding to the number of indices that give rise tomax⁡\{Ψ\}\\max\\\{\\Psi\\\}, and BF finds one of these indices as long as the accumulatedℓ2\\ell\_\{2\}error per step at theppth self\-multiplication step satisfies

ε⁡\(Ψ⊙2p,\(Ψ⊙2p\)SVD\)≤\(max⁡\{\(Ψ⊙2p\)\[j\]\}−maxsub⁡\{\(Ψ⊙2p\)\[j\]\}\)/2\\displaystyle\\varepsilon\(\\Psi^\{\\odot\{2^\{p\}\}\},\(\\Psi^\{\\odot\{2^\{p\}\}\}\)\_\{\\mathrm\{SVD\}\}\)\\leq\(\\max\\\{\(\\Psi^\{\\odot 2^\{p\}\}\)^\{\[j\]\}\\\}\-\\max\_\{\\mathrm\{sub\}\}\\\{\(\\Psi^\{\\odot 2^\{p\}\}\)^\{\[j\]\}\\\}\)/\\sqrt\{2\}\(2\)for the gap betweenmax⁡\{\(Ψ⊙2p\)\[j\]\}\\max\\\{\(\\Psi^\{\\odot 2^\{p\}\}\)^\{\[j\]\}\\\}andmaxsub⁡\{\(Ψ⊙2p\)\[j\]\}\\max\_\{\\mathrm\{sub\}\}\\\{\(\\Psi^\{\\odot 2^\{p\}\}\)^\{\[j\]\}\\\}the largest element strictly below the maximum \(cf\. SM[II](https://arxiv.org/html/2608.17135#S2a)for details\)\.

When the required number of self\-multiplications is intractably large, theℓ2\\ell\_\{2\}error also becomes large due to the limited bond dimension\. In such a case, the maximum\-finding no longer becomes exact, and one obtains a heuristic, approximate solution, which represents the best estimate ofJJattainable within the available computational resources\. The key is to perform BF in the randomized order, referred to as random binary fixation \(RBF\), and remove the identified peak from the landscape \(deflation\)\. The RBF procedure applies BF iterations to randomly selected sequences of local tensors\. The subsequent deflation step,Ψ′​→defl​Ψ′′\\Psi^\{\\prime\}\\overset\{\\text\{defl\}\}\{\\rightarrow\}\\Psi^\{\\prime\\prime\}, removes the identified configuration fromΨ′\\Psi^\{\\prime\}such thatΨ′′\[j\]=Ψ′\[j\]−δj,rΨ′\[r\]\\Psi^\{\\prime\\prime\[j\]\}=\\Psi^\{\\prime\[j\]\}\-\\delta\_\{j,r\}\\Psi^\{\\prime\[r\]\}\(see Methods[subsection V\.5](https://arxiv.org/html/2608.17135#S5.SS5)\), whererris the index found in a given round of RBF\. After the assigned amount of RBF resource is exhausted, the configuration found that yields the largest value ofΨ\\Psiis accepted as the approximate solution to the global maximum\.

For the rest of this article, we focus on studying the numerical performance of nonlinear functions applied to the most difficult structured data arising from both continuous and discrete problems\. As nonlinear elementary functions, filtering, and extremum\-finding represent the most demanding nonlinearities for TT\-based numerical methods, they serve as a direct benchmark for the overall robustness and precision of the ITNT framework\.

#### Application I: Data processing of a 3D reactive flow\.

We first demonstrate the utility of the ITNT framework for large\-scale data processing by applying it to high\-resolution grid data\. The TT algorithms for extremum\-finding and sign\-function evaluation can be directly applied to field\-processing tasks, as illustrated in Fig\.[1](https://arxiv.org/html/2608.17135#S1.F1)\(b\)\. Here, we use them to post\-process the temperature field of a precomputed methane/air jet flame\([26](https://arxiv.org/html/2608.17135#bib.bib26)\)\. This demonstrates how ITNT implements nonlinear functions by computing the Arrhenius reaction rate as an example\.

![Refer to caption](https://arxiv.org/html/2608.17135v1/figures/usecase_jet.png)Figure 2:Nonlinear post processing with ITNTs\.\(a\) Convergence comparison between ITNT and TCI for the Arrhenius equation \(Eq\. \([3](https://arxiv.org/html/2608.17135#S2.E3)\)\) computation\. While TCI suffers from instability aroundχ=200\\chi=200, ITNT demonstrates a drop in error in the same region, consistently converging towards the exact result with increasingχ\\chi\. \(b\) Relativeℓ2\\ell\_\{2\}error of the1/T1/Tfield computed exactly \(Ψexact\\Psi\_\{\\mathrm\{exact\}\}\) from the original field and by using the ITNT algorithm \(ΨTT\\Psi\_\{\\mathrm\{TT\}\}\) forχ=50\\chi=50,100100,200200,300300, and500500\(light to dark, left panel\), and the same error at the 10th iteration for different bond dimensions\. In the right panel, the blue dashed line indicates1/χ1/\\chidecay as a guide to the eye\. \(c\) The position of the maximum temperature is indicated by a red dot and the orange region represents the isothermal surface at9090% of the maximum temperature \(0\.9​max⁡\{T\}0\.9\\max\\\{T\\\}\)\. \(d\) Relative error of the maximum temperature,1−max⁡\{TTT\}/max⁡\{Texact\}1\-\\max\\\{T\_\{\\mathrm\{TT\}\}\\\}/\\max\\\{T\_\{\\mathrm\{exact\}\}\\\}, found using ITNT forχ=50\\chi=50,100100,200200,300300, and500500\(light to dark\)\. \(e\) Relative error of the volume of the region whereT\>0\.9​max⁡\{T\}T\>0\.9\\max\\\{T\\\}forχ=50\\chi=50,100100,200200,300300, and500500\(light to dark\)\. Here, the value of0\.9​max⁡\{T\}0\.9\\max\\\{T\\\}is computed with the maximum\-finding ITNT and the volume is determined by constructing the filter with sign\-transformation ITNT\.Formally, the Arrhenius reaction rate is defined as\([38](https://arxiv.org/html/2608.17135#bib.bib38)\),

Ar=A0exp\(−Ec/T\),\\displaystyle A\_\{r\}=A\_\{0\}\\exp\(\-E\_\{c\}/T\),\(3\)whereTTis the temperature andEc=17865\.2E\_\{c\}=17865\.2\[K\] is the activation energy of methane combustion\. For simplicity, we set the pre\-exponential factorA0=1A\_\{0\}=1\. This rate quantifies the chemical kinetics within the flame, which are critical factors in engine design\.

To perform the computation, we encode the temperature field into a TT representation with a maximum bond dimensionχ\\chi\. Before encoding, we trim the outer grid points so that each grid dimension is discretized to the nearest power of22\. We then apply the elementary\-function ITNT \(see Methods[V\.3](https://arxiv.org/html/2608.17135#S5.SS3)\) in two stages: the first involves computing the reciprocal of the temperature field at each point, and the second involves exponentiating the result from the previous step\. Throughout these operations, we cap the bond dimension of all intermediate tensors atχ\\chi, where we obtain the TT representation ofArA\_\{r\}with a final bond dimension ofχ\\chi\. Alternatively, for comparison, a Tensor Cross Interpolation \(TCI\) approach is employed via thexfaclibrary\([18](https://arxiv.org/html/2608.17135#bib.bib18)\)to map the nonlinear function directly\. This TCI procedure utilizes an iterative pivot\-based search, starting from the maximum temperature location to ensure a stable initial estimate, refining the tensor cores until the pivot error falls below a tolerance ofϵ=10−12\\epsilon=10^\{\-12\}or the maximum bond dimension ofχ\\chiis reached\.

As shown in Fig\.[1](https://arxiv.org/html/2608.17135#S1.F1)\(b\), ITNT provides a good approximation to the true reaction rate features already atχ=200\\chi=200, which corresponds to 0\.4% of the original size\. In contrast, TCI fails to reproduce most of the features and creates artifacts\. This difference in behavior is further observed in the convergence of the relativeℓ2\\ell\_\{2\}error

ε~​\(Ψexact,ΨTT\)=‖Ψexact−ΨTT‖2/‖Ψexact‖2\\displaystyle\\tilde\{\\varepsilon\}\(\\Psi\_\{\\mathrm\{exact\}\},\\Psi\_\{\\mathrm\{TT\}\}\)=\\\|\\Psi\_\{\\mathrm\{exact\}\}\-\\Psi\_\{\\mathrm\{TT\}\}\\\|\_\{2\}/\\\|\\Psi\_\{\\mathrm\{exact\}\}\\\|\_\{2\}\(4\)between the exact representationΨexact\\Psi\_\{\\mathrm\{exact\}\}and the TT representationΨTT\\Psi\_\{\\mathrm\{TT\}\}\. As shown in Fig\.[2](https://arxiv.org/html/2608.17135#S2.F2)\(a\), TCI exhibits a region of instability around bond dimensionχ=200\\chi=200, while ITNT shows a steady decrease in error without encountering this instability, demonstrating that ITNT succeeds in regimes where the interpolation method fails\.

This strong numerical stability and convergence with the bond dimensionχ\\chiare further supported by the rapid quadratic convergence enabled by the Newton\-Raphson method, as shown in Fig\.[2](https://arxiv.org/html/2608.17135#S2.F2)\(b\) for the first step, which computes the reciprocal field \(1/T1/Tfield\)\. Here, even atχ=200\\chi=200and after only 5 iterations, ITNT converges rapidly and approximates the true1/T1/Tfield to within a relative error below 2% \(left panel, red\)\. Furthermore, we observe a1/χ1/\\chidecay in the error with respect to the bond dimension \(right panel\)\.

In Fig\.[2](https://arxiv.org/html/2608.17135#S2.F2)\(c\)\-\(e\), we further perform a filtering operation on the temperature field and identify the high\-temperature region shown in \(c\)\. To this end, we identify the maximum temperature of the field via the maximum\-finding ITNT and compute the volume of the region with a temperature greater than 90% of the identified maximum\. Under the naive assumption of underlying smoothness and continuity of the encoded field, Bayesian optimization\([39](https://arxiv.org/html/2608.17135#bib.bib39)\)is known to perform well\. However, as shown in \(c\), due to chemical reactions and the resulting draft in the temperature field the maximum temperature occurs in a small, isolated region away from the main gas stream and its high\-temperature surroundings, a typical feature of combustion flows\([40](https://arxiv.org/html/2608.17135#bib.bib40)\)\. This feature makes standard sampling or learning algorithms inefficient for detecting the maximum\.

To find the peak temperature withouta prioriknowledge of the field’s exact representation, we perform the maximum\-finding ITNT, executing five rounds of self\-multiplication with a controlled bond dimensionχ\\chi\. The subsequent deflation steps \(Methods Sec\.[V\.5](https://arxiv.org/html/2608.17135#S5.SS5)\) immediately find a candidate value within 1% of the true maximum \(Fig\.[2](https://arxiv.org/html/2608.17135#S2.F2)\(d\)\)\. Then, using the sign\-transformation ITNT, we compute the volume of the region exceeding 90% of the found maximum temperature for a givenχ\\chi\. To evaluate accuracy, we plot the relative error of the volume against the true volume in Fig\.[2](https://arxiv.org/html/2608.17135#S2.F2)\(e\)\. Even for a bond dimension as small asχ=200\\chi=200, a relative error of less than 1% is achieved after 10 ITNT rounds\. The continued convergence with increasingχ\\chidemonstrates that this method can perform nonlinear data processing with high accuracy on complex data structures encoded in the TT format \(see SM[III](https://arxiv.org/html/2608.17135#S3a)for further analysis on ITNT transformations acting on continuous data\)\.

#### Application II: Solving NP\-hard optimization problems

Beyond continuous fields, we study combinatorial optimization problems as a use case for extremum\-finding on TT\-encoded energy landscapes of discrete configurations\. We focus on the class of combinatorial optimization problems that can be mapped to the energy minimization of a classical Ising model with long\-range interactions\. Specifically, in the Max\-kk\-SAT problem withnnvariables andMMclauses, we want to find a bit string𝝈=\(σ0,σ1,…,σn−1\)\\bm\{\\sigma\}=\(\\sigma\_\{0\},\\sigma\_\{1\},\\ldots,\\sigma\_\{n\-1\}\)with a corresponding indexσ\\sigmathat maximizes the followingkk\-body Ising\-type cost function,

E\[σ\]=M−∑𝒞=1M∏m=1klm𝒞,E^\{\[\\sigma\]\}=M\-\\sum\_\{\\mathcal\{C\}=1\}^\{M\}\\prod\_\{m=1\}^\{k\}l\_\{m\}^\{\\mathcal\{C\}\},\(5\)wherelm𝒞∈\{σj\(𝒞,m\),1−σj\(𝒞,m\)\}l\_\{m\}^\{\\mathcal\{C\}\}\\in\\\{\\sigma\_\{j\_\{\(\\mathcal\{C\},m\)\}\},1\-\\sigma\_\{j\_\{\(\\mathcal\{C\},m\)\}\}\\\}is chosen based on whether the corresponding literal is negated\. Here,𝒞∈\{1,…,M\}\\mathcal\{C\}\\in\\\{1,\\ldots,M\\\}is the index of the logical clause, andj\(𝒞,m\)∈\[0,n−1\]j\_\{\(\\mathcal\{C\},m\)\}\\in\[0,n\-1\]denotes the index of the bit appearing in themmthterm of𝒞\\mathcal\{C\}\. Concretely,lm𝒞l\_\{m\}^\{\\mathcal\{C\}\}is defined as a*literal\-falsity indicator*: we setlm𝒞=1−σj\(𝒞,m\)l\_\{m\}^\{\\mathcal\{C\}\}=1\-\\sigma\_\{j\_\{\(\\mathcal\{C\},m\)\}\}for an unnegated literal, and setlm𝒞=σj\(𝒞,m\)l\_\{m\}^\{\\mathcal\{C\}\}=\\sigma\_\{j\_\{\(\\mathcal\{C\},m\)\}\}for a negated literal\. With this convention,∏m=1klm𝒞=1\\prod\_\{m=1\}^\{k\}l\_\{m\}^\{\\mathcal\{C\}\}=1if and only if clause𝒞\\mathcal\{C\}is violated, meaning thatE\[σ\]E^\{\[\\sigma\]\}counts the number of satisfied clauses\.

BecauseE\[σ\]E^\{\[\\sigma\]\}counts the number of satisfied clauses, it effectively defines the energy landscape of an Ising model with variablesσj\\sigma\_\{j\}\. Consequently, the Max\-kk\-SAT problem translates into finding the global maximum of this landscape\. We create a TT that encodes the energy landscape of the Max\-33\-SAT problem, which is known to be an NP\-hard problem\([41](https://arxiv.org/html/2608.17135#bib.bib41)\), utilizing the Ising\-to\-TT encoding \(Methods[subsection V\.6](https://arxiv.org/html/2608.17135#S5.SS6)\)\. We select an instance withn=70n=70variables andM=700M=700clauses from the “Max\-SAT 2016 competition” dataset\([42](https://arxiv.org/html/2608.17135#bib.bib42)\), specifically, instances3v70c700\-1\. While an array representation of the energy landscape has2702^\{70\}elements, the mapped TT consists of7070local tensors and has a bond dimension ofχ=701\\chi=701\.

![Refer to caption](https://arxiv.org/html/2608.17135v1/figures/usecase_maxsat.png)Figure 3:Solving Max\-3\-SAT problems with ITNTs\.\(a\) Truncation error at each bond in the zip\-up self\-multiplication algorithm\. In order to solve 70 variables 700 clauses Max\-3SAT problem, 10 self\-multiplications are performed, and for each self\-multiplication is computed with the zip\-up algorithm with on\-site truncation, while fixingχ=700\\chi=700\. Each color corresponds to a different self\-multiplication iteration \(iterations 1 to 10, light to dark\) \(b\) Deflation result for the max 3\-SAT problems3v70c700\-1\. Three configurations with an energy ofEEwithE=679E=679are found after thei=1655i=1655th, thei=6442i=6442nd, and thei=30186i=30186thextremum\-finding ITNT steps\. \(c\) Scaling of the minimum bond dimension for successful optimization and solution certification\. Here, we generate 100 randomnn\-variable,M=α​nM=\\alpha n\-clause Max\-3SAT problems forn=10n=10to2020, at clause densitiesα=8\\alpha=8andα=3\\alpha=3\. Orange dashed lines show the minimum bond dimension required for the maximum search to return the correct optimum afterP=⌈log2⁡\(n​M\)⌉\+1P=\\left\\lceil\\log\_\{2\}\(nM\)\\right\\rceil\+1self\-multiplication steps\. The blue line shows the minimum bond dimension required for the truncation budget in Method[V\.7](https://arxiv.org/html/2608.17135#S5.SS7)to be obeyed afterPPsteps\. The required bond dimension grows withnn, saturating the bond dimension∼2n/2\\sim 2^\{n/2\}\(black dotted line\)\.To find the maximum number of satisfied clauses and their truth values, the landscape TT is first self\-multiplied ten times with the bond dimension fixed atχ=701\\chi=701\. Subsequently, the extremum\-locating ITNT \(RBF iterations\) is performed10510^\{5\}times\. During these RBF iterations, the maximum bond dimension is set toχmax=1000\\chi\_\{\\max\}=1000; once the TT exceeds this threshold, it triggers a fallback procedure, truncating the bond dimension back to900900\. In Fig\.[3](https://arxiv.org/html/2608.17135#S2.F3)\(a\), we plot the truncation error\. Starting from the first iteration of self\-multiplication \(lightest\), aggressive truncation of the bond dimension causes the truncation errorϵtr=∑χ<iλi2\\epsilon\_\{\\mathrm\{tr\}\}=\\sum\_\{\\chi<i\}\\lambda\_\{i\}^\{2\}, where the singular valuesλi\\lambda\_\{i\}are sorted in descending order, to reach103210^\{32\}\. Meanwhile, in the later stages of the self\-multiplication iterations, the truncation error gradually decreases, showing that the encoded energy landscape converges towards a few dominant configurations \(compared to the full2702^\{70\}\)\.

Through these iterative steps, we identify three configurations with the extremal valueE=679E=679\(Fig\.[3](https://arxiv.org/html/2608.17135#S2.F3)\(b\)\)\. These match the best\-known solution for this instance and are separated by Hamming distances of11,1212, and1313:

`1001001000110000011001111011010101001000000111010010000000110110010101``1001001000110000011001111011010101001010000111010010000000110110010101``1001001110110000010001111011000101011000011001110010000000111100010101`

Importantly, ITNT identifies these solutions while probing only a tiny fraction of the full configuration space\. Instead of exhaustively evaluating all270≈10212^\{70\}\\approx 10^\{21\}elements, the algorithm visits only𝒪⁡\(104\)\\mathcal\{O\}\(10^\{4\}\)peaks in the landscape\. The slightly suboptimal extremum \(the next\-best solution withE=678E=678\) is identified early in the process after 11 steps, providing a high\-quality, approximate solution\. This suggests an early\-stopping mode for the algorithm; if a near\-optimal solution is sufficient, the number of deflation steps can be significantly reduced\. This example illustrates that ITNT has the capability to reach state\-of\-the\-art solutions and is effective in solving discrete NP\-hard optimization problems directly within the TT framework\.

Furthermore, a fundamental distinction between the ITNT approach and traditional heuristic methods is its inherently global search mechanism\. Unlike traditional solvers, which are sequential and stochastic, the algebraic framework aggregates information across the entire ensemble of configurations, sidestepping the limitations of incremental, local spin flips\. By comparing partial sums over an exponential number of configurations in a single operation, the RBF procedure identifies global extrema with guaranteed convergence in exact arithmetic, given a sufficient number of self\-multiplication steps\. This global comparison is performed efficiently within the TT manifold, offering an exponential advantage over brute\-force numerics\.

Beyond merely finding a solution, ITNT offers a method to verify whether an identified solution is the true maximum\. Under truncation, the convergence guarantee above survives only if the error injected at each step remains within a budget set by the margin separating the optimum from the rest of the landscape\. The certification of the global maximum is done by monitoring the truncation error injected at each self\-multiplication step by assuring the dominance of the true maximum afterPPself\-multiplication steps \(Methods[V\.7](https://arxiv.org/html/2608.17135#S5.SS7)\)\. However, the unconditional success of such a verification method would imply that NP\-hard problems can be solved in polynomial time, which would contradict the P≠\\neqNP conjecture\([41](https://arxiv.org/html/2608.17135#bib.bib41)\)\. In Fig\.[3](https://arxiv.org/html/2608.17135#S2.F3)\(c\), we show the minimum bond dimension required for the truncation\-error budget to be satisfied \(blue\) and that required for the BF maximum search to succeed throughoutP=⌈log2⁡\(n​M\)⌉\+1P=\\lceil\\log\_\{2\}\(nM\)\\rceil\+1steps for Max\-33\-SAT problems withM=α​nM=\\alpha nclauses for different values of the hardness parameterα\\alpha\. As expected from the P≠\\neqNP conjecture, the minimum bond dimension for certification grows exponentially with the number of variables, saturating the maximum bond dimension bound \(black dotted\)2n/22^\{n/2\}\.

Finally, we observe that minimum bond dimensions required for successful BF are far smaller than the certification threshold, which explains why then=70n=70Max\-3\-SAT problem could be solved successfully despite the certification threshold being out of reach\. These two results presented in Fig\.[3](https://arxiv.org/html/2608.17135#S2.F3)\(a\) and \(c\) hint at a global optimization that requires far fewer computational resources than brute\-force search with a cost that is, for all practical purposes, polynomial innn\. The resource advantage applies to finding solutions, while verification of optimality may still require exponential resources\. The presented ITNT\-based methods thus shed light on a new computational paradigm offering a novel approach to combinatorial optimization\.

## IIIDiscussion

In this article, we have established iterative tensor network transformations \(ITNTs\) as a general framework for performing element\-wise computations on data compressed as a tensor train \(TT\)\. This work fundamentally broadens the scope of tensor network applications, introducing a method to evaluate arbitrary nonlinear functions entirely within the compressed domain\. Beyond field processing, ITNT introduces a novel optimization paradigm distinct from both \(quantum\) annealing and conventional local heuristic algorithms\. By comparing partial sums over exponentially large configuration subspaces at each step, the extremum\-finding algorithm performs a truly global search for optimization tasks where exact and heuristic methods often struggle with rugged, multi\-modal energy landscapes\. This demonstrates that ITNT provides high\-quality approximate solutions to NP\-hard problems within polynomial time under a controlled maximum bond dimension, while remaining informed by the solution landscape as a whole\.

The performance and ultimate limits of ITNTs represent a rich area for future investigation\. While the zip\-up method employed in this article is robust, alternative approaches may offer superior scaling in different regimes of system size and bond dimension\([43](https://arxiv.org/html/2608.17135#bib.bib43);[20](https://arxiv.org/html/2608.17135#bib.bib20);[44](https://arxiv.org/html/2608.17135#bib.bib44)\)\. Finally, the convergence of iterative routines could be significantly accelerated through tailored preprocessing steps, such as initial normalization or the application of a Fourier transform to the TT\([45](https://arxiv.org/html/2608.17135#bib.bib45)\)\. We anticipate that combining these algorithmic refinements with high\-precision arithmetic and GPU acceleration\([15](https://arxiv.org/html/2608.17135#bib.bib15);[46](https://arxiv.org/html/2608.17135#bib.bib46)\)will enable the application of ITNTs to problems on a much larger scale\.

This work opens a pathway for applying tensor network methods to a broader class of computationally intensive optimization problems, including those mappable to Ising Hamiltonians, such as prime factorization\([47](https://arxiv.org/html/2608.17135#bib.bib47)\), genetic haplotype reconstruction\([48](https://arxiv.org/html/2608.17135#bib.bib48)\), and large\-scale financial simulations\([49](https://arxiv.org/html/2608.17135#bib.bib49)\)\. One long\-term direction is to investigate the mapping of ITNT algorithms onto quantum hardware, where TT\-based methods are already showing significant promise\([50](https://arxiv.org/html/2608.17135#bib.bib50);[51](https://arxiv.org/html/2608.17135#bib.bib51);[52](https://arxiv.org/html/2608.17135#bib.bib52)\)\. Furthermore, the ability to encode symmetries directly within the TT formalism\([53](https://arxiv.org/html/2608.17135#bib.bib53)\)could be integrated with ITNTs to enforce physical or combinatorial constraints in complex optimization tasks, further enhancing the power and reach of this computational paradigm\.

## IVAcknowledgements

XW appreciates the insightful comments from Pan Zhang and Li You\. DJ and TH acknowledges support by the European Union’s Horizon Programme \(HORIZON\-CL42021DIGITALEMERGING\-02\-10\) Grant Agreement 101080085 QCFD\. DJ and TH are partially funded by the Cluster of Excellence ‘Advanced Imaging of Matter’ of the Deutsche Forschungsgemeinschaft \(DFG\)–EXC 2056\-project ID 390715994\. PS acknowledges financial support by the DLR\-Quantum\-Fellowship Program\. DJ acknowledges support by DFG project “Quantencomputing mit neutralen Atomen” \(JA 1793/1\-1, Japan\-JST\-DFG\-ASPIRE 2024\) and the Hamburg Quantum Computing Initiative \(HQIC\) project EFRE\. The EFRE project is co\-financed by ERDF of the European Union and by the “Fonds of the Hamburg Ministry of Science, Research, Equalities and Districts \(BWFGB\)”\.

## VMethods

In this section, we provide additional details on the foundations and applications of ITNTs discussed in the main text\. First,[subsection V\.1](https://arxiv.org/html/2608.17135#S5.SS1)describes the SVD\-based truncation procedure and sweeping protocol used to maintain low\-rank tensor representations\. Second,[subsection V\.2](https://arxiv.org/html/2608.17135#S5.SS2)outlines the partial integration formalism for tensor trains\. Third,[subsection V\.3](https://arxiv.org/html/2608.17135#S5.SS3)explains how to compute elementary functions with ITNT using the Newton\-Raphson method\. Fourth,[subsection V\.4](https://arxiv.org/html/2608.17135#S5.SS4)describes preconditioning techniques that accelerate the convergence of element\-wise sign function evaluation\. Fifth,[subsection V\.5](https://arxiv.org/html/2608.17135#S5.SS5)details the extremum\-finding ITNT scheme using self\-multiplication, binary fixation, and deflation\. Finally,[subsection V\.6](https://arxiv.org/html/2608.17135#S5.SS6)explains how the energy landscape of a classical Ising model is encoded in a compressed TT representation\.

### V\.1The SVD truncation

The origin of compression in the TT representation lies in the singular value decomposition \(SVD\), which identifies the optimal low\-rank approximation of a tensor by minimizing theℓ2\\ell\_\{2\}norm of the truncation error\([2](https://arxiv.org/html/2608.17135#bib.bib2)\)\. By retaining only the subspace that corresponds to the most significant singular values, SVD truncation ensures the highest\-fidelity approximation for a fixed bond dimensionχ\\chi\. The careful choice ofχ\\chiensures numerical stability and efficient convergence\. Furthermore, when performing two\-site SVD truncation on a TT representing a smooth function, it is crucial to start the SVD sweep from the smallest scale, i\.e\., truncating from the least significant bit, to achieve efficient compression\. Notably, reversing the sequence of this sweeping process \(i\.e\., fixing the bond dimension but starting the truncation from the large\-scale side\) leads to a significant increase in the truncation error, emphasizing the critical role of the chosen order\.

### V\.2Partial Integration

The partial integrationℐ⁡\(Ψ,q\)\\mathcal\{I\}\(\\Psi,q\)gives a two\-element array, whereℐ​\(Ψ,q\)\[0\]\\mathcal\{I\}\(\\Psi,q\)^\{\[0\]\}is the sum of all elements ofΨ\[j\]\\Psi^\{\[j\]\}with theqqthbit of index𝒋\\bm\{j\}equal to00, andℐ​\(Ψ,q\)\[1\]\\mathcal\{I\}\(\\Psi,q\)^\{\[1\]\}is the corresponding sum for the bit equal to11\. Formally,ℐ⁡\(Ψ,q\)\\mathcal\{I\}\(\\Psi,q\)is written as

ℐ\(Ψ,q\)\[jq\]=∑𝒋\\\{jq\}∑𝜶Ψ\[0\]α0,α1\[j0\]⋯Ψ\[q\]αq,αq\+1\[jq\]⋯Ψ\[n−1\]αn−1,αn\[jn−1\]\\displaystyle\\mathcal\{I\}\(\\Psi,q\)^\{\[j\_\{q\}\]\}=\\sum\_\{\\bm\{j\}\\backslash\\\{j\_\{q\}\\\}\}\\sum\_\{\\bm\{\\alpha\}\}\\Psi\[0\]\_\{\\alpha\_\{0\},\\alpha\_\{1\}\}^\{\[j\_\{0\}\]\}\\cdots\\Psi\[q\]\_\{\\alpha\_\{q\},\\alpha\_\{q\+1\}\}^\{\[j\_\{q\}\]\}\\cdots\\Psi\[n\-1\]\_\{\\alpha\_\{n\-1\},\\alpha\_\{n\}\}^\{\[j\_\{n\-1\}\]\}\(6\)where𝒋\\\{jq\}\\bm\{j\}\\backslash\\\{j\_\{q\}\\\}denotes the bit string𝒋\\bm\{j\}without theqqthbit\.

### V\.3ITNT for realizing elementary functions

We define the elementary functions as those specified by Bourchteinet al\.\([31](https://arxiv.org/html/2608.17135#bib.bib31)\), comprising polynomial, rational, exponential, trigonometric, and irrational functions \(i\.e\., non\-integer powers and logarithms\)\. First, given the TTffrepresenting the discretized functionf⁡\(x\)f\(x\)overN=2nN=2^\{n\}grid points, we construct the TTR\(k\)R^\{\(k\)\}which represents thekkthinteger power off⁡\(x\)f\(x\),\(f⁡\(x\)\)k\\left\(f\(x\)\\right\)^\{k\}, using the following iterative element\-wise product

R\(k\+1\)=f⊙R\(k\),\\displaystyle R^\{\(k\+1\)\}=f\\odot R^\{\(k\)\},\(7\)whereR\(1\)=fR^\{\(1\)\}=f, and the operationh=f⊙gh=f\\odot gon TTsffandgg, denotes an element\-wise product, yielding a resulting TThhwith elementsh\[j\]=f\[j\]​g\[j\]h^\{\[j\]\}=f^\{\[j\]\}g^\{\[j\]\}\. With this choice of transform, we obtainR\(k\)R^\{\(k\)\}, which is equivalent tof​\(x\)⊙kf\(x\)^\{\\odot k\}\. We assume that the bond dimensions of bothffandR\(k\)R^\{\(k\)\}areχ\\chi, and this element\-wise product is realized using the zip\-up method by transforming one of them into a diagonal operator using a copy tensor\([29](https://arxiv.org/html/2608.17135#bib.bib29);[30](https://arxiv.org/html/2608.17135#bib.bib30)\)\. Here, the maximum bond dimension is truncated toχ\\chi\. Under this assumption, the element\-wise product has a time complexity of𝒪⁡\(n​χ4\)\\mathcal\{O\}\(n\\chi^\{4\}\)and a memory complexity of𝒪⁡\(χ3\)\\mathcal\{O\}\(\\chi^\{3\}\)\.

To compute the element\-wise exponentiation of the field,f⁡\(x\)→ef⁡\(x\)f\(x\)\\to e^\{f\(x\)\}, we rely on the definition of exponentiation as a limit

ef⁡\(x\)=liml→∞\(1\+f⁡\(x\)2l\)⊙2l,\\displaystyle e^\{f\(x\)\}=\\lim\_\{l\\to\\infty\}\\left\(1\+\\frac\{f\(x\)\}\{2^\{l\}\}\\right\)^\{\\odot 2^\{l\}\},\(8\)wherel∈ℤ\+l\\in\\mathbb\{Z\}^\{\+\}\. We can construct this iteration in the TT representation by applying the following transformationlltimes,

R\(k\+1\)=F⁡\(R\(k\)\)=R\(k\)⊙R\(k\),\\displaystyle R^\{\(k\+1\)\}=F\\left\(R^\{\(k\)\}\\right\)=R^\{\(k\)\}\\odot R^\{\(k\)\},\(9\)starting fromR\(0\)=𝕀\+f/2lR^\{\(0\)\}=\\mathbb\{I\}\+f/2^\{l\}\. Here, the unit\-element TT𝕀\\mathbb\{I\}\(i\.e\., representing a vector with all its elements equal to11\) is constructed from the tensor product of rank\-1 unit tensors, and thus has a bond dimension of11\. Asllincreases, the final TT approaches the exponentiated TT,liml→∞R\(l\)=exp⁡\(f\)\\lim\_\{l\\to\\infty\}R^\{\(l\)\}=\\exp\\left\(f\\right\), which representsef⁡\(x\)e^\{f\(x\)\}\. Although the error in Eq\. \([8](https://arxiv.org/html/2608.17135#S5.E8)\) becomes negligible when2l≫fmax≡maxx⁡f⁡\(x\)2^\{l\}\\gg f\_\{\\text\{max\}\}\\equiv\\max\_\{x\}f\(x\), we still need to restrict the number of iteration stepsllto avoid accumulating large truncation errors from the element\-wise products in the ITNT framework\.

This exponentiation ITNT enables the evaluation of thermal path integrals for multi\-particle systems, ranging from discrete spin models to continuous\-space chemical dynamics, directly in the TT representation\. Specifically, we evaluate the sum ofeA⁡\(xt\)e^\{A\(x\_\{t\}\)\}over all configurationsxtx\_\{t\}, where each configuration represents a distinct imaginary\-time\-discretized path\. Using similar iterative self\-multiplication techniques, we can also construct the unitary phase transformationf⁡\(x\)→ei​f​\(x\)f\(x\)\\to e^\{\\mathrm\{i\}f\(x\)\}for complex\-valued TT representations, which serves as a foundation for realizing sine, cosine, and arcsine transforms\. While the convergence of the self\-multiplication method may become slow for highly oscillatory functions, the process can be significantly accelerated using the CORDIC algorithm\([54](https://arxiv.org/html/2608.17135#bib.bib54)\)\. The CORDIC framework decomposes complex functional evaluations into a sequence of elementary, pre\-computed rotations\. By utilizing only bit shifts \(scalar multiplications by negative integer powers of 2\), additions, and comparisons \(as discussed in this article\), CORDIC avoids the computational overhead of high\-order polynomial expansions\. This allows us to perform the element\-wise evaluation of trigonometric and hyperbolic functions within the TT representation \(cf\. Ref\.\([54](https://arxiv.org/html/2608.17135#bib.bib54)\)for a modern implementation\)\.

We further expand the ITNT framework by incorporating the Newton\-Raphson method\. The method obtains the solutionRRto the equationg⁡\(R\)=0g\(R\)=0iteratively according to

R\(k\+1\)=R\(k\)−g⁡\(R\(k\)\)g′​\(R\(k\)\),R^\{\(k\+1\)\}=R^\{\(k\)\}\-\\frac\{g\(R^\{\(k\)\}\)\}\{g^\{\\prime\}\(R^\{\(k\)\}\)\},\(10\)whereg′​\(R\)g^\{\\prime\}\(R\)is the derivative ofggwith respect toRR\. The Newton\-Raphson method exhibits quadratic convergence as long as the functiong⁡\(R\)g\(R\)is differentiable at its root\.

This Newton\-Raphson method can be used to perform an elementary nonlinear transformationf⁡\(x\)→𝒯⁡\(f⁡\(x\)\)f\(x\)\\to\\mathcal\{T\}\(f\(x\)\)at each point of the initial functionf⁡\(x\)f\(x\)\. To achieve this for an arbitrary functionR⁡\(x\)R\(x\), we define

g⁡\(R⁡\(x\)\)≡𝒯−1​\(R⁡\(x\)\)−f⁡\(x\),g\(R\(x\)\)\\equiv\\mathcal\{T\}^\{\-1\}\(R\(x\)\)\-f\(x\),\(11\)where𝒯−1​\(⋅\)\\mathcal\{T\}^\{\-1\}\(\\cdot\)denotes the inverse transformation of𝒯⁡\(⋅\)\\mathcal\{T\}\(\\cdot\)\. Thus, the solution tog⁡\(R⁡\(x\)\)=0g\(R\(x\)\)=0is given byR⁡\(x\)=𝒯⁡\(f⁡\(x\)\)R\(x\)=\\mathcal\{T\}\(f\(x\)\)\. However, this solutionR⁡\(x\)R\(x\)can also be determined as the fixed point of the Newton\-Raphson iterative routine in Eq\. \([10](https://arxiv.org/html/2608.17135#S5.E10)\)\. The iteration is given by

R\(k\+1\)​\(x\)=R\(k\)​\(x\)−𝒯−1​\(R\(k\)​\(x\)\)−f⁡\(x\)\(𝒯−1\)′​\(R\(k\)​\(x\)\),R^\{\(k\+1\)\}\(x\)=R^\{\(k\)\}\(x\)\-\\frac\{\\mathcal\{T\}^\{\-1\}\(R^\{\(k\)\}\(x\)\)\-f\(x\)\}\{\(\\mathcal\{T\}^\{\-1\}\)^\{\\prime\}\(R^\{\(k\)\}\(x\)\)\},\(12\)where\(𝒯−1\)′\(\\mathcal\{T\}^\{\-1\}\)^\{\\prime\}denotes the derivative of𝒯−1\\mathcal\{T\}^\{\-1\}with respect toR\(k\)​\(x\)R^\{\(k\)\}\(x\)\. Thus, we find that our targeted transformation is given by𝒯⁡\(f⁡\(x\)\)=limk→∞R\(k\)​\(x\)\\mathcal\{T\}\(f\(x\)\)=\\lim\_\{k\\to\\infty\}R^\{\(k\)\}\(x\), provided that the initial functionR\(0\)​\(x\)R^\{\(0\)\}\(x\)is properly chosen to ensure convergence\.

The above Newton\-Raphson equation \(Eq\. \([12](https://arxiv.org/html/2608.17135#S5.E12)\)\) can be simplified and then iteratively computed in the TT representation from the knownf⁡\(x\)f\(x\)and a properly chosenR\(0\)R^\{\(0\)\}\. We explain this procedure in the following examples, where the target nonlinear transformation𝒯⁡\(⋅\)\\mathcal\{T\}\(\\cdot\)respectively represents the inverse, square\-root, and logarithmic transformations\.

We use the Newton\-Raphson method to compute the element\-wise reciprocal,R⁡\(x\)=1/f⁡\(x\)R\(x\)=1/f\(x\), of the functionf⁡\(x\)f\(x\)using ITNT\. To this end, the equation is reformulated as a root\-finding procedure forg⁡\(R⁡\(x\)\)=\(R⁡\(x\)\)−1−f⁡\(x\)=0g\(R\(x\)\)=\(R\(x\)\)^\{\-1\}\-f\(x\)=0, defined in Eq\. \([11](https://arxiv.org/html/2608.17135#S5.E11)\)\. The resulting Newton\-Raphson iteration becomesR\(k\+1\)​\(x\)=R\(k\)​\(x\)​\(−f⁡\(x\)​R\(k\)​\(x\)\+2\)R^\{\(k\+1\)\}\(x\)=R^\{\(k\)\}\(x\)\(\-f\(x\)R^\{\(k\)\}\(x\)\+2\)\. This can be rewritten in the TT representation, using the element\-wise product and the unit\-element TT𝕀\\mathbb\{I\}, as follows

R\(k\+1\)=R\(k\)⊙\(−f⊙R\(k\)\+2𝕀\)\.R^\{\(k\+1\)\}=R^\{\(k\)\}\\odot\\left\(\-f\\odot R^\{\(k\)\}\+2\\mathbb\{I\}\\right\)\.\(13\)Here, we chooseR\(0\)=f/cR^\{\(0\)\}=f/cas an initial guess, whereccis a suppression constant that ensures that the iterations converge\.

To compute the square root,R⁡\(x\)=f⁡\(x\)R\(x\)=\\sqrt\{f\(x\)\}, of a functionf⁡\(x\)f\(x\)using ITNT, we again follow Eq\. \([11](https://arxiv.org/html/2608.17135#S5.E11)\) and defineg⁡\(R\)=R2−f⁡\(x\)=0g\(R\)=R^\{2\}\-f\(x\)=0\. The corresponding Newton\-Raphson iteration becomesR\(k\+1\)​\(x\)=12​\(R\(k\)​\(x\)\+f⁡\(x\)R\(k\)​\(x\)\)R^\{\(k\+1\)\}\(x\)=\\frac\{1\}\{2\}\\left\(R^\{\(k\)\}\(x\)\+\\frac\{f\(x\)\}\{R^\{\(k\)\}\(x\)\}\\right\), which translates into the TT representation as

R\(k\+1\)=12​\(R\(k\)\+f⊙\(1/R\(k\)\)\)\.R^\{\(k\+1\)\}=\\frac\{1\}\{2\}\\left\(R^\{\(k\)\}\+f\\odot\\left\(1/R^\{\(k\)\}\\right\)\\right\)\.\(14\)This is also known as Heron’s method\. These two operations enable the element\-wise evaluation of fractional and negative powers of a field\. Thus, by iteratively applying these operations, we can achieve the element\-wise evaluation of arbitrary real powers\.

The natural logarithm,R⁡\(x\)=ln⁡f⁡\(x\)R\(x\)=\\ln f\(x\), can be computed by definingg⁡\(R\)=eR−f⁡\(x\)g\(R\)=e^\{R\}\-f\(x\), resulting in the iterative transformationR\(k\+1\)​\(x\)=R\(k\)​\(x\)−1\+f⁡\(x\)​e−R\(k\)​\(x\)R^\{\(k\+1\)\}\(x\)=R^\{\(k\)\}\(x\)\-1\+f\(x\)e^\{\-R^\{\(k\)\}\(x\)\}, which translates into the TT representation as follows

R\(k\+1\)=R\(k\)\+f⊙exp⁡\(−R\(k\)\)−𝕀\.R^\{\(k\+1\)\}=R^\{\(k\)\}\+f\\odot\\exp\\left\(\-R^\{\(k\)\}\\right\)\-\\mathbb\{I\}\.\(15\)We start fromR\(0\)=0R^\{\(0\)\}=0\. This ensures thatln⁡\(f⁡\(x\)\)=R\(∞\)​\(x\)\\ln\(f\(x\)\)=R^\{\(\\infty\)\}\(x\)\. We then compute the exponentiated TT,exp⁡\(−R\(k\)\)\\exp\\left\(\-R^\{\(k\)\}\\right\), using the ITNT algorithm defined in Eq\. \([8](https://arxiv.org/html/2608.17135#S5.E8)\)\. Convergence is guaranteed whenf⁡\(x\)\>0f\(x\)\>0for allxx\.

### V\.4Accelerating element\-wise sign function evaluation

It was shown for the sign transformation\([36](https://arxiv.org/html/2608.17135#bib.bib36)\)that applying the transformation

G⁡\(R\(k\)\)=−1210​\(1359​\(\(R\(k\)\)⊙3\)−2126​R\(k\)\)\\displaystyle G\(R^\{\(k\)\}\)=\-\\frac\{1\}\{2^\{10\}\}\\left\(1359\(\(R^\{\(k\)\}\)^\{\\odot 3\}\)\-2126R^\{\(k\)\}\\right\)\(16\)before iterating withFF\(Eq\. \([1](https://arxiv.org/html/2608.17135#S2.E1)\)\) further accelerates convergence\. In practice, applying the auxiliary transformationGGtwice serves as an optimal preconditioning step, allowing the TTffto converge more rapidly toSign​\{f\}\\mathrm\{Sign\}\\\{f\\\}\. The numerical performance of this iterative composition is illustrated in Fig\.[4](https://arxiv.org/html/2608.17135#S5.F4)\(b\)\. Starting from the identity functionf⁡\(x\)=xf\(x\)=x, the compositionF\(3\)∘G\(2\)F^\{\(3\)\}\\circ G^\{\(2\)\}rapidly develops a sharp, step\-like profile near the origin, effectively approximating the sign function within the compressed TT manifold\.

As depicted in the figure, the transformation maintains high fidelity for input values within the stable domain; however, for\|x\|\|x\|exceeding the radius of convergence \(indicated by the red shaded regions\), the iterative map exhibits unstable oscillatory behavior and eventual divergence\. This divergence occurs because the polynomial compositionF∘GF\\circ Gis an iterative map with fixed points at\{−1,0,1\}\\\{\-1,0,1\\\}\. Input values outside the basin of attraction of these points are driven towards infinity, necessitating the global rescaling protocol discussed previously to ensure that all elements of the field remain within the convergence regime\.

### V\.5Extremum\-finding ITNT with a limited number of self\-multiplications

![Refer to caption](https://arxiv.org/html/2608.17135v1/Fig2.png)Figure 4:ITNTs for element\-wise nonlinear operations in TT representation\.\(a\) Extrema detection with TT\. i Starting from a TT representation ofΨ\\Psi\(black curve\), repeated element\-wise self\-multiplication sharpens the state into a multi\-peak profileΦ\\Phi, where its integral is dominated by extremal configurations \(red curve\)\. ii A peak is then located by a binary fixation procedure\. This involves integrating out all unfixed sites and choosing the local configuration that maximizes the summed value\. Starting from every index not being fixed, in the leftmost panel, a comparison betweenj2=0j\_\{2\}=0\(orange\) andj2=1j\_\{2\}=1\(cyan\) is performed, identifyingj2=1j\_\{2\}=1as a candidate for the optimal configuration\. In the subsequent next panel,j1=0j\_\{1\}=0\(orange\) andj1=1j\_\{1\}=1\(cyan\) are compared andj1=0j\_\{1\}=0is selected, while the integrated out region from the fixation ofj2=0j\_\{2\}=0is indicated in gray\. Continuing these steps with a random sequence of local tensors identifies the indexIIof the maximum \(orange, right most panel\)\. iii The deflation process then removes the identified peak to ensure that subsequent applications find the true maximum\. \(b\) Element\-wise evaluation of sign function\. A highly nonlinear, non\-analytic element\-wise filter is realized through the composition of iterative applications of the functionsFF\(Eq\. \([1](https://arxiv.org/html/2608.17135#S2.E1)\)\) andGG\(Eq\. \([16](https://arxiv.org/html/2608.17135#S5.E16)\)\) that produce an element\-wise sign functionSign​\{f\}\[j\]=Sign⁡\{f\[j\]\}\\mathrm\{Sign\}\\\{f\\\}^\{\[j\]\}=\\mathrm\{Sign\}\\\{f^\{\[j\]\}\\\}\. As an illustrative example, the blue curve shows the output of ITNT for the input fieldf⁡\(x\)=xf\(x\)=x\(black dashed line\) after applyingGGtwice followed by repeated iterations ofFF\. The red shaded region marks the regime of input values where the composed iterationF⁡\(…​F​\(G⁡\(G⁡\(f\)\)\)\)F\(\\ldots F\(G\(G\(f\)\)\)\)does not converge\. The curves are evaluated in the full\-array representation solely to visualize the mechanism of the procedure\.Following the main text and Fig\.[4](https://arxiv.org/html/2608.17135#S5.F4)\(a\), in this subsection, we describe the ITNT algorithm for extrema location in more detail\. Given annn\-tensor TTffencoding a real\-valued field, we amplify large elements by iterative element\-wise self\-multiplication,

R\(k\+1\)\\displaystyle R^\{\(k\+1\)\}=R\(k\)⊙R\(k\)\\displaystyle=R^\{\(k\)\}\\odot R^\{\(k\)\}\(R\(0\)=f\)\.\\displaystyle\(R^\{\(0\)\}=f\)\.\(17\)We then apply normalization after each step,

R\(k\)←R\(k\)‖R\(k\)‖\.R^\{\(k\)\}\\leftarrow\\frac\{R^\{\(k\)\}\}\{\\\|R^\{\(k\)\}\\\|\}\.\(18\)Afterkkiterations,R\(k\)R^\{\(k\)\}representsf⊙2kf^\{\\odot 2^\{k\}\}up to normalization and typically develops a small set of dominant peaks \(cf\. Fig\.[4](https://arxiv.org/html/2608.17135#S5.F4)\(a\) i\.\)\. Since self\-multiplication squares the elements, the peaks appear at indices where\|f\[j\]\|\|f^\{\[j\]\}\|is large\.

#### Peak location with binary fixation

To extract a candidate indexIIfrom a multi\-peak TTRR, we randomly choose a site ordering\(p0,…,pn−1\)\(p\_\{0\},\\dots,p\_\{n\-1\}\)\. At a binary fixation steptt, with sites\{p0,…,pt−1\}\\\{p\_\{0\},\\dots,p\_\{t\-1\}\\\}already fixed, we compare the magnitudes of the two partial integrations forpt=0p\_\{t\}=0andpt=1p\_\{t\}=1, after all other free tensors have been summed out\. We then setptp\_\{t\}to the branch with the larger magnitude, and continue this process until all sites are fixed, yielding a candidate indexII\(Fig\.[4](https://arxiv.org/html/2608.17135#S5.F4)\(a\)ii\)\.

#### Deflation

Once a candidate indexIIyielding a maximum element ofRRis obtained, we remove the corresponding peak fromRRby subtracting its contribution \(Fig\.[4](https://arxiv.org/html/2608.17135#S5.F4)\(a\)iii\)

R←R−R\[I\]​δI,R\\leftarrow R\-R^\{\[I\]\}\\delta\_\{I\},\(19\)whereδI\\delta\_\{I\}is aχ=1\\chi=1TT whose local tensors are given byδ​\[q\]0,0\[i\]=δiq,Iq\\delta\[q\]^\{\[i\]\}\_\{0,0\}=\\delta\_\{i\_\{q\},I\_\{q\}\}\. Hereδa,b\\delta\_\{a,b\}is a Kronecker delta, defined asδa,b=0\\delta\_\{a,b\}=0\(a≠ba\\neq b\) andδa,b=1\\delta\_\{a,b\}=1\(a=ba=b\)\. Iterating binary fixation→\\todeflation produces a ranked list of candidate maxima, from which the largest is accepted as the true maximum\.

In principle, choosing a sufficient number of self\-multiplications followed by BF allows one to identify the true maximum with certainty\. However, the iterative self\-multiplication step leads to an exponential growth in the bond dimension\. Thus, intermediate truncations of the bond dimension are necessary but will introduce truncation errors in most cases\. We find that the accumulation of truncation errors and round\-off errors in double\-precision arithmetic can prevent the ITNT algorithm from reliably locating the global extremum, especially if the energy landscape features near\-degenerate extrema\.

This can be mitigated by randomly traversing the site order across different binary fixation steps and performing deflation\. The randomized binary fixation sequence at each iteration prevents the algorithm from getting stuck in local extrema\. While a fixed sequence like\(0,…,n−1\)\(0,\\dots,n\-1\)tends to yield suboptimal peaks with redundant leading configurations, the stochastic selection of the fixation sites ensures that the global maximum eventually emerges, overcoming earlier truncation and round\-off errors \(see SM[IV](https://arxiv.org/html/2608.17135#S4a)for an application in energy matching in the Ising model\)\.

### V\.6Encoding the energy landscape of Ising Hamiltonians

In this section, we explain how the energy landscape of the classical Ising model can be encoded into a TT\. We define the classical Ising model onnnsites with the Hamiltonian

H=∑i=1n∑j=1iJi​j​σi​σj,H=\\sum\_\{i=1\}^\{n\}\\sum\_\{j=1\}^\{i\}J\_\{ij\}\\sigma\_\{i\}\\sigma\_\{j\},\(20\)whereσi∈\{0,1\}\\sigma\_\{i\}\\in\\\{0,1\\\}are the spin variables\. Here, we consider Hamiltonians with only two\-body interactions; however, the result trivially generalizes to Hamiltonians with arbitrarymm\-body interactions\. For a given configuration𝝈=\(σ0,σ1,…,σn−1\)\\bm\{\\sigma\}=\(\\sigma\_\{0\},\\sigma\_\{1\},\\ldots,\\sigma\_\{n\-1\}\), the energy is given by

E⁡\(𝝈\)=∑i=1n∑j=1iJi​j​σi​σj\.E\(\\bm\{\\sigma\}\)=\\sum\_\{i=1\}^\{n\}\\sum\_\{j=1\}^\{i\}J\_\{ij\}\\sigma\_\{i\}\\sigma\_\{j\}\.\(21\)Here,Ji​jJ\_\{ij\}denotes the Ising couplings fori≠ji\\neq j, andJi​iJ\_\{ii\}represents the on\-site terms, sinceσi2=σi\\sigma\_\{i\}^\{2\}=\\sigma\_\{i\}\.

To find the configuration that maximizesE⁡\(𝝈\)E\(\\bm\{\\sigma\}\)by performing the extremum\-finding ITNT, we encode the energy landscape into the TTEEby setting the values ofEEat indexσ\\sigmato the energyE⁡\(𝝈\)E\(\\bm\{\\sigma\}\), i\.e\.,E\[σ\]=E⁡\(𝝈\)E^\{\[\\sigma\]\}=E\(\\bm\{\\sigma\}\), whereσ\\sigmais the index obtained by interpreting the configuration as a bit string\. The TT representation of the energy landscapeE⁡\(𝝈\)E\(\\bm\{\\sigma\}\)is sparse\. The bond dimensionχ\\chiofEEis bounded by the number of non\-zero coupling terms inJi​jJ\_\{ij\}, because we can decomposeEEasE=∑i=1n∑j=1iE⁡\[i,j\]E=\\sum\_\{i=1\}^\{n\}\\sum\_\{j=1\}^\{i\}E\[i,j\], whereE⁡\[i,j\]E\[i,j\]is a TT encoding the coupling energy corresponding to the termJi​jJ\_\{ij\}across all configurations\. Crucially,E⁡\[i,j\]E\[i,j\]can be recast into aχ=1\\chi=1TT as follows, fori≠ji\\neq j,

E\[i,j\]\[k\]=Ji​j∑𝜶\\displaystyle E\[i,j\]^\{\[k\]\}=J\_\{ij\}\\sum\_\{\\bm\{\\alpha\}\}A\[0\]α0,α1\[k0\]⋯A\[i−1\]αi−1,αi\[ki−1\]B\[i\]αi,αi\+1\[ki\]A\[i\+1\]αi\+1,αi\+2\[ki\+1\]⋯\\displaystyle A\[0\]\_\{\\alpha\_\{0\},\\alpha\_\{1\}\}^\{\[k\_\{0\}\]\}\\cdots A\[i\-1\]\_\{\\alpha\_\{i\-1\},\\alpha\_\{i\}\}^\{\[k\_\{i\-1\}\]\}B\[i\]\_\{\\alpha\_\{i\},\\alpha\_\{i\+1\}\}^\{\[k\_\{i\}\]\}A\[i\+1\]\_\{\\alpha\_\{i\+1\},\\alpha\_\{i\+2\}\}^\{\[k\_\{i\+1\}\]\}\\cdotsA\[j−1\]αj−1,αj\[kj−1\]B\[j\]αj,αj\+1\[kj\]A\[j\+1\]αj\+1,αj\+2\[kj\+1\]⋯A\[n−1\]αn−1,αn\[kn−1\],\\displaystyle A\[j\-1\]\_\{\\alpha\_\{j\-1\},\\alpha\_\{j\}\}^\{\[k\_\{j\-1\}\]\}B\[j\]\_\{\\alpha\_\{j\},\\alpha\_\{j\+1\}\}^\{\[k\_\{j\}\]\}A\[j\+1\]\_\{\\alpha\_\{j\+1\},\\alpha\_\{j\+2\}\}^\{\[k\_\{j\+1\}\]\}\\cdots A\[n\-1\]\_\{\\alpha\_\{n\-1\},\\alpha\_\{n\}\}^\{\[k\_\{n\-1\}\]\},\(22\)whereA⁡\[q\]A\[q\]is the unit tensor with all elements equal to11, andB⁡\[q\]B\[q\]is aχ=1\\chi=1local tensor which takes the formB​\[q\]0,0\[kq\]=kqB\[q\]^\{\[k\_\{q\}\]\}\_\{0,0\}=k\_\{q\}, for a given siteqq\. To encode the energy landscape of anymm\-body interaction term,E⁡\[i1,i2,i3,⋯,im\]E\[i\_\{1\},i\_\{2\},i\_\{3\},\\cdots,i\_\{m\}\], one creates the unit TT𝕀\\mathbb\{I\}and then replaces the local tensors corresponding to the sites that appear in the interaction term with theBBtensor\.

According to Ref\.\([2](https://arxiv.org/html/2608.17135#bib.bib2)\), the sum of two TTs yields another TT with a summed bond dimension\. Thus, summing overE⁡\[i,j\]E\[i,j\]results in the TTEEhaving a bond dimension equivalent to the number of terms inE⁡\(𝝈\)E\(\\bm\{\\sigma\}\), including the constant offset\. Because of this sparsity, we perform ITNTs directly onEEto extract various features of the energy landscape: the lowest\- and highest\-energy configurations, as well as configurations closest to a specified energy value\. We view this as a global transformation of the energy landscape in its TT representation, which acts simultaneously across all spin configurations\. This approach differs from tensor network simulations of quantum annealing\([55](https://arxiv.org/html/2608.17135#bib.bib55)\), where the TT encodes the quantum state of the annealer rather than the energy landscape considered here\.

### V\.7Truncation error budget for certified optimality of Max\-3\-SAT solutions

Finally, we derive a truncation\-error bound extending Eq\. \([2](https://arxiv.org/html/2608.17135#S2.E2)\) that certifies the maximum found by the maximum\-finding ITNT\. First we derive the sufficient numberPPof self\-multiplication steps for the fieldEEthat encodes the energy landscape of theMM\-clause Max\-3\-SAT problem as defined in Eq\. \([5](https://arxiv.org/html/2608.17135#S2.E5)\)\. We then define the normalized fieldΨ=E/A\\Psi=E/A, whereA=max⁡\{E\}A=\\max\\\{E\\\}is the maximum number of satisfied clauses\.

A value ofPPsufficient for BF to retain an optimal branch is obtained by considering the worst case in which every nonoptimal element ofΨ\\Psitakes its largest possible value,1−1/A1\-1/A\. LetJJbe the index withΨ\[J\]=1\\Psi^\{\[J\]\}=1, assumed unique for simplicity; degeneracy of the optimum only increases with the margin below\. The sum of all elements ofΨ⊙2P\\Psi^\{\\odot 2^\{P\}\}other thanJJthen obeys

∑j≠J\(Ψ⊙2P\)\[j\]≤2n\(1−1/A\)2P≤2nexp\(−2P/M\),\\displaystyle\\sum\_\{j\\neq J\}\(\\Psi^\{\\odot 2^\{P\}\}\)^\{\[j\]\}\\leq 2^\{n\}\(1\-1/A\)^\{2^\{P\}\}\\leq 2^\{n\}\\exp\\left\(\-2^\{P\}/M\\right\),\(23\)where the right\-hand inequality usesA≤MA\\leq M, valid since at mostMMclauses can be satisfied\. It is therefore sufficient to choosePPsuch that2nexp\(−2P/M\)<12^\{n\}\\exp\\left\(\-2^\{P\}/M\\right\)<1, and hence we obtainP\>log2⁡\(n​M​ln⁡2\)P\>\\log\_\{2\}\(nM\\ln 2\)\. For the result in Fig\.[3](https://arxiv.org/html/2608.17135#S2.F3)\(c\) we setP=⌈log2⁡\(n​M\)⌉\+1P=\\left\\lceil\\log\_\{2\}\(nM\)\\right\\rceil\+1, which obeys the bound\.

GivenPPsatisfies the bound, we now derive the truncation budget for this dominance condition to be met at the end ofPPself\-multiplication steps\. LetR\(p\)R^\{\(p\)\}denote the computed field afterp≤Pp\\leq Psteps, withR\(0\)=ΨR^\{\(0\)\}=\\Psiand

R\(p\)=\(R\(p−1\)\)⊙2\+τp\\displaystyle R^\{\(p\)\}=\(R^\{\(p\-1\)\}\)^\{\\odot 2\}\+\\tau\_\{p\}\(24\)whereτp\\tau\_\{p\}is the total error injected by the truncation at steppp\.

Lettp=ε⁡\(\(R\(p−1\)\)⊙2,R\(p\)\)=‖τp‖2t\_\{p\}=\\varepsilon\(\(R^\{\(p\-1\)\}\)^\{\\odot 2\},R^\{\(p\)\}\)=\\left\\\|\\tau\_\{p\}\\right\\\|\_\{2\}be theℓ2\\ell\_\{2\}norm of the error as defined in Eq\. \([2](https://arxiv.org/html/2608.17135#S2.E2)\) and letΔp=R\(p\)−Ψ⊙2p\\Delta\_\{p\}=R^\{\(p\)\}\-\\Psi^\{\\odot 2^\{p\}\}be the error tensor, the accumulated deviation from the true representation of the self\-multiplied field\.

AfterPPiterations, at each step BF compares two partial integrals over the sites not yet fixed and keeps the larger\. Suppose all the previous steps have identified the bits ofJJcorrectly\. The difference between the partial integral at theqqth site with the correct bit \(jq=Jqj\_\{q\}=J\_\{q\}\) and with the incorrect bit \(jq≠Jqj\_\{q\}\\neq J\_\{q\}\) is then bounded from below by

Λ\(Ψ⊙2P,q\)≥1−∑i≠J\(Ψ⊙2P\)\[i\]≥1−2nexp\(−2P/M\),\\displaystyle\\Lambda\(\\Psi^\{\\odot 2^\{P\}\},q\)\\geq 1\-\\sum\_\{i\\neq J\}\(\\Psi^\{\\odot 2^\{P\}\}\)^\{\[i\]\}\\geq 1\-2^\{n\}\\exp\(\-2^\{P\}/M\),\(25\)whereΛ⁡\(X,q\)=ℐ​\(X,q\)jq=Jq−ℐ​\(X,q\)jq≠Jq\\Lambda\(X,q\)=\\mathcal\{I\}\(X,q\)^\{j\_\{q\}=J\_\{q\}\}\-\\mathcal\{I\}\(X,q\)^\{j\_\{q\}\\neq J\_\{q\}\}for a given fieldXXis the difference between the branches, one with the solution and the other without, and the second inequality holds sinceJJcontributes11to the first integral and every remaining configuration contributes non\-negatively to the second\. Invoking the triangular inequality, forR\(p\)R^\{\(p\)\}to identifyJJcorrectly it therefore suffices that

Λ\(R\(P\),q\)≥1−2nexp\(−2P/M\)−\|Λ\(ΔP,q\)\|\>0\.\\displaystyle\\Lambda\(R^\{\(P\)\},q\)\\geq 1\-2^\{n\}\\exp\(\-2^\{P\}/M\)\-\\left\|\\Lambda\(\\Delta\_\{P\},q\)\\right\|\>0\.\(26\)
From the definition,Δp\\Delta\_\{p\}obeys the recursion

Δp=2​Ψ⊙2p−1⊙Δp−1\+Δp−1⊙2\+τp\.\\displaystyle\\Delta\_\{p\}=2\\Psi^\{\\odot 2^\{p\-1\}\}\\odot\\Delta\_\{p\-1\}\+\\Delta^\{\\odot 2\}\_\{p\-1\}\+\\tau\_\{p\}\.\(27\)From this, provided that the initial truncation is00, we derive

ΔP\\displaystyle\\Delta\_\{P\}=2P−p​Ψ⊙\(2P−2p\)⊙Δp\+∑m=p\+1P2P−m​Ψ⊙\(2P−2m\)⊙Δm−1⊙2\+∑m=p\+1P−12P−m​Ψ⊙\(2P−2m\)⊙τm\\displaystyle=2^\{P\-p\}\\Psi^\{\\odot\(2^\{P\}\-2^\{p\}\)\}\\odot\\Delta\_\{p\}\+\\sum\_\{m=p\+1\}^\{P\}2^\{P\-m\}\\Psi^\{\\odot\(2^\{P\}\-2^\{m\}\)\}\\odot\\Delta\_\{m\-1\}^\{\\odot 2\}\+\\sum\_\{m=p\+1\}^\{P\-1\}2^\{P\-m\}\\Psi^\{\\odot\(2^\{P\}\-2^\{m\}\)\}\\odot\\tau\_\{m\}=∑m=1P2P−m​Ψ⊙\(2P−2m\)⊙Δm−1⊙2\+∑m=1P−12P−m​Ψ⊙\(2P−2m\)⊙τm,\\displaystyle=\\sum\_\{m=1\}^\{P\}2^\{P\-m\}\\Psi^\{\\odot\(2^\{P\}\-2^\{m\}\)\}\\odot\\Delta\_\{m\-1\}^\{\\odot 2\}\+\\sum\_\{m=1\}^\{P\-1\}2^\{P\-m\}\\Psi^\{\\odot\(2^\{P\}\-2^\{m\}\)\}\\odot\\tau\_\{m\},\(28\)which is merely a telescoping error propagation\.

Finally, letup=2​up−1\+up−12\+tpu\_\{p\}=2u\_\{p\-1\}\+u^\{2\}\_\{p\-1\}\+t\_\{p\}withu0=0u\_\{0\}=0, combining the linearity of the partial integration with inequality Eq\. \([28](https://arxiv.org/html/2608.17135#S5.E28)\) and applying the bounds

\|Λ⁡\(Ψ⊙\(2P−2m\)⊙τm,q\)\|≤‖Ψ⊙\(2P−2m\)‖2​‖τm‖2≤2​tm\\displaystyle\|\\Lambda\(\\Psi^\{\\odot\(2^\{P\}\-2^\{m\}\)\}\\odot\\tau\_\{m\},q\)\|\\leq\\\|\\Psi^\{\\odot\(2^\{P\}\-2^\{m\}\)\}\\\|\_\{2\}\\\|\\tau\_\{m\}\\\|\_\{2\}\\leq\\sqrt\{2\}t\_\{m\}\(29\)and

\|Λ⁡\(Ψ⊙\(2P−2m\)⊙Δm−12,q\)\|≤‖Ψ⊙\(2P−2m\)‖∞​‖Δm−1‖2≤um−12,\\displaystyle\|\\Lambda\(\\Psi^\{\\odot\(2^\{P\}\-2^\{m\}\)\}\\odot\\Delta\_\{m\-1\}^\{2\},q\)\|\\leq\\\|\\Psi^\{\\odot\(2^\{P\}\-2^\{m\}\)\}\\\|\_\{\\infty\}\\\|\\Delta\_\{m\-1\}\\\|\_\{2\}\\leq u\_\{m\-1\}^\{2\},\(30\)we obtain

\|Λ⁡\(ΔP,q\)\|≤∑m=1P−12P−m​2​tm\+∑m=1P2P−m​um−12\.\\displaystyle\\left\|\\Lambda\\\!\\left\(\\Delta\_\{P\},q\\right\)\\right\|\\leq\\sum\_\{m=1\}^\{P\-1\}2^\{P\-m\}\\sqrt\{2\}\\,t\_\{m\}\+\\sum\_\{m=1\}^\{P\}2^\{P\-m\}u\_\{m\-1\}^\{2\}\.\(31\)Here, the last truncation is omitted since fusing the final squaring into the BF contraction requires no compression \(tP=0t\_\{P\}=0\)\. Together with Eq\. \([26](https://arxiv.org/html/2608.17135#S5.E26)\), the optimality of the solution is certified provided

∑m=1P−12P−m​2​tm\+∑m=1P2P−m​um−12<1−2n​exp⁡\(−2PM\)\.\\displaystyle\\sum\_\{m=1\}^\{P\-1\}2^\{P\-m\}\\sqrt\{2\}\\,t\_\{m\}\+\\sum\_\{m=1\}^\{P\}2^\{P\-m\}u\_\{m\-1\}^\{2\}<1\-2^\{n\}\\exp\\left\(\-\\frac\{2^\{P\}\}\{M\}\\right\)\.\(32\)This is exactly the bound that the truncation errorstmt\_\{m\}must obey over the fullPPsteps for the optimality of the solution to be certified\.

## References

- McCulloch \(2007\)I\. P\. McCulloch, From density\-matrix renormalization group to matrix product states,[Journal of Statistical Mechanics: Theory and Experiment2007, P10014 \(2007\)](https://doi.org/10.1088/1742-5468/2007/10/P10014)\.
- Schollwöck \(2011\)U\. Schollwöck, The density\-matrix renormalization group in the age of matrix product states,[Annals of Physics326, 96 \(2011\)](https://www.sciencedirect.com/science/article/pii/S0003491610001752)\.
- Orús \(2014\)R\. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states,[Annals of Physics349, 117 \(2014\)](https://doi.org/10.1016/j.aop.2014.06.013),[arXiv:1306\.2164](https://arxiv.org/abs/1306.2164)\.
- Paeckel*et al\.*\(2019\)S\. Paeckel, T\. Köhler, A\. Swoboda, S\. R\. Manmana, U\. Schollwöck, and C\. Hubig, Time\-evolution methods for matrix\-product states,[Annals of Physics411, 167998 \(2019\)](https://doi.org/10.1016/j.aop.2019.167998),[arXiv:1901\.05824](https://arxiv.org/abs/1901.05824)\.
- Verstraete*et al\.*\(2023\)F\. Verstraete, T\. Nishino, U\. Schollwöck, M\. C\. Bañuls, G\. K\. Chan, and M\. E\. Stoudenmire, Density matrix renormalization group, 30 years on,[Nature Reviews Physics5, 273 \(2023\)](https://doi.org/10.1038/s42254-023-00572-5)\.
- Lubasch*et al\.*\(2018\)M\. Lubasch, P\. Moinier, and D\. Jaksch, Multigrid renormalization,[Journal of Computational Physics372, 587 \(2018\)](https://doi.org/10.1016/j.jcp.2018.06.065)\.
- Gourianov*et al\.*\(2022\)N\. Gourianov, M\. Lubasch, S\. Dolgov, Q\. Y\. van den Berg, H\. Babaee, P\. Givi, M\. Kiffner, and D\. Jaksch, A Quantum Inspired Approach to Exploit Turbulence Structures,[Nature Computational Science2, 30 \(2022\)](https://doi.org/10.1038/s43588-021-00181-1),[arXiv:2106\.05782 \[physics, physics:quant\-ph\]](https://arxiv.org/abs/2106.05782)\.
- Kiffner and Jaksch \(2023\)M\. Kiffner and D\. Jaksch, Tensor network reduced order models for wall\-bounded flows,[Physical Review Fluids8, 124101 \(2023\)](https://doi.org/10.1103/PhysRevFluids.8.124101)\.
- Gourianov*et al\.*\(2025\)N\. Gourianov, P\. Givi, D\. Jaksch, and S\. B\. Pope, Tensor networks enable the calculation of turbulence probability distributions,[Science Advances11, eads5990 \(2025\)](https://doi.org/10.1126/sciadv.ads5990)\.
- van Hülst*et al\.*\(2025\)N\.\-L\. van Hülst, P\. Siegl, P\. Over, S\. Bengoechea, T\. Hashizume, M\. G\. Cecile, T\. Rung, and D\. Jaksch,[Quantum\-inspired tensor\-network fractional\-step method for incompressible flow in curvilinear coordinates](https://arxiv.org/abs/2507.05222)\(2025\),[arXiv:2507\.05222 \[physics\.flu\-dyn\]](https://arxiv.org/abs/2507.05222)\.
- van Hülst*et al\.*\(2026\)N\.\-L\. van Hülst, M\. G\. Cecile, H\.\-Y\. Van, T\. Hashizume, E\. de Villiers, and D\. Jaksch,[Quantum\-Inspired Simulation of 2D Turbulent Rayleigh\-Bénard Convection](https://doi.org/10.48550/arXiv.2604.16179)\(2026\),[arXiv:2604\.16179 \[physics\]](https://arxiv.org/abs/2604.16179)\.
- Huynh*et al\.*\(2023\)L\. Huynh, J\. Hong, A\. Mian, H\. Suzuki, Y\. Wu, and S\. Camtepe,[Quantum\-Inspired Machine Learning: A Survey](https://doi.org/10.48550/arXiv.2308.11269)\(2023\),[arXiv:2308\.11269 \[cs\]](https://arxiv.org/abs/2308.11269)\.
- Han*et al\.*\(2018\)Z\.\-Y\. Han, J\. Wang, H\. Fan, L\. Wang, and P\. Zhang, Unsupervised Generative Modeling Using Matrix Product States,[Physical Review X8, 031012 \(2018\)](https://doi.org/10.1103/PhysRevX.8.031012)\.
- Chen*et al\.*\(2025\)Q\.\-C\. Chen, I\. Liu, J\.\-W\. Li, C\.\-M\. Chung,*et al\.*, Solving the gross\-pitaevskii equation with quantic tensor trains: Ground states and nonlinear dynamics, arXiv preprint arXiv:2507\.04279 \(2025\)\.
- Hölscher*et al\.*\(2025\)L\. Hölscher, P\. Rao, L\. Müller, J\. Klepsch, A\. Luckow, T\. Stollenwerk, and F\. K\. Wilhelm, Quantum\-inspired fluid simulation of two\-dimensional turbulence with GPU acceleration,[Physical Review Research7, 013112 \(2025\)](https://doi.org/10.1103/PhysRevResearch.7.013112)\.
- Niedermeier*et al\.*\(2026\)M\. Niedermeier, A\. Moulinas, T\. Louvet, J\. L\. Lado, and X\. Waintal, Solving the gross\-pitaevskii equation on multiple different scales using the quantics tensor train representation, Physical Review Research8, 023006 \(2026\)\.
- Michailidis*et al\.*\(2024\)A\. A\. Michailidis, C\. Fenton, and M\. Kiffner,[Tensor Train Multiplication](https://doi.org/10.48550/arXiv.2410.19747)\(2024\),[arXiv:2410\.19747 \[physics\]](https://arxiv.org/abs/2410.19747)\.
- Fernández*et al\.*\(2025\)Y\. N\. Fernández, M\. K\. Ritter, M\. Jeannin, J\.\-W\. Li, T\. Kloss, T\. Louvet, S\. Terasaki, O\. Parcollet, J\. von Delft, H\. Shinaoka, and X\. Waintal, Learning tensor networks with tensor cross interpolation: New algorithms and libraries,[SciPost Physics18, 104 \(2025\)](https://doi.org/10.21468/SciPostPhys.18.3.104),[arXiv:2407\.02454 \[physics\]](https://arxiv.org/abs/2407.02454)\.
- Danis*et al\.*\(2025\)M\. E\. Danis, D\. Truong, I\. Boureima, O\. Korobkin, K\. Ø\. Rasmussen, and B\. S\. Alexandrov, Tensor\-train WENO scheme for compressible flows,[Journal of Computational Physics529, 113891 \(2025\)](https://doi.org/10.1016/j.jcp.2025.113891)\.
- Meng*et al\.*\(2026\)Z\. Meng, Y\. Khoo, J\. Li, and E\. M\. Stoudenmire, Recursive sketched interpolation: Efficient hadamard products of tensor trains, arXiv preprint arXiv:2602\.17974 \(2026\)\.
- Needham \(1993\)T\. Needham, A Visual Explanation of Jensen’s Inequality,[The American Mathematical Monthly100, 768 \(1993\)](https://doi.org/10.2307/2324783),[2324783](https://arxiv.org/abs/2324783)\.
- Bullen \(2003\)P\. S\. Bullen,[*Handbook of Means and Their Inequalities*](https://doi.org/10.1007/978-94-017-0399-4)\(Springer Netherlands, Dordrecht, 2003\)\.
- Jarzynski \(1997\)C\. Jarzynski, Nonequilibrium equality for free energy differences,[Physical Review Letters78, 2690 \(1997\)](https://doi.org/10.1103/PhysRevLett.78.2690)\.
- Ruel and Ayres \(1999\)J\. J\. Ruel and M\. P\. Ayres, Jensen’s inequality predicts effects of environmental variation,[Trends in Ecology & Evolution14, 361 \(1999\)](https://doi.org/10.1016/S0169-5347(99)01664-X)\.
- Jumper*et al\.*\(2021\)J\. Jumper, R\. Evans, A\. Pritzel, T\. Green, M\. Figurnov, O\. Ronneberger, K\. Tunyasuvunakool, R\. Bates, A\. Žídek, A\. Potapenko, A\. Bridgland, C\. Meyer, S\. A\. A\. Kohl, A\. J\. Ballard, A\. Cowie, B\. Romera\-Paredes, S\. Nikolov, R\. Jain, J\. Adler, T\. Back, S\. Petersen, D\. Reiman, E\. Clancy, M\. Zielinski, M\. Steinegger, M\. Pacholska, T\. Berghammer, S\. Bodenstein, D\. Silver, O\. Vinyals, A\. W\. Senior, K\. Kavukcuoglu, P\. Kohli, and D\. Hassabis, Highly accurate protein structure prediction with AlphaFold,[Nature596, 583 \(2021\)](https://doi.org/10.1038/s41586-021-03819-2)\.
- Brouzet*et al\.*\(2021\)D\. Brouzet, M\. Talei, M\. J\. Brear, and B\. Cuenot, The impact of chemical modelling on turbulent premixed flame acoustics, J\. Fluid Mech\.915\(2021\)\.
- Dobrynin*et al\.*\(2024\)D\. Dobrynin, A\. Renaudineau, M\. Hizzani, D\. Strukov, M\. Mohseni, and J\. P\. Strachan, Energy landscapes of combinatorial optimization in Ising machines,[Physical Review E110, 045308 \(2024\)](https://doi.org/10.1103/PhysRevE.110.045308)\.
- Shen*et al\.*\(2025\)Z\.\-S\. Shen, F\. Pan, Y\. Wang, Y\.\-D\. Men, W\.\-B\. Xu, M\.\-H\. Yung, and P\. Zhang, Free\-energy machine for combinatorial optimization, Nature Computational Science , 1 \(2025\)\.
- Biamonte*et al\.*\(2011\)J\. D\. Biamonte, S\. R\. Clark, and D\. Jaksch, Categorical Tensor Network States,[AIP Advances1, 042172 \(2011\)](https://doi.org/10.1063/1.3672009)\.
- Stoudenmire and White \(2010\)E\. M\. Stoudenmire and S\. R\. White, Minimally Entangled Typical Thermal State Algorithms,[New Journal of Physics12, 055026 \(2010\)](https://doi.org/10.1088/1367-2630/12/5/055026),[arXiv:1002\.1305 \[cond\-mat\]](https://arxiv.org/abs/1002.1305)\.
- Bourchtein and Bourchtein \(2023\)A\. Bourchtein and L\. Bourchtein,*Elementary Functions*\(Springer, 2023\)\.
- Iacono and Boyd \(2015\)R\. Iacono and J\. P\. Boyd, The Kidder Equation:,[Studies in Applied Mathematics135, 63 \(2015\)](https://doi.org/10.1111/sapm.12073)\.
- Parand*et al\.*\(2017\)K\. Parand, P\. Mazaheri, M\. Delkhosh, and A\. Ghaderi, New numerical solutions for solving kidder equation by using the rational jacobi functions, SeMA Journal74, 569 \(2017\)\.
- Parand*et al\.*\(2018\)K\. Parand, S\. Latifi, M\. Delkhosh, and M\. M\. Moayeri, Generalized Lagrangian Jacobi Gauss collocation method for solving unsteady isothermal gas through a micro\-nano porous medium,[The European Physical Journal Plus133, 28 \(2018\)](https://doi.org/10.1140/epjp/i2018-11859-5)\.
- Cheon*et al\.*\(2020a\)J\. H\. Cheon, D\. Kim, and D\. Kim, Efficient Homomorphic Comparison Methods with Optimal Complexity, in[*Advances in Cryptology – ASIACRYPT 2020*](https://doi.org/10.1007/978-3-030-64834-3_8), edited by S\. Moriai and H\. Wang \(Springer International Publishing, Cham, 2020\) pp\. 221–256\.
- Cheon*et al\.*\(2020b\)J\. H\. Cheon, D\. Kim, and D\. Kim, Efficient homomorphic comparison methods with optimal complexity, in*Advances in Cryptology – ASIACRYPT 2020*, edited by S\. Moriai and H\. Wang \(Springer International Publishing, Cham, 2020\) pp\. 221–256\.
- Oseledets \(2011\)I\. V\. Oseledets, Tensor\-Train decomposition,[SIAM Journal on Scientific Computing33, 2295 \(2011\)](https://doi.org/10.1137/090752286)\.
- Connors \(1990\)K\. Connors,*Chemical Kinetics: The Study of Reaction Rates in Solution*\(VCH, 1990\)\.
- Greenhill*et al\.*\(2020\)S\. Greenhill, S\. Rana, S\. Gupta, P\. Vellanki, and S\. Venkatesh, Bayesian Optimization for Adaptive Experimental Design: A Review,[IEEE Access8, 13937 \(2020\)](https://doi.org/10.1109/ACCESS.2020.2966228)\.
- Turns \(2012\)S\. Turns,*An Introduction to Combustion: Concepts and Applications*, McGraw\-Hill Series in Mechanical Engineering \(McGraw\-Hill, 2012\)\.
- Arora and Barak \(2009\)S\. Arora and B\. Barak,*Computational Complexity: A Modern Approach*, 1st ed\. \(Cambridge University Press, USA, 2009\)\.
- max \(2016\)[Eleventh max\-sat evaluation](http://www.maxsat.udl.cat/16/benchmarks/index.html)\(2016\)\.
- Michailidis*et al\.*\(2025\)A\. A\. Michailidis, C\. Fenton, and M\. Kiffner, Element\-wise multiplication of tensor trains,[SIAM Journal on Scientific Computing47, B1158 \(2025\)](https://doi.org/10.1137/24M1714149)\.
- Sun and Chan \(2026\)J\. Sun and G\. K\. Chan, Stochastic tensor contraction for quantum chemistry, arXiv preprint arXiv:2602\.17158 \(2026\)\.
- Chen*et al\.*\(2023\)J\. Chen, E\. Stoudenmire, and S\. R\. White, Quantum Fourier Transform Has Small Entanglement,[PRX Quantum4, 040318 \(2023\)](https://doi.org/10.1103/PRXQuantum.4.040318)\.
- Hauck*et al\.*\(2025\)S\. H\. Hauck, M\. Kabel, and N\. R\. Gauger, Performance benchmarking of tensor trains for accelerated quantum\-inspired homogenization on tpu, gpu and cpu architectures \(2025\),[arXiv:2512\.07811](https://arxiv.org/abs/arXiv:2512.07811)\.
- Jiang*et al\.*\(2018\)S\. Jiang, K\. A\. Britt, A\. J\. McCaskey, T\. S\. Humble, and S\. Kais, Quantum Annealing for Prime Factorization,[Scientific Reports8, 17667 \(2018\)](https://doi.org/10.1038/s41598-018-36058-z)\.
- Zhang*et al\.*\(2026\)R\. Zhang, X\.\-Z\. Tao, Y\. Chen, J\. Zhang, L\. He, D\. Fang, L\. Yang, Y\. Sun, Q\. Zheng, X\. Shi,*et al\.*,[Qhap: Quantum\-inspired haplotype phasing](https://arxiv.org/abs/2603.25762)\(2026\),[arXiv:2603\.25762](https://arxiv.org/abs/2603.25762)\.
- Mansini and Speranza \(1999\)R\. Mansini and M\. G\. Speranza, Heuristic algorithms for the portfolio selection problem with minimum transaction lots,[European Journal of Operational Research114, 219 \(1999\)](https://doi.org/10.1016/S0377-2217(98)00252-5)\.
- Malz*et al\.*\(2024\)D\. Malz, G\. Styliaris, Z\.\-Y\. Wei, and J\. I\. Cirac, Preparation of matrix product states with log\-depth quantum circuits,[Physical Review Letters132, 040404 \(2024\)](https://doi.org/10.1103/PhysRevLett.132.040404),[arXiv:2307\.01696 \[quant\-ph\]](https://arxiv.org/abs/2307.01696)\.
- Termanova*et al\.*\(2024\)A\. Termanova, A\. Melnikov, E\. Mamenchikov, N\. Belokonev, S\. Dolgov, A\. Berezutskii, R\. Ellerbrock, C\. Mansell, and M\. Perelshtein, Tensor Quantum Programming, https://arxiv\.org/abs/2403\.13486v1 \(2024\)\.
- Siegl*et al\.*\(2025\)P\. Siegl, G\. S\. Reese, T\. Hashizume, N\.\-L\. van Hülst, and D\. Jaksch,[Tensor\-Programmable Quantum Circuits for Solving Differential Equations](https://doi.org/10.48550/arXiv.2502.04425)\(2025\),[arXiv:2502\.04425 \[quant\-ph\]](https://arxiv.org/abs/2502.04425)\.
- McCulloch and Gulácsi \(2002\)I\. P\. McCulloch and M\. Gulácsi, The non\-Abelian density matrix renormalization group algorithm,[Europhysics Letters57, 852 \(2002\)](https://doi.org/10.1209/epl/i2002-00393-0)\.
- Muller \(2006\)J\.\-M\. Muller,*Elementary functions: algorithms and implementation*\(Springer, 2006\)\.
- Tindall*et al\.*\(2025\)J\. Tindall, A\. Mello, M\. Fishman, M\. Stoudenmire, and D\. Sels,[Dynamics of disordered quantum systems with two\- and three\-dimensional tensor networks](https://arxiv.org/abs/2503.05693)\(2025\)\.

Supplemental Material for Iterative Tensor Network transformations for Element\-Wise Evaluation of Elementary and Filtering Functions

Xiao Wang1, Tomohiro Hashizume2,3, Pia Siegl3,4, Dieter Jaksch1,2,3

1Clarendon Laboratory, University of Oxford, Parks Road, Oxford OX1 3PU, United Kingdom 2The Hamburg Centre for Ultrafast Imaging, Hamburg 22761, Germany\. 3Institute for Quantum Physics, University of Hamburg, Luruper Chaussee 149, Hamburg 22761, Germany\. 4Institute of Software Methods for Product Virtualization, German Aerospace Center \(DLR\), Nöthnitzer Straße 46b, 01187 Dresden, Germany

\\justifying

## IApplication in nonlinear PDE

In this section, we demonstrate the application of the Newton\-Raphson\-based ITNT method in[subsection V\.3](https://arxiv.org/html/2608.17135#S5.SS3)to solve nonlinear PDEs in the TT representation, with particular focus on the Kidder equation, which models flow pressure in porous tubes\.

### I\.AUsing the Newton\-Raphson Method to realize the nonlinearity in Kidder equation

To solve the Kidder equation, we need to realize the element\-wise transformationw⁡\(x\)→\(1−κ​w​\(x\)\)−12w\(x\)\\to\\bigl\(1\-\\kappa w\(x\)\\bigr\)^\{\-\\frac\{1\}\{2\}\}in TT representation\. Following[subsection V\.3](https://arxiv.org/html/2608.17135#S5.SS3), we obtain the solution ofg⁡\(R⁡\(x\)\)=R​\(x\)−2−1−κ−w⁡\(x\)=0g\(R\(x\)\)=\\frac\{R\(x\)^\{\-2\}\-1\}\{\-\\kappa\}\-w\(x\)=0, via the iteration:

R\(k\+1\)​\(x\)=R\(k\)​\(x\)−R\(k\)​\(x\)−2−\(1−κ​w​\(x\)\)−2​R\(k\)​\(x\)−3,R^\{\(k\+1\)\}\(x\)=R^\{\(k\)\}\(x\)\-\\frac\{R^\{\(k\)\}\(x\)^\{\-2\}\-\\bigl\(1\-\\kappa w\(x\)\\bigr\)\}\{\-2R^\{\(k\)\}\(x\)^\{\-3\}\},Simplifying this expression, we obtain the following iteration in TT representation:

R\(k\+1\)=12​R\(k\)⊙\(3​𝕀−\(1−κ​w\)⊙R\(k\)⊙R\(k\)\)\.R^\{\(k\+1\)\}=\\frac\{1\}\{2\}R^\{\(k\)\}\\odot\\left\(3\\mathbb\{I\}\-\\bigl\(1\-\\kappa w\\bigr\)\\odot R^\{\(k\)\}\\odot R^\{\(k\)\}\\right\)\.\(S1\)Here,R\(k\)R^\{\(k\)\}andwware TT encoding the functionR\(k\)R^\{\(k\)\}andw⁡\(x\)w\(x\)on the equally discretizedxx\-grid\. The term𝕀\\mathbb\{I\}is aχ=1\\chi=1unit TT where all the element of its local tensors are11as explained in the main text\. For\|κ\|<1\|\\kappa\|<1and\|w⁡\(x\)\|<1\|w\(x\)\|<1, the iteration converges when starting from a unit\-element TTR\(0\)=𝕀R^\{\(0\)\}=\\mathbb\{I\}which representsR\(0\)​\(x\)=1R^\{\(0\)\}\(x\)=1\. Then, the iteration converges rapidly due to the quadratic convergence of the Newton\-Raphson method, such that the final TT,limk→∞\|R\(k\)⟩\\lim\_\{k\\to\\infty\}\\ket\{R^\{\(k\)\}\}, accurately represents\(1−κ​w​\(x\)\)−12\(1\-\\kappa w\(x\)\)^\{\-\\frac\{1\}\{2\}\}\.

### I\.BConducting the self\-consistent Picard iteration

The Kidder equation is given by:

w′′​\(x\)\+2​x​w′​\(x\)​\(1−κ​w​\(x\)\)−12=0,w^\{\\prime\\prime\}\(x\)\\;\+\\;2\\,x\\,w^\{\\prime\}\(x\)\\,\\bigl\(1\-\\kappa\\,w\(x\)\\bigr\)^\{\-\\tfrac\{1\}\{2\}\}\\;=\\;0,subject to the boundary conditions

w⁡\(0\)=1andlimx→∞w⁡\(x\)=0\.w\(0\)\\;=\\;1\\quad\\text\{and\}\\quad\\lim\_\{x\\to\\infty\}w\(x\)\\;=\\;0\.To solve it numerically, we discretize the interval\[0,xmax\]\[0,x\_\{\\max\}\]into2n2^\{n\}grid pointsxix\_\{i\}, with inter\-grid point spacingΔ​x=xmax/2n\\Delta x=x\_\{\\max\}/2^\{n\}\. Here,nnrepresents the number of tensors required to store the discretized field in the TT format\. For the Kidder equation withκ=0\.5\\kappa=0\.5, it is sufficient to setxmax=16x\_\{\\max\}=16to approximatexmax≈∞x\_\{\\max\}\\approx\\infty\([33](https://arxiv.org/html/2608.17135#bib.bib33)\), leading to the grid resolution ofΔ​x=2−12\\Delta x=2^\{\-12\}\. We label the approximate solution bywi≈w⁡\(xi\)w\_\{i\}\\approx w\(x\_\{i\}\)and use finite differences to approximatew′​\(x\)w^\{\\prime\}\(x\)andw′′​\(x\)w^\{\\prime\\prime\}\(x\)\. For the interior pointsi=1,…,2ni=1,\\dots,2^\{n\}, the discrete system typically reads as follows:

Fi​\(𝐰\)=wi\+1−2​wi\+wi−1\(Δ​x\)2\+2​xi​wi\+1−wi−12​Δ​x​\(1−κ​wi\)−12=0\.\\begin\{split\}F\_\{i\}\(\\mathbf\{w\}\)&\\;=\\;\\frac\{w\_\{i\+1\}\-2\\,w\_\{i\}\+w\_\{i\-1\}\}\{\(\\Delta x\)^\{2\}\}\\\\ &~~\\;\+\\;2\\,x\_\{i\}\\,\\frac\{w\_\{i\+1\}\-w\_\{i\-1\}\}\{2\\,\\Delta x\}\\,\\bigl\(1\-\\kappa\\,w\_\{i\}\\bigr\)^\{\-\\tfrac\{1\}\{2\}\}\\;=\\;0\.\\end\{split\}\(S2\)We keepw0=1w\_\{0\}=1andw2n=0w\_\{2^\{n\}\}=0fixed to enforce the boundary conditions\.

A Picard fixed\-point iteration solves this nonlinear system by rewriting each equation to isolatewiw\_\{i\}on one side\. Concretely, one arranges:

wi=Gi​\(wi−1,wi,wi\+1\),w\_\{i\}\\;=\\;G\_\{i\}\\bigl\(w\_\{i\-1\},\\,w\_\{i\},\\,w\_\{i\+1\}\\bigr\),so that each interior variablewiw\_\{i\}is expressed in terms of the old iteration values\. The explicit update formula, after algebraic rearrangement, reads

Gi​\(𝐰\(k\)\)=wi\+1\(k\)\+wi−1\(k\)2\+xi​Δ​x2​\[wi\+1\(k\)−wi−1\(k\)\]​\(1−κ​wi\(k\)\)−12,\\begin\{split\}G\_\{i\}\\bigl\(\\mathbf\{w\}^\{\(k\)\}\\bigr\)&=\\frac\{w\_\{i\+1\}^\{\(k\)\}\+w\_\{i\-1\}^\{\(k\)\}\}\{2\}\\\\ &\+\\frac\{x\_\{i\}\\,\\Delta x\}\{2\}\\,\\Bigl\[w\_\{i\+1\}^\{\(k\)\}\-w\_\{i\-1\}^\{\(k\)\}\\Bigr\]\\bigl\(1\-\\kappa\\,w\_\{i\}^\{\(k\)\}\\bigr\)^\{\-\\tfrac\{1\}\{2\}\},\\end\{split\}\(S3\)wherewi−1\(k\),wi\+1\(k\)w\_\{i\-1\}^\{\(k\)\},\\,w\_\{i\+1\}^\{\(k\)\}, andwi\(k\)w\_\{i\}^\{\(k\)\}come from the previous iteration\. In TT representation, this Picard iteration becomes

G⁡\(w\(k\)\)=12​\(L^\+R^\)⋅w\(k\)\+Δ​x2X⊙\(1−κw\(k\)\)−12⊙\(L^−R^\)⋅w\(k\)\.\.\\begin\{split\}&G\\bigl\(w^\{\(k\)\}\\bigr\)=\\frac\{1\}\{2\}\\left\(\\hat\{L\}\+\\hat\{R\}\\right\)\\cdot w^\{\(k\)\}\\\\ &\+\\frac\{\\Delta x\}\{2\}X\\odot\\bigl\(1\-\\kappa\\,w^\{\(k\)\}\\bigr\)^\{\-\\tfrac\{1\}\{2\}\}\\odot\\left\(\\hat\{L\}\-\\hat\{R\}\\right\)\\cdot w^\{\(k\)\}\.\\end\{split\}\.\(S4\)Here,XXis a TT that encodes the coordinates of the discretized grid \(X\[i\]=i​Δ​xX^\{\[i\]\}=i\\Delta x\)k,L^\\hat\{L\}andR^\\hat\{R\}are linear operators that shift the index by one, i\.e\.Ψ′=L^​Ψ\\Psi^\{\\prime\}=\\hat\{L\}\\Psis\.t\.Ψ′\[j\]=Ψ\[j\+1\]\\Psi^\{\\prime\[j\]\}=\\Psi^\{\[j\+1\]\}andΨ′′=R^​Ψ\\Psi^\{\\prime\\prime\}=\\hat\{R\}\\Psis\.t\.Ψ′′\[j\]=Ψ\[j−1\]\\Psi^\{\\prime\\prime\[j\]\}=\\Psi^\{\[j\-1\]\}\([6](https://arxiv.org/html/2608.17135#bib.bib6)\)\. The TT representation of the term\(1−κ​wi\(k\)\)−12\(1\-\\kappa w\_\{i\}^\{\(k\)\}\)^\{\-\\tfrac\{1\}\{2\}\}, denoted by\(1−κ​wOPEN\(k\)\)−12CLOSE\\bigl\(1\-\\kappa\\,w^\{\(k\)\\bigr\)^\{\-\\tfrac\{1\}\{2\}\}\}, is given by ITNTs in Eq\. \([S1](https://arxiv.org/html/2608.17135#S1.E1)\)\. We start from an initial guessw\(0\)w^\{\(0\)\}and perform the iteration

w\(k\+1\)=G\(w\(k\)\),k=0,1,2,…w^\{\(k\+1\)\}=G\\bigl\(w^\{\(k\)\}\\bigr\),\\quad k=0,1,2,\\dotsuntil successive solutions are sufficiently close\. This fixed\-point iteration converges to the unique TTwwsatisfying the discretized Kidder equation and all of its boundary conditions\.

### I\.CAnderson Acceleration

While a basic Picard iteration is often sufficient for moderate grids or small parameter regimes, it can converge slowly or encounter numerical difficulties for larger or stiffer systems\. Anderson Acceleration \(also known as Anderson mixing\) enhances the fixed\-point iteration by forming and solving a small least\-squares problem at each step, combining information from multiple past iterates and residuals to produce a more effective update\.

Concretely, suppose we have at iterationnn:

w\(k\+1\)=G⁡\(w\(k\)\)\.w^\{\(k\+1\)\}=G\\bigl\(w^\{\(k\)\}\\bigr\)\.Define the residual as

f\(k\)=G⁡\(w\(k\)\)−w\(k\)\.f^\{\(k\)\}=G\\bigl\(w^\{\(k\)\}\\bigr\)\-w^\{\(k\)\}\.Rather than simply setting

w\(k\+1\)=w\(k\)\+f\(k\),w^\{\(k\+1\)\}=w^\{\(k\)\}\+f^\{\(k\)\},we store several past pairs\{w\(k−i\),f\(k−i\)\}\\\{w^\{\(k\-i\)\},f^\{\(k\-i\)\}\\\}and solve a small linear or least\-squares system to find combination coefficientsαi\\alpha\_\{i\}that minimize the new residual\. A commonly used variant enforces∑iαi=1\\sum\_\{i\}\\alpha\_\{i\}=1, preserving affine invariance\. Then the next iterate can be written, for example, as

w\(k\+1\)=∑iαi​\(w\(k−i\)\+f\(k−i\)\),with​∑iαi=1\.w^\{\(k\+1\)\}\\;=\\;\\sum\_\{i\}\\alpha\_\{i\}\\,\\bigl\(w^\{\(k\-i\)\}\+f^\{\(k\-i\)\}\\bigr\),~\\text\{with\}~\\sum\_\{i\}\\alpha\_\{i\}=1\.The small system to find a set of scalars\{αi\}\\\{\\alpha\_\{i\}\\\}by computing dot products among the stored residuals\. Once thoseαi\\alpha\_\{i\}are determined, one combines the corresponding old solutions,w\(k−i\)\+f\(k−i\)=G⁡\(w\(k−i\)\)w^\{\(k\-i\)\}\+f^\{\(k\-i\)\}=G\(w^\{\(k\-i\)\}\), to producew\(k\+1\)w^\{\(k\+1\)\}\. This linear combination step can drastically reduce the number of iterations required for convergence, especially in stiff or nonlinear regimes\.

In practice, we limit the memory size \(e\.g\., storingmmprevious iterates\) to keep the least\-squares problem small\. In a TT\-based Kidder solver, the main computational effort lies in evaluating inner products of residuals𝐟\(k\)\\mathbf\{f\}^\{\(k\)\}stored in tensor form and carefully truncating the TT rank to control memory usage\. Anderson Acceleration thus provides a flexible, robust approach to speeding up Picard iteration for the discretized Kidder equation, particularly when direct methods \(like Newton’s method\) are more difficult to implement or too expensive at a large scale\.

### I\.DResults

We set the SVD tolerance to10−1210^\{\-12\}withχ=400\\chi=400, and conduct the Picard iteration for10510^\{5\}times with an Anderson acceleration with a memory ofm=10m=10\. After this, we prolong the TT to a finer grid \(see Ref\.\([6](https://arxiv.org/html/2608.17135#bib.bib6)\)\) on which we re\-run the above iteration\. We start from a grid finesse ofNres=4N\_\{\\text\{res\}\}=4and end atNres=17N\_\{\\text\{res\}\}=17, the final result provides an initial slope of \-1\.1917912520439131, which shows an error of∼6×10−7\\sim 6\\times 10^\{\-7\}compared to the standard result given in Ref\.\([33](https://arxiv.org/html/2608.17135#bib.bib33)\)\.

Here, we comment on the limitations of Picard iterations\. To reach a high accuracy in the result \(e\.g\., to find the initial slope to 10 digits\), we need to improve the SVD tolerance to∼10−14\\sim 10^\{\-14\}\. In this case, each iteration step becomes slower as the bond dimension increases\. Besides, as the grid finesseNresN\_\{\\text\{res\}\}reaches∼18\\sim 18, even the Anderson acceleration becomes slow to reach convergence; thus, we need a more efficient scheme to accelerate the Picard iteration, such as Newton’s method\. Neverthless, our ITNT enables the step of finding\(1−κ​w​\(x\)\)−0\.5\(1\-\\kappa w\(x\)\)^\{\-0\.5\}via the Newton\-Raphson method in the TT representation, and this approach is highly effective\.

## IITruncation error bounding

For a positive fieldΨ\\Psi, with\(⋯\)SVD\(\\cdots\)\_\{\\mathrm\{SVD\}\}denoting the field after SVD truncation, the maximum and the sub\-maximum ofΨ⊙2p\\Psi^\{\\odot 2^\{p\}\}retain their order provided

ε⁡\(Ψ⊙2p,\(Ψ⊙2p\)SVD\)≤\(max⁡\{\(Ψ⊙2p\)\[j\]\}−maxsub⁡\{\(Ψ⊙2p\)\[j\]\}\)/2\\displaystyle\\varepsilon\(\\Psi^\{\\odot\{2^\{p\}\}\},\(\\Psi^\{\\odot\{2^\{p\}\}\}\)\_\{\\mathrm\{SVD\}\}\)\\leq\(\\max\\\{\(\\Psi^\{\\odot 2^\{p\}\}\)^\{\[j\]\}\\\}\-\\max\_\{\\mathrm\{sub\}\}\\\{\(\\Psi^\{\\odot 2^\{p\}\}\)^\{\[j\]\}\\\}\)/\\sqrt\{2\}\(S5\)as presented in Eq\. \([2](https://arxiv.org/html/2608.17135#S2.E2)\)\. In this appendix, we derive the right\-hand side of this inequality by showing that1/21/\\sqrt\{2\}is the largest factor for which order preservation is guaranteed\.

LetΦ=Ψ⊙2p\\Phi=\\Psi^\{\\odot 2^\{p\}\}andΦ′=\(Ψ⊙2p\)SVD\\Phi^\{\\prime\}=\(\\Psi^\{\\odot 2^\{p\}\}\)\_\{\\mathrm\{SVD\}\}for the exact and truncated arrays\. For simplicity we writeε=‖Φ−Φ′‖2\\varepsilon=\\left\\\|\\Phi\-\\Phi^\{\\prime\}\\right\\\|\_\{2\}for the truncation error andΔ=max⁡\{Φ\}−maxsub⁡\{Φ\}\>0\\Delta=\\max\\\{\\Phi\\\}\-\\max\_\{\\mathrm\{sub\}\}\\\{\\Phi\\\}\>0for the gap between maximum and the second maximal value of the untruncated field\.

Truncation displaces these two valus by

δ1=max⁡\{Φ\}−max⁡\{Φ′\},δ2=maxsub⁡\{Φ′\}−maxsub⁡\{Φ\}\\displaystyle\\delta\_\{1\}=\\max\\\{\\Phi\\\}\-\\max\\\{\\Phi^\{\\prime\}\\\},\\quad\\delta\_\{2\}=\\max\_\{\\mathrm\{sub\}\}\\\{\\Phi^\{\\prime\}\\\}\-\\max\_\{\\mathrm\{sub\}\}\\\{\\Phi\\\}\(S6\)which measure the depression of the maximum \(δ1\\delta\_\{1\}\) and the elevation of the second maximum \(δ2\\delta\_\{2\}\)\.

In order for the order of the maximum and second maximum to be conserved,δ1\\delta\_\{1\}andδ2\\delta\_\{2\}must meetδ1\+δ2≤Δ\\delta\_\{1\}\+\\delta\_\{2\}\\leq\\Delta, while the truncation error irself constraints them byδ12\+δ22≤ε2\\delta\_\{1\}^\{2\}\+\\delta\_\{2\}^\{2\}\\leq\\varepsilon^\{2\}\. Invoking Cauchy\-Schwarz inequality, we obtain

δ1\+δ2≤2​δ12\+δ22≤2​ε≤Δ\\displaystyle\\delta\_\{1\}\+\\delta\_\{2\}\\leq\\sqrt\{2\}\\sqrt\{\\delta\_\{1\}^\{2\}\+\\delta\_\{2\}^\{2\}\}\\leq\\sqrt\{2\}\\varepsilon\\leq\\Delta\(S7\)where the first inequality is saturated whenδ1=δ2\\delta\_\{1\}=\\delta\_\{2\}, and the second whenδ1=δ2=ε/2\\delta\_\{1\}=\\delta\_\{2\}=\\varepsilon/\\sqrt\{2\}\. With this, Eq\. \([2](https://arxiv.org/html/2608.17135#S2.E2)\) is shown\.

## IIIApplication of ITNT as Sign transformation \- Accuracy and Scaling

To test the accuracy and the convergence properties of the sign transform, we consider the one\-dimensional continuous function

![Refer to caption](https://arxiv.org/html/2608.17135v1/figures/N14_Chi20_Step20_roots7.png)
![Refer to caption](https://arxiv.org/html/2608.17135v1/figures/N14_Chi25_Step20_roots7.png)
![Refer to caption](https://arxiv.org/html/2608.17135v1/figures/N18_Chi25_Step20_roots7.png)
![Refer to caption](https://arxiv.org/html/2608.17135v1/figures/N18_Chi75_Step20_roots23.png)
f1​\(x\)f\_\{1\}\(x\)ITNT​\(f1​\(x\)\)\\mathrm\{ITNT\}\(f\_\{1\}\(x\)\)

Figure S1:Result of the Sign\-function transformation via ITNT of the function defined in Eq\.[S8](https://arxiv.org/html/2608.17135#S3.E8)for different number of rootsNrootN\_\{\\text\{root\}\}and different qubit numbersnn\. For all examples we considerNiter=20N\_\{\\text\{iter\}\}=20and an SVD tolerance of10−810^\{\-8\}\. Sign transformation of the function with the following properties: \(a\)n=14n=14,Nroot=7N\_\{\\text\{root\}\}=7,χmax=20\\chi\_\{\\text\{max\}\}=20\. \(b\)n=14n=14,Nroot=7N\_\{\\text\{root\}\}=7,χmax=25\\chi\_\{\\text\{max\}\}=25\. \(c\)n=18n=18,Nroot=7N\_\{\\text\{root\}\}=7,χmax=25\\chi\_\{\\text\{max\}\}=25\. \(d\)n=18n=18,Nroot=23N\_\{\\text\{root\}\}=23,χmax=75\\chi\_\{\\text\{max\}\}=75\. ForNroot=23N\_\{\\text\{root\}\}=23, a noiseless sign transformation requiresχmax=75\\chi\_\{\\text\{max\}\}=75\.![Refer to caption](https://arxiv.org/html/2608.17135v1/figures/Figure_sign_transform_and_scaling.png)Figure S2:Scaling of the element\-wise evaluation of sign function via ITNT\(a\) Element\-wise evaluation of sign functiony=Sign​\{fω​\(x\)\}y=\\mathrm\{Sign\}\\\{f\_\{\\omega\}\(x\)\\\}\(red\) using ITNT of a functiony=fω​\(x\)y=f\_\{\\omega\}\(x\)\(black\)ω=1\.0\\omega=1\.0\. \(b\) Required bond dimension in dependence on the number of tensorsnnforω=0\.5\\omega=0\.5\. Considered are the maximal intermediate bond\-dimensionχpeak\\chi\_\{\\mathrm\{peak\}\}\(crosses\) used by the algorithm if the bond dimension is not restricted and the SVD tolerance is fixed at10−810^\{\-8\}\. Furthermore, we consider the final bond dimensionχfinal\\chi\_\{\\mathrm\{final\}\}\(circles\) required to representyyaccurate up to the order of the SVD tolerance of10−810^\{\-8\}with respect to theL2L^\{2\}\-error\. Finally, there is the cutoff bond dimensionχmin\\chi\_\{\\mathrm\{min\}\}\(squares\) which defines a lower cutoff bound, required to reach the result with ITNT up to this accuracy while reaching the bond dimensionχfinal\\chi\_\{\\mathrm\{final\}\}within200200iteration steps\. \(c\) Required bond dimensions for different numbers of rootsNrootsN\_\{\\mathrm\{roots\}\}by usingω=\[0\.5,1\.0,1\.5,2\.0\]\\omega=\[0\.5,1\.0,1\.5,2\.0\]corresponding toNroots=\[7,13,19,25\]N\_\{\\mathrm\{roots\}\}=\[7,13,19,25\]andn=18n=18to converge within a maximum of200200iteration steps\. The bond dimensions are defined as for \(b\)\. \(d\) The number of iteration stepsNiterN\_\{\\mathrm\{iter\}\}needed to reach the final compressed result with respect to the chosen maximal cutoff bond dimensionχmax\\chi\_\{\\rm\{max\}\}forω=0\.5\\omega=0\.5forn=14,18,22n=14,18,22\(light to dark\)\. All simulations were implemented using the zip\-up algorithm for point\-wise multiplication\.![Refer to caption](https://arxiv.org/html/2608.17135v1/figures/MaxBondDim_vs_IterStep_small_ChiMax.png)Figure S3:The bond dimension of the sign\-transformed TT afterNiterN\_\{\\text\{iter\}\}steps of iteration\. The qubit number isn=18n=18with root numberNroot=7N\_\{\\text\{root\}\}=7, andχmax\\chi\_\{\\text\{max\}\}denotes the maximum bond dimension allowed during TT operations \(e\.g\., element\-wise products\), with the SVD tolerance fixed at10−810^\{\-8\}\. Forχmax=28\\chi\_\{\\text\{max\}\}=28, the convergence is reached atNiter≈64N\_\{\\text\{iter\}\}\\approx 64\. Forχmax≳35\\chi\_\{\\text\{max\}\}\\gtrsim 35, the required iteration count reduces toNiter≈34N\_\{\\text\{iter\}\}\\approx 34\.fω​\(x\)=cos⁡\(ω​x\)−0\.000035​\(x−2\)3\+0\.001​x2−0\.02f\_\{\\omega\}\(x\)=\\cos\(\\omega x\)\-0\.000035\(x\-2\)^\{3\}\+0\.001x^\{2\}\-0\.02\(S8\)encoded acrossnntensors, wherex∈\[0,64\)x\\in\[0,64\)and the frequencyω\\omegadetermines the number of roots on this interval\. Varying the number of roots, we study the impact of the chosen cutoff bond dimensionχmax\\chi\_\{\\rm\{max\}\}on the convergence of the algorithm\. The result of the ITNT\-based sign\-transformation offω​\(x\)f\_\{\\omega\}\(x\)is depicted in Fig\.[S1](https://arxiv.org/html/2608.17135#S3.F1)for different values ofω\\omega, and hence different number of rootsNrootN\_\{\\rm\{root\}\}, and with a changing bond\-dimension\. For the simulation shown in Fig\.[S1](https://arxiv.org/html/2608.17135#S3.F1)\(a\), the chosen cutoff bond dimension is visibly not large enough to ensure the full convergence of the simulation, leading to additional oscillations around the correct result\. For the simulation results shown in Fig\.[S1](https://arxiv.org/html/2608.17135#S3.F1)\(b\)\-\(d\), the chosenχmax\\chi\_\{\\rm\{max\}\}is sufficient for convergence\. Hereby, the necessaryχmax\\chi\_\{\\rm\{max\}\}seems to depend less on the number of grid points than on the number of roots\. A more detailed analysis of these dependencies is shown in Fig\.[S2](https://arxiv.org/html/2608.17135#S3.F2)\(a\) and \(b\), where we study the required bond dimension over the number of tensors and number of roots\. Interestingly, the bond dimension required to represent the final resultχfinal\\chi\_\{\\rm\{final\}\}up to an accuracy of order10−810^\{\-8\}, saturates quickly with the number of tensors while showing a linear dependence onNrootN\_\{\\rm\{root\}\}\. Furthermore, we consider the minimum cutoff bond dimensionχmin\\chi\_\{\\rm\{min\}\}required during ITNT to allow convergence within200200steps\. Thisχmin\\chi\_\{\\rm\{min\}\}is significantly lower than the bond dimensionχpeak\\chi\_\{\\rm\{peak\}\}, the bond dimension used by the algorithm when no cutoff bond dimension is set, and only the SVD tolerance of10−810^\{\-8\}is defined\. This indicates that choosing a cutoff bond dimension can be beneficial to bound the resource requirements\. However, it can be beneficial to choose it slightly larger thanχmin\\chi\_\{\\rm\{min\}\}to allow for a quicker convergence \(c\.f\.[S2](https://arxiv.org/html/2608.17135#S3.F2)\(c\)\)\. In our simulations, we observe that selecting a sufficiently large bond dimensionχmax≳5​Nroot\\chi\_\{\\text\{max\}\}\\gtrsim 5N\_\{\\text\{root\}\}is crucial to minimize the convergent iteration step,NiterN\_\{\\text\{iter\}\}\. Furthermore, we find that when the maximum bond dimension is large enough \(i\.e\.,χmax≳5​Nroot\\chi\_\{\\text\{max\}\}\\gtrsim 5N\_\{\\text\{root\}\}\), the sign transformation ITNT fully converges afterNiter≈2​NN\_\{\\text\{iter\}\}\\approx 2Nsteps of iteration\. For such a large enoughχmax\\chi\_\{\\text\{max\}\}, after the ITNT reaches convergence, the bond dimension of the TT shrinks back roughly to the root numberNrootN\_\{\\text\{root\}\}\. We can truncate the bond dimension toχmax\\chi\_\{\\text\{max\}\}up toχmax∼4​Nroot\\chi\_\{\\text\{max\}\}\\sim 4N\_\{\\text\{root\}\}\(c\.f\. Fig\.[S3](https://arxiv.org/html/2608.17135#S3.F3)\)\. In this case, the convergence can still be reached, although it is reached at a larger iteration stepNiterN\_\{\\text\{iter\}\}\. As we truncate to an even smallerχmax\\chi\_\{\\text\{max\}\}, the iteration converges very slowly, and the whole iteration may eventually diverge due to numerical errors\.

The runtime scalability with qubit numberNNis evident: as we increaseNN, the bond dimension and runtime remain manageable\. The runtimes required to perform the sign transformation with ITNT are summarized in Tab\.[1](https://arxiv.org/html/2608.17135#S3.T1)and compared with a standard sign transformation using full\-grid point\-wise evaluation in Python/NumPy\. Both computations were performed on the same laptop with an Intel Core Ultra 9 275HX CPU\.

nnNrootN\_\{\\text\{root\}\}χmax\\chi\_\{\\text\{max\}\}time \(sec\) of ITNT\-based sign tranformstime \(sec\) of standard sign transform\.2274531\.6221910071\.62251250651\.6247453\.56\.524191008\.56\.52451250846\.52674552826191001128265125011028Table 1:Comparison of runtime performance across different qubit numbersnn, root numbersNrootN\_\{\\rm\{root\}\}, and chosen cutoff bond dimensionsχmax\\chi\_\{\\rm\{max\}\}\.A larger SVD tolerance can further boost the convergence speed, so we need to choose this tolerance properly\. In our above simulations, we find that increasing the tolerance from10−810^\{\-8\}to10−610^\{\-6\}will only slightly increase the error of the ITNT result \(i\.e\., theL2L^\{2\}\-norm of the difference compared with the element\-wise transformed array\), from∼10−8\\sim 10^\{\-8\}to∼10−6\\sim 10^\{\-6\}\.

At the core of the sign\-transformation lays the the homomorphic comparison iteration,F⁡\(f\)=0\.5​f​\(3−f2\)F\(f\)=0\.5f\(3\-f^\{2\}\)\. In the following we will consider its iteration speed for different distances to the root\.

#### Nearf=±1f=\\pm 1: Double\-exponential Convergence

The stable fixed points ofFFare atf∗=±1f^\{\*\}=\\pm 1\. Nearf=1f=1, letϵ\(k\)=1−f\(k\)\\epsilon^\{\(k\)\}=1\-f^\{\(k\)\}represent the small deviation from the fixed point\. The recurrence is approximatelyϵ\(k\+1\)≈1\.5​\(ϵ\(k\)\)2\\epsilon^\{\(k\+1\)\}\\approx 1\.5\\,\(\\epsilon^\{\(k\)\}\)^\{2\}\. This leads to double\-exponential convergenceϵ\(k\)=\(1\.5​ϵ\(0\)\)2k/1\.5\\epsilon^\{\(k\)\}=\(1\.5\\,\\epsilon^\{\(0\)\}\)^\{2^\{k\}\}/\{1\.5\}\. If1\.5​ϵ\(0\)<11\.5\\epsilon^\{\(0\)\}<1,ϵ\(k\)\\epsilon^\{\(k\)\}decays extremely rapidly, ensuring that the ITNT quickly approachesf∗=1f^\{\*\}=1\. A similar argument applies nearf=−1f=\-1\.

#### Nearf=0f=0: Exponential Departure and Resolution Threshold

In contrast to the stable fixed points atf=±1f=\\pm 1, the pointf=0f=0is an unstable fixed point\. To see this, we linearizeFFaroundf=0f=0:

F⁡\(f\)=0\.5​f​\(3−f2\)≈1\.5​ffor small​f\.F\(f\)=0\.5f\(3\-f^\{2\}\)\\approx 1\.5f\\quad\\text\{for small \}f\.Thus, any small nonzero value offfnear zero grows by approximately 1\.5 per iteration,f\(k\)≈f\(0\)×\(1\.5\)kf^\{\(k\)\}\\approx f^\{\(0\)\}\\times\(1\.5\)^\{k\}\. This results in exponential growth away from zero, repelling fromf=0f=0\. Suppose that afterSSsteps of iteration, the iterative process realizes a sign\-like function, i\.e\.,f\(S\)≈1f^\{\(S\)\}\\approx 1, we can estimate the “resolution threshold” for the initial value,f\(0\)≈\(1\.5\)−Sf^\{\(0\)\}\\approx\(1\.5\)^\{\-S\}\.

This threshold,\(1\.5\)−S\(1\.5\)^\{\-S\}, sets up a resolution limit in approximating the sign function\. AfterSSsteps, any initial valuesf\(0\)f^\{\(0\)\}\(and hence portions of the functionf⁡\(x\)f\(x\)\) larger than\(1\.5\)−S\(1\.5\)^\{\-S\}are effectively sign\-transformed and pushed toward11, while those smaller than\(1\.5\)−S\(1\.5\)^\{\-S\}remain evidently far from 1\. In other words, after sufficientlySSsteps of iteration, the iterative map cannot distinguish initial value differences larger than\(1\.5\)−S\(1\.5\)^\{\-S\}from 1, thus achieving a sign\-like output up to that resolution level\.

Finally, we provide details of applying the ITNT\-based sign function to locate the roots of smooth 1D functions\. For a smooth functionf⁡\(x\)f\(x\), the intermediate value theorem implies that whenever there is a change in the sign off⁡\(x\)f\(x\), there exists at least one root of a function\. When a curve is represented as an TT, roots can be found by first applying the sign\-transformation ITNT, then, we apply the differential differential linear operatorR^−L^\\hat\{R\}\-\\hat\{L\}for the linear operatorsR^\\hat\{R\}andL^\\hat\{L\}introduced in Eq\. \([S4](https://arxiv.org/html/2608.17135#S1.E4)\) to obtain the TT that marks the root positions as delta peaks\. From the resulting multi\-delta\-peak TT, we use the deflation method in Sec\.[V\.5](https://arxiv.org/html/2608.17135#S5.SS5)to iteratively locate each peak position in the TT representation\. This procedure outputs the set of root positions on the discretized grid\.

## IVBenchmark: performance of extremum finding on discrete problem

Here, we focus on a classical Ising model on a square lattice of size12×1212\\times 12with nearest\-neighbor interactions and on\-site terms, where all parameters are drawn randomly from\(−0\.5,0\.5\)\(\-0\.5,0\.5\)\(c\.f\. Methods[V\.6](https://arxiv.org/html/2608.17135#S5.SS6)for details of the model\)\. The energy landscape of the classical Ising model can be mapped onto energy\-landscape TTEEwith the bond dimension that is proportional to the number of interaction terms in the Hamiltonian\. For this example, with the open boundary condition, we foundχ=409\\chi=409achieves the exact encoding of the energy landscape\. To identify the maximum energy configuration, we fix the Hamiltonian and perform self\-multiplicationp=8p=8times, then perform 800 independent deflation runs using 100 different stochastic seeds\.

In Fig\.[S4](https://arxiv.org/html/2608.17135#S4.F4)\(a\), we present the statistics of the located extremum values with respect to the iteration steps for all 100 independent runs with the same Hamiltonian but different random site sequences for binary fixation\. The percentage on the right of each data point represents the fraction of runs that identified the value as the maximum energy configuration within those steps\. While the subsequent ITNT steps identifies and converges to the true optimum, as increasingly many trajectories find higher energy configurations over the iterations\. In practice, each deflation step identifies a pronounced peak close to the optimal solution; after 800 stochastic deflation steps, all 100 runs successfully converge to the same maximum energy configuration\. The detailed deflation process for a representative run is shown in Fig\.[S4](https://arxiv.org/html/2608.17135#S4.F4)\(b\), illustrating how ITNT identifies the extremum within TT framework\.

This task involves a 144\-site system spanning a configuration space of2144∼10432^\{144\}\\sim 10^\{43\}\. Performance\-wise, a single self\-multiplication of aχ=409\\chi=409TT takes∼2000\\sim 2000seconds on an Intel Core Ultra 9 275HX CPU\. which has been used for all computational benchmarks\. The subsequent deflation steps take between 3 and 30 seconds as the bond dimension grows from 409 to 1209\. Finally, this extremum\-location ITNT can be generalized to find configurations with energies close to arbitrary target values\. This capability is discussed in detail in the following, where we sample theE=0E=0configurations of the Ising model and find a configuration𝝈∈\{0,1\}n\\bm\{\\sigma\}\\in\\\{0,1\\\}^\{n\}, with the corresponding indexσ\\sigma, that satisfies\|E\[σ\]\|<10−4\|E^\{\[\\sigma\]\}\|<10^\{\-4\}\.

![Refer to caption](https://arxiv.org/html/2608.17135v1/figures/12by12Ising.png)Figure S4:Accuracy and scaling benchmarks\.\(a\) Binary\-discard deflation for locating extrema of a12×1212\\times 12classical Ising model\. The algorithm is performed for 100 runs \(each including 800 extremum\-finding ITNT stepsii\), and the cumulative best value at each step is shown, with the percentages indicating the fraction of runs that attain the corresponding value\. \(b\) A single run of the extremum\-finding ITNT steps\. We fixχ=409\\chi=409and iteratively self\-multiplyEEthat encodes the energy of all the configurations 8 times\. In the resulting TT, we identified800800sub\-leading peaks, and found the extremum\|E\|m​a​x∼33\.35\|E\|\_\{max\}\\sim 33\.35after thei=324i=324thstep\. Here, the configuration returned by theiithdeflation step is denoted as𝝈⁡\(i\)\\bm\{\\sigma\}\(i\)with corresponding indexσ⁡\(i\)\\sigma\(i\), andE\[σ⁡\(i\)\]E^\{\[\\sigma\(i\)\]\}is the corresponding energy of the found configuration\.### IV\.AImportance of self\-multiplication steps

![Refer to caption](https://arxiv.org/html/2608.17135v1/figures/square_12by12_find_max_deflation_without_self-multiply.png)Figure S5:Extremum\-finding ITNT without preprocessing\.Shown is a single instance of the results obtained from extremum\-finding ITNT steps without performing the self\-multiplication preprocessing procedure\. Clearly, in comparison to Fig\.[S4](https://arxiv.org/html/2608.17135#S4.F4)b\., the values of the found configurations are significantly lower, indicating a failure to detect the global extremum\.To resolve tightly packed spectral features near the maximum in continuous data or degenerate configurations, incorporating self\-multiplication steps before the main extremum\-finding ITNT routine is crucial\. In Fig\.[S5](https://arxiv.org/html/2608.17135#S4.F5), we illustrate the configuration found throughout the extremum\-finding ITNT steps\. As shown, without these self\-multiplication preprocessing steps, the algorithm fails to identify the correct extremum due to the model’s rugged energy landscape\.

### IV\.BApplication in energy matching

![Refer to caption](https://arxiv.org/html/2608.17135v1/figures/square_12by12_energy_matching_iter_1e4.png)Figure S6:Energy matching\.Energy matchingE⁡\(𝝈\)=0E\(\{\\bm\{\\sigma\}\}\)=0for the144144\-qubit model is conducted and10410^\{4\}candidate are identified\. The best match,\|E⁡\(𝝈best\)\|∼1\.38×10−4\|E\(\{\\bm\{\\sigma\}\}\_\{\\text\{best\}\}\)\|\\sim 1\.38\\times 10^\{\-4\}, is found after thei=784i=784thstep\. The second\-best match,\|E⁡\(𝝈next\-best\)\|∼1\.81×10−4\|E\(\{\\bm\{\\sigma\}\}\_\{\\text\{next\-best\}\}\)\|\\sim 1\.81\\times 10^\{\-4\}, is found after thei=5813i=5813rddeflation\. During the self\-multiplication, the bond dimension is fixed at818818for1414iterations \(each taking∼104\\sim 10^\{4\}seconds\)\. During the following deflation, the bond dimension is truncated to19001900whenever it reaches20002000, yielding a per\-extremum\-finding ITNT runtime of∼60\\sim 60seconds\.In this section we detail the*energy\-matching*protocol referenced in the main text \(Results\), whose aim is to find configurations𝝈\\bm\{\\sigma\}withE⁡\(𝝈\)≈0E\(\\bm\{\\sigma\}\)\\approx 0for the same12×1212\\times 12random\-field Ising instance\. Following the notation from Methods[V\.6](https://arxiv.org/html/2608.17135#S5.SS6)\), We encode the energy landscape asEσ=E⁡\(𝝈\)E^\{\\sigma\}=E\(\\bm\{\\sigma\}\), and form the transformed landscapeE′​\(𝝈\)=E​\(𝝈\)2−\|E\|max2E^\{\\prime\}\(\\bm\{\\sigma\}\)=E\(\\bm\{\\sigma\}\)^\{2\}\-\\lvert E\\rvert\_\{\\max\}^\{2\}with\|E\|max=max𝝈⁡\|E⁡\(𝝈\)\|\\lvert E\\rvert\_\{\\max\}=\\max\_\{\\bm\{\\sigma\}\}\\lvert E\(\\bm\{\\sigma\}\)\\rvertestimated via the extrema\-locating ITNT\. SinceE′​\(𝝈\)≤0E^\{\\prime\}\(\\bm\{\\sigma\}\)\\leq 0, its extremum \(most negative\) value identifies the configuration withE⁡\(𝝈\)E\(\\bm\{\\sigma\}\)closest to zero, which we then locate using the same extremum\-finding ITNT steps\.

The landscape ofE′​\(𝝈\)E^\{\\prime\}\(\{\\bm\{\\sigma\}\}\)can be highly degenerate and rugged near its extrema\. Thus, to locate the extremum ofE′​\(𝝈\)E^\{\\prime\}\(\{\\bm\{\\sigma\}\}\), the deflation method is again very useful to mitigate the error arising from the bond\-dimension truncation, allowing us to find a string matchingE⁡\(𝝈\)E\(\{\\bm\{\\sigma\}\}\)closer to00\. We find that after performing the iterative self\-multiplication for 14 times \(with bond dimension fixed at a larger valueχ=818\\chi=818to ensure a better accuracy of the ITNT algorithms\), the random binary fixation finds matching accuracy atE⁡\(𝝈first\-try\)∼10−1E\(\{\\bm\{\\sigma\}\}\_\{\\text\{first\-try\}\}\)\\sim 10^\{\-1\}, while the subsequent iterative deflation \(locating10410^\{4\}sub\-leading peaks\) improves the matching accuracy toE⁡\(𝝈best\)∼1\.38×10−4E\(\{\\bm\{\\sigma\}\}\_\{\\text\{best\}\}\)\\sim 1\.38\\times 10^\{\-4\}\. Notably, in Fig\.[S6](https://arxiv.org/html/2608.17135#S4.F6), the best matching configuration separates from the next\-best configuration,E⁡\(𝝈next\-best\)∼1\.81×10−4E\(\{\\bm\{\\sigma\}\}\_\{\\text\{next\-best\}\}\)\\sim 1\.81\\times 10^\{\-4\}, by a Hamming distance of4545\. This rugged, multi\-modal landscape severely challenges both gradient descent and sampling\-based methods such as Parallel Tempering, which often get trapped in local minima or fail to access degenerate solutions, where our TT deflation method inherently bypasses\.