Paper deep dive

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

Xiao Wang, Tomohiro Hashizume, Pia Siegl, Dieter Jaksch

Year: 2026Venue: arXiv preprintArea: cs.LGType: PreprintEmbeddings: 107

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%

Last extracted: 8/19/2026, 4:17:25 AM

Summary

The paper introduces Iterative Tensor Network Transformations (ITNTs), a framework for performing element-wise evaluation of elementary and nonlinear filtering functions on data encoded as Tensor Trains (TTs). This method operates in the compressed domain, enabling efficient computation on exponentially large datasets. Key applications demonstrated include computing reaction rates in 3D reactive flow fields with high fidelity and solving Max-SAT optimization problems on spaces up to 2^70 configurations by finding extrema through iterative self-multiplication and truncation.

Entities (8)

Tomohiro Hashizume · person · 99%Xiao Wang · person · 99%ITNT · algorithm · 98%Tensor Train · datastructure · 97%Max-SAT · problem · 95%Arrhenius reaction rate · function · 92%Newton-Raphson Method · algorithm · 90%Tensor Cross Interpolation · algorithm · 85%

Relation Signals (6)

ITNT → uses → Tensor Train

confidence 98% · ITNTs... for the element-wise evaluation... on data encoded as tensor trains (TTs)

Xiao Wang → authored → ITNT

confidence 95% · Xiao Wang... These authors contributed equally to this work... Here, we introduce iterative tensor network transformations

Tomohiro Hashizume → authored → ITNT

confidence 95% · Tomohiro Hashizume... These authors contributed equally to this work... Here, we introduce iterative tensor network transformations

ITNT → solves → Max-SAT

confidence 95% · finding extrema in complex optimization problems, such as solving Max-SAT instances

ITNT → computes → Arrhenius reaction rate

confidence 92% · evaluating highly nonlinear elementary and filtering functions... enabling high-fidelity reaction rate computation

ITNT → employs → Newton-Raphson Method

confidence 90% · This is achieved by leveraging... the iteratively converging Newton-Raphson method

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

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 to $2^{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.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

104,880 characters extracted from source content.

Expand or collapse full text

Iterative tensor network transformations for element-wise evaluation of elementary and filtering functions Xiao Wang Affiliation: Clarendon Laboratory, University of Oxford, Parks Road, Oxford OX1 3PU, United Kingdom Email: xw970921@gmail.com Thanks: These authors contributed equally to this work. Tomohiro Hashizume 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, Germany Email: tomohiro.hashizume@uni-hamburg.de Thanks: These authors contributed equally to this work. Pia Siegl Affiliation: Institute for Quantum Physics, University of Hamburg, Luruper Chaussee 149, Hamburg 22761, Germany Affiliation: Institute of Software Methods for Product Virtualization, German Aerospace Center (DLR), Nöthnitzer Straße 46b, 01187 Dresden, Germany Dieter Jaksch Affiliation: Clarendon Laboratory, University of Oxford, Parks Road, Oxford OX1 3PU, United Kingdom 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, 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 (I) finding extrema in complex optimization problems, such as solving Max-SAT instances on spaces up to 2702^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 I Introduction 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 the curse of dimensionality, which has been successfully addressed by tensor network algorithms. Among these, the tensor train (T) is a one-dimensional tensor network that represents a data-encoded array of size N as a product of smaller tensors (1; 2; 3; 4; 5). In many cases, the T structure enables an exponential reduction of the computational and memory requirements to ⁡(log⁡N)O( 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, T algorithms have been extended beyond their quantum origins to classical problems, including fluid dynamics (6; 7; 8; 9; 10; 11), plasma physics (12), and machine learning (13). 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; 15; 16). 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 T representation. This requirement, however, runs against a structural limitation of the format. A T 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 T implementation. One route around this limitation is to construct the transformed data by interpolation (17; 18; 19; 20); 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; 22) illustrates, loss of fine structure due to compression or interpolation introduces systematic biases when highly nonlinear observables are evaluated (23; 24). 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), turbulent flow simulations (26; 9; 11), and statistical models encoding combinatorial optimization problems (27). 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 T. ITNTs leverage the iterative application of a small set of efficient T operations: linear operations, (partial) integrations, and element-wise multiplications (Fig. 1 (a)). This enables the computation of functions ranging from elementary functions to highly nonlinear operations, such as filtering, to be applied over N points 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 T form, giving rise to an order-of-magnitude higher fidelity in comparison to the interpolation method (Fig. 1 (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 (c)). Notably, we identify a solution matching those found by state-of-the-art heuristic algorithms (28) from a configuration space of N=270≈1021N=2^70≈ 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. 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, T, directly on T-encoded data Ψ defined on array indices ∈0,1,…,N−1=2n−1 j∈\0,1,…,N-1=2^n-1\. This is achieved by leveraging 2n2^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). The ITNT method is used to compute the highly nonlinear Arrhenius reaction rate (exp(−Ec/T) (-E_c/T)). The plot shows the iso-rate surface at 10% of the maximum reaction rate, computed directly from a compressed temperature field T [K] using a normalized reaction energy of Ec=17865.2E_c=17865.2 [K]. The comparison uses a T bond dimension of χ=200χ=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_opt) and sub-optimal peaks (EsubE_sub), within the problem-encoded energy landscape. Simultaneously, T 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 T format. I Results Let Ψ denote an array of length N, and Ψ[j] ^[j] (j∈0,⋯,N−1j∈\0,·s,N-1\) be the jjth element of this array. Here, square brackets are used for superscripts to distinguish them from exponents. The T representation compresses the element-wise numerical representation of Ψ into a product of n local rank-33 tensors Ψ⁡[q] [q], such that Ψ[j]≈∑Ψ[0]α0,α1[j0]Ψ[1]α1,α2[j1]⋯Ψ[n−1]αn−1,αn[jn−1] ^[j]≈ _ α [0]^[j_0]_ _0, _1 [1]^[j_1]_ _1, _2·s [n-1]^[j_n-1]_ _n-1, _n. In this formulation, jq∈0,1j_q∈\0,1\ corresponds to the qqth bit of the binary representation of index j, denoted by the sequence ∈j0,j1,⋯,jn−1 j∈\j_0,j_1,·s,j_n-1\, and αr _r corresponds to the bond index of the rrth bond. The summation runs over all the bond indices =α0,α1,⋯,αn α=\ _0, _1,·s, _n\, each of which may have a size up to χ, the maximum bond dimension that controls the level of compression (Methods V.1). While a T can exactly represent any array Ψ given an exponentially large maximum bond dimension χ=2n/2χ=2^n/2, in practice, χ is kept at a much smaller constant, e.g., χ=200χ=200 for Fig. 1 (b) and provides an approximation of Ψ[j] ^[j]. This value allows us to balance the minimization of compression loss against the reduced memory requirements ⁡(nχ2)O(nχ^2), and the computational scaling ⁡(poly⁡(n)poly(χ))O(poly(n)poly(χ)). To proceed, we introduce the elementary T operations that constitute the ITNT framework, as illustrated in Fig. 1 (a): scalar multiplication, variational addition, element-wise products, and partial integration. Specifically, scalar multiplication is defined as κΨκ , which is implemented by applying the scalar κ to a single local tensor of the T. Addition of two TTs Ψ1+Ψ2 _1+ _2 is a well-established, stable routine and requires a computational cost scaling as ⁡(nχ3)O(nχ^3) (4). The element-wise (Hadamard) product, Ψ1⊙Ψ2 _1 _2, is realized with a zip-up algorithm using copy tensors (29), resulting in a scaling of ⁡(nχ4)O(nχ^4) (30) or allowing inflation in the intermediate bond dimension, ⁡(nχ3)O(nχ^3) (20). Finally, the partial integration operator ℐ(Ψ,q)=∑\jq∑Ψ[0]α0,α1[j0]Ψ[1]α1,α2[j1]⋯Ψ[n−1]αn−1,αn[jn−1]I( ,q)= _ j \j_q\ _ α [0]^[j_0]_ _0, _1 [1]^[j_1]_ _1, _2·s [n-1]^[j_n-1]_ _n-1, _n constructs a two-element array by summing all elements for which jq=0j_q=0 and jq=1j_q=1, respectively (Methods V.2). 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 T (31), particularly those whose inverses are readily computable from the aforementioned primitives, we employ the Newton-Raphson method (Fig. 1 (a), bottom right). This involves computing the root Ψ′ that satisfies −1(Ψ′)−Ψ=0T^-1( )- =0, where −1(Ψ)T^-1( ) is the inverse transformation satisfying −1((Ψ))=ΨT^-1(T( ))= . Thus, the elementary functions defined by Bourchtein et al. (31) can be fully implemented within the T framework (for details, see Methods V.3). To demonstrate the robustness and accuracy of this approach, in addition to the methane/air jet flame in Fig. 1 (b), we apply this method to solve the Kidder equation (32; 33; 34) in SM I. 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) 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), which implements comparisons using only addition, scalar multiplication, and element-wise multiplication. The T Ψ is compared to a zero field ∅ (∅[j]=0 ^[j]=0 for all j), yielding an element-wise evaluation of the sign function with the outcome SignΨ[j]=Sign⁡Ψ[j]Sign\ \^[j]=Sign\ ^[j]\, where Signx=1Sign\x\=1 (x>0x>0), Signx=0Sign\x\=0 for x=0x=0, and Signx=−1Sign\x\=-1 (x<0x<0). The implementation is efficiently achieved via the iteration R(k+1)=F⁡(R(k))R^(k+1)=F (R^(k) ) with R(0)=ΨR^(0)= where F⁡(R(k))=−12((R(k))⊙3−3R(k)) F(R^(k))=- 12 ((R^(k)) 3-3R^(k) ) (1) and Ψ⊙θ θ denotes the θ element-wise power. Convergence of this iteration to SignΨSign\ \ requires |Ψ[j]|<3| ^[j]|< 3 (∀j∀ j). Therefore, Ψ might first need to be rescaled to ensure that the iterations remain within the radius of convergence. Preconditioning Ψ further may speed up convergence (for details, see Methods subsection V.4). 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 Ψ with Ψ[j]≥0 ^[j]≥ 0, the toal integral obeys ∑jq∈0,1ℐ(Ψ,q)[jq]=∑jΨ[j] _j_q∈\0,1\I( ,q)^[j_q]= _j ^[j]; if this sum is dominated by the maximum element max⁡Ψ[j] \ ^[j]\, then argmaxℐ[jq](Ψ,q)argmax\I^[j_q]( ,q)\ coincides with the qqth component of the global maximum index J=argmaxΨ[j]J=argmax\ ^[j]\. We therefore define the binary fixation (BF) as sequentially fixing the qqth index of a T to mq=argmax⁡ℐ[jq](Ψ,q)m_q=argmax\I^[j_q]( ,q)\, starting from the 00th tensor. BF recovers the full index of the maximum =(m0,m1,⋯,mn−1)= m=(m_0,m_1,·s,m_n-1)= J after n iterations. In practice, a raw field rarely holds such a dominance condition. Therefore, as a preprocessing, we propose to raise the field components to a power 2P2^P via self-multiplications by performing P steps of self-multiplications. This concentrates the field onto its largest elements while preserving the location of the global maximum, and for sufficiently large P, the dominance condition holds for the exact Ψ⊙2P 2^P. However, this self-multiplication causes the bond dimension to grow as χ2Pχ^2^P without bond dimension truncation. To keep the bond dimension of the output T manageable, we use singular value decomposition (SVD) for truncation (37). SVD, by design, minimizes the ℓ2 _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 _2 norm. The ℓ2 _2 error, defined as ε⁡(Ψ,ΨSVD)=‖Ψ−ΨSVD‖2 ( , _SVD)=\| - _SVD\|_2, where ΨSVD _SVD is the Ψ T obtained after SVD truncation, bounds the ℓ∞ _∞ error (ϵ∞=‖Ψ−ΨSVD‖∞ _∞=\| - _SVD\|_∞), 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 _f corresponding to the number of indices that give rise to max⁡Ψ \ \, and BF finds one of these indices as long as the accumulated ℓ2 _2 error per step at the ppth self-multiplication step satisfies ε⁡(Ψ⊙2p,(Ψ⊙2p)SVD)≤(max⁡(Ψ⊙2p)[j]−maxsub⁡(Ψ⊙2p)[j])/2 ( 2^p,( 2^p)_SVD)≤( $ 2^p)^[j]\- _sub\( 2^p)^[j]$/ 2 (2) for the gap between max⁡(Ψ⊙2p)[j] \( 2^p)^[j]\ and maxsub⁡(Ψ⊙2p)[j] _sub\( 2^p)^[j]\ the largest element strictly below the maximum (cf. SM I for details). When the required number of self-multiplications is intractably large, the ℓ2 _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 of J attainable 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Ψ′ defl→ , removes the identified configuration from Ψ′ such that Ψ′[j]=Ψ′[j]−δj,rΨ′[r] [j]= [j]- _j,r [r] (see Methods subsection V.5), where r is 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 Ψ is 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 T-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 T algorithms for extremum-finding and sign-function evaluation can be directly applied to field-processing tasks, as illustrated in Fig. 1 (b). Here, we use them to post-process the temperature field of a precomputed methane/air jet flame (26). This demonstrates how ITNT implements nonlinear functions by computing the Arrhenius reaction rate as an example. Figure 2: Nonlinear post processing with ITNTs. (a) Convergence comparison between ITNT and TCI for the Arrhenius equation (Eq. (3)) computation. While TCI suffers from instability around χ=200χ=200, ITNT demonstrates a drop in error in the same region, consistently converging towards the exact result with increasing χ. (b) Relative ℓ2 _2 error of the 1/T1/T field computed exactly (Ψexact _exact) from the original field and by using the ITNT algorithm (ΨT _T) for χ=50χ=50, 100100, 200200, 300300, and 500500 (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 indicates 1/χ1/χ decay 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 at 9090% of the maximum temperature (0.9max⁡T0.9 \T\). (d) Relative error of the maximum temperature, 1−max⁡T/max⁡Texact1- \T_T\/ \T_exact\, found using ITNT for χ=50χ=50, 100100, 200200, 300300, and 500500 (light to dark). (e) Relative error of the volume of the region where T>0.9max⁡TT>0.9 \T\ for χ=50χ=50, 100100, 200200, 300300, and 500500 (light to dark). Here, the value of 0.9max⁡T0.9 \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), Ar=A0exp(−Ec/T), A_r=A_0 (-E_c/T), (3) where T is the temperature and Ec=17865.2E_c=17865.2 [K] is the activation energy of methane combustion. For simplicity, we set the pre-exponential factor A0=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 T representation with a maximum bond dimension χ. Before encoding, we trim the outer grid points so that each grid dimension is discretized to the nearest power of 22. We then apply the elementary-function ITNT (see Methods V.3) 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 χ, where we obtain the T representation of ArA_r with a final bond dimension of χ. Alternatively, for comparison, a Tensor Cross Interpolation (TCI) approach is employed via the xfac library (18) 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ε=10^-12 or the maximum bond dimension of χ is reached. As shown in Fig. 1 (b), ITNT provides a good approximation to the true reaction rate features already at χ=200χ=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 _2 error ε~(Ψexact,ΨT)=‖Ψexact−ΨT‖2/‖Ψexact‖2 ( _exact, _T)=\| _exact- _T\|_2/\| _exact\|_2 (4) between the exact representation Ψexact _exact and the T representation ΨT _T. As shown in Fig. 2 (a), TCI exhibits a region of instability around bond dimension χ=200χ=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 χ are further supported by the rapid quadratic convergence enabled by the Newton-Raphson method, as shown in Fig. 2 (b) for the first step, which computes the reciprocal field (1/T1/T field). Here, even at χ=200χ=200 and after only 5 iterations, ITNT converges rapidly and approximates the true 1/T1/T field to within a relative error below 2% (left panel, red). Furthermore, we observe a 1/χ1/χ decay in the error with respect to the bond dimension (right panel). In Fig. 2 (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) 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). This feature makes standard sampling or learning algorithms inefficient for detecting the maximum. To find the peak temperature without a priori knowledge of the field’s exact representation, we perform the maximum-finding ITNT, executing five rounds of self-multiplication with a controlled bond dimension χ. The subsequent deflation steps (Methods Sec. V.5) immediately find a candidate value within 1% of the true maximum (Fig. 2 (d)). Then, using the sign-transformation ITNT, we compute the volume of the region exceeding 90% of the found maximum temperature for a given χ. To evaluate accuracy, we plot the relative error of the volume against the true volume in Fig. 2 (e). Even for a bond dimension as small as χ=200χ=200, a relative error of less than 1% is achieved after 10 ITNT rounds. The continued convergence with increasing χ demonstrates that this method can perform nonlinear data processing with high accuracy on complex data structures encoded in the T format (see SM I for further analysis on ITNT transformations acting on continuous data). Application I: Solving NP-hard optimization problems Beyond continuous fields, we study combinatorial optimization problems as a use case for extremum-finding on T-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-k-SAT problem with n variables and M clauses, we want to find a bit string =(σ0,σ1,…,σn−1) σ=( _0, _1,…, _n-1) with a corresponding index σ that maximizes the following k-body Ising-type cost function, E[σ]=M−∑=1M∏m=1klm,E^[σ]=M- _C=1^M _m=1^kl_m^C, (5) where lm∈σj(,m),1−σj(,m)l_m^C∈\ _j_(C,m),1- _j_(C,m)\ is chosen based on whether the corresponding literal is negated. Here, ∈1,…,MC∈\1,…,M\ is the index of the logical clause, and j(,m)∈[0,n−1]j_(C,m)∈[0,n-1] denotes the index of the bit appearing in the mmth term of C. Concretely, lml_m^C is defined as a literal-falsity indicator: we set lm=1−σj(,m)l_m^C=1- _j_(C,m) for an unnegated literal, and set lm=σj(,m)l_m^C= _j_(C,m) for a negated literal. With this convention, ∏m=1klm=1 _m=1^kl_m^C=1 if and only if clause C is violated, meaning that E[σ]E^[σ] counts the number of satisfied clauses. Because E[σ]E^[σ] counts the number of satisfied clauses, it effectively defines the energy landscape of an Ising model with variables σj _j. Consequently, the Max-k-SAT problem translates into finding the global maximum of this landscape. We create a T that encodes the energy landscape of the Max-33-SAT problem, which is known to be an NP-hard problem (41), utilizing the Ising-to-T encoding (Methods subsection V.6). We select an instance with n=70n=70 variables and M=700M=700 clauses from the “Max-SAT 2016 competition” dataset (42), specifically, instance s3v70c700-1. While an array representation of the energy landscape has 2702^70 elements, the mapped T consists of 7070 local tensors and has a bond dimension of χ=701χ=701. 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χ=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 problem s3v70c700-1. Three configurations with an energy of E with E=679E=679 are found after the i=1655i=1655th, the i=6442i=6442nd, and the i=30186i=30186th extremum-finding ITNT steps. (c) Scaling of the minimum bond dimension for successful optimization and solution certification. Here, we generate 100 random n-variable, M=αnM=α n-clause Max-3SAT problems for n=10n=10 to 2020, at clause densities α=8α=8 and α=3α=3. Orange dashed lines show the minimum bond dimension required for the maximum search to return the correct optimum after P=⌈log2⁡(nM)⌉+1P= _2(nM) +1 self-multiplication steps. The blue line shows the minimum bond dimension required for the truncation budget in Method V.7 to be obeyed after P steps. The required bond dimension grows with n, saturating the bond dimension ∼2n/2 2^n/2 (black dotted line). To find the maximum number of satisfied clauses and their truth values, the landscape T is first self-multiplied ten times with the bond dimension fixed at χ=701χ=701. Subsequently, the extremum-locating ITNT (RBF iterations) is performed 10510^5 times. During these RBF iterations, the maximum bond dimension is set to χmax=1000 _ =1000; once the T exceeds this threshold, it triggers a fallback procedure, truncating the bond dimension back to 900900. In Fig. 3 (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 _tr= _χ<i _i^2, where the singular values λi _i are sorted in descending order, to reach 103210^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 full 2702^70). Through these iterative steps, we identify three configurations with the extremal value E=679E=679 (Fig. 3 (b)). These match the best-known solution for this instance and are separated by Hamming distances of 11, 1212, and 1313: 1001001000110000011001111011010101001000000111010010000000110110010101 1001001000110000011001111011010101001010000111010010000000110110010101 1001001110110000010001111011000101011000011001110010000000111100010101 Importantly, ITNT identifies these solutions while probing only a tiny fraction of the full configuration space. Instead of exhaustively evaluating all 270≈10212^70≈ 10^21 elements, the algorithm visits only ⁡(104)O(10^4) peaks in the landscape. The slightly suboptimal extremum (the next-best solution with E=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 T 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 T 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 after P self-multiplication steps (Methods V.7). 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≠ conjecture (41). In Fig. 3 (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 throughout P=⌈log2⁡(nM)⌉+1P= _2(nM) +1 steps for Max-33-SAT problems with M=αnM=α n clauses for different values of the hardness parameter α. As expected from the P≠ 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 the n=70n=70 Max-3-SAT problem could be solved successfully despite the certification threshold being out of reach. These two results presented in Fig. 3 (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 in n. 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. I Discussion 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 (T). 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; 20; 44). 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 T (45). We anticipate that combining these algorithmic refinements with high-precision arithmetic and GPU acceleration (15; 46) 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), genetic haplotype reconstruction (48), and large-scale financial simulations (49). One long-term direction is to investigate the mapping of ITNT algorithms onto quantum hardware, where T-based methods are already showing significant promise (50; 51; 52). Furthermore, the ability to encode symmetries directly within the T formalism (53) 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. IV Acknowledgements 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)”. V Methods In this section, we provide additional details on the foundations and applications of ITNTs discussed in the main text. First, subsection V.1 describes the SVD-based truncation procedure and sweeping protocol used to maintain low-rank tensor representations. Second, subsection V.2 outlines the partial integration formalism for tensor trains. Third, subsection V.3 explains how to compute elementary functions with ITNT using the Newton-Raphson method. Fourth, subsection V.4 describes preconditioning techniques that accelerate the convergence of element-wise sign function evaluation. Fifth, subsection V.5 details the extremum-finding ITNT scheme using self-multiplication, binary fixation, and deflation. Finally, subsection V.6 explains how the energy landscape of a classical Ising model is encoded in a compressed T representation. V.1 The SVD truncation The origin of compression in the T representation lies in the singular value decomposition (SVD), which identifies the optimal low-rank approximation of a tensor by minimizing the ℓ2 _2 norm of the truncation error (2). 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 χ. The careful choice of χ ensures numerical stability and efficient convergence. Furthermore, when performing two-site SVD truncation on a T 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.2 Partial Integration The partial integration ℐ⁡(Ψ,q)I( ,q) gives a two-element array, where ℐ(Ψ,q)[0]I( ,q)^[0] is the sum of all elements of Ψ[j] ^[j] with the qqth bit of index j equal to 00, and ℐ(Ψ,q)[1]I( ,q)^[1] is the corresponding sum for the bit equal to 11. Formally, ℐ⁡(Ψ,q)I( ,q) is written as ℐ(Ψ,q)[jq]=∑\jq∑Ψ[0]α0,α1[j0]⋯Ψ[q]αq,αq+1[jq]⋯Ψ[n−1]αn−1,αn[jn−1] ( ,q)^[j_q]= _ j \j_q\ _ α [0]_ _0, _1^[j_0]·s [q]_ _q, _q+1^[j_q]·s [n-1]_ _n-1, _n^[j_n-1] (6) where \jq j \j_q\ denotes the bit string j without the qqth bit. V.3 ITNT for realizing elementary functions We define the elementary functions as those specified by Bourchtein et al. (31), comprising polynomial, rational, exponential, trigonometric, and irrational functions (i.e., non-integer powers and logarithms). First, given the T f representing the discretized function f⁡(x)f(x) over N=2nN=2^n grid points, we construct the T R(k)R^(k) which represents the kkth integer power of f⁡(x)f(x), (f⁡(x))k (f(x) )^k, using the following iterative element-wise product R(k+1)=f⊙R(k), R^(k+1)=f R^(k), (7) where R(1)=fR^(1)=f, and the operation h=f⊙gh=f g on TTs f and g, denotes an element-wise product, yielding a resulting T h with elements h[j]=f[j]g[j]h^[j]=f^[j]g^[j]. With this choice of transform, we obtain R(k)R^(k), which is equivalent to f(x)⊙kf(x) k. We assume that the bond dimensions of both f and R(k)R^(k) are χ, 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; 30). Here, the maximum bond dimension is truncated to χ. Under this assumption, the element-wise product has a time complexity of ⁡(nχ4)O(nχ^4) and a memory complexity of ⁡(χ3)O(χ^3). To compute the element-wise exponentiation of the field, f⁡(x)→ef⁡(x)f(x)→ e^f(x), we rely on the definition of exponentiation as a limit ef⁡(x)=liml→∞(1+f⁡(x)2l)⊙2l, e^f(x)= _l→∞ (1+ f(x)2^l ) 2^l, (8) where l∈ℤ+l ^+. We can construct this iteration in the T representation by applying the following transformation l times, R(k+1)=F⁡(R(k))=R(k)⊙R(k), R^(k+1)=F (R^(k) )=R^(k) R^(k), (9) starting from R(0)=+f/2lR^(0)=I+f/2^l. Here, the unit-element T I (i.e., representing a vector with all its elements equal to 11) is constructed from the tensor product of rank-1 unit tensors, and thus has a bond dimension of 11. As l increases, the final T approaches the exponentiated T, liml→∞R(l)=exp⁡(f) _l→∞R^(l)= (f ), which represents ef⁡(x)e^f(x). Although the error in Eq. (8) becomes negligible when 2l≫fmax≡maxx⁡f⁡(x)2^l f_max≡ _xf(x), we still need to restrict the number of iteration steps l to 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 T representation. Specifically, we evaluate the sum of eA⁡(xt)e^A(x_t) over all configurations xtx_t, where each configuration represents a distinct imaginary-time-discretized path. Using similar iterative self-multiplication techniques, we can also construct the unitary phase transformation f⁡(x)→eif(x)f(x)→ e^if(x) for complex-valued T 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). 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 T representation (cf. Ref. (54) for a modern implementation). We further expand the ITNT framework by incorporating the Newton-Raphson method. The method obtains the solution R to the equation g⁡(R)=0g(R)=0 iteratively according to R(k+1)=R(k)−g⁡(R(k))g′(R(k)),R^(k+1)=R^(k)- g(R^(k))g (R^(k)), (10) where g′(R)g (R) is the derivative of g with respect to R. The Newton-Raphson method exhibits quadratic convergence as long as the function g⁡(R)g(R) is differentiable at its root. This Newton-Raphson method can be used to perform an elementary nonlinear transformation f⁡(x)→⁡(f⁡(x))f(x) (f(x)) at each point of the initial function f⁡(x)f(x). To achieve this for an arbitrary function R⁡(x)R(x), we define g⁡(R⁡(x))≡−1(R⁡(x))−f⁡(x),g(R(x)) ^-1(R(x))-f(x), (11) where −1(⋅)T^-1(·) denotes the inverse transformation of ⁡(⋅)T(·). Thus, the solution to g⁡(R⁡(x))=0g(R(x))=0 is given by R⁡(x)=⁡(f⁡(x))R(x)=T(f(x)). However, this solution R⁡(x)R(x) can also be determined as the fixed point of the Newton-Raphson iterative routine in Eq. (10). 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)- T^-1(R^(k)(x))-f(x)(T^-1) (R^(k)(x)), (12) where (−1)′(T^-1) denotes the derivative of −1T^-1 with respect to R(k)(x)R^(k)(x). Thus, we find that our targeted transformation is given by ⁡(f⁡(x))=limk→∞R(k)(x)T(f(x))= _k→∞R^(k)(x), provided that the initial function R(0)(x)R^(0)(x) is properly chosen to ensure convergence. The above Newton-Raphson equation (Eq. (12)) can be simplified and then iteratively computed in the T representation from the known f⁡(x)f(x) and a properly chosen R(0)R^(0). We explain this procedure in the following examples, where the target nonlinear transformation ⁡(⋅)T(·) 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 function f⁡(x)f(x) using ITNT. To this end, the equation is reformulated as a root-finding procedure for g⁡(R⁡(x))=(R⁡(x))−1−f⁡(x)=0g(R(x))=(R(x))^-1-f(x)=0, defined in Eq. (11). The resulting Newton-Raphson iteration becomes R(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 T representation, using the element-wise product and the unit-element T I, as follows R(k+1)=R(k)⊙(−f⊙R(k)+2).R^(k+1)=R^(k) (-f R^(k)+2I ). (13) Here, we choose R(0)=f/cR^(0)=f/c as an initial guess, where c is a suppression constant that ensures that the iterations converge. To compute the square root, R⁡(x)=f⁡(x)R(x)= f(x), of a function f⁡(x)f(x) using ITNT, we again follow Eq. (11) and define g⁡(R)=R2−f⁡(x)=0g(R)=R^2-f(x)=0. The corresponding Newton-Raphson iteration becomes R(k+1)(x)=12(R(k)(x)+f⁡(x)R(k)(x))R^(k+1)(x)= 12 (R^(k)(x)+ f(x)R^(k)(x) ), which translates into the T representation as R(k+1)=12(R(k)+f⊙(1/R(k))).R^(k+1)= 12 (R^(k)+f (1/R^(k) ) ). (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)= f(x), can be computed by defining g⁡(R)=eR−f⁡(x)g(R)=e^R-f(x), resulting in the iterative transformation R(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 T representation as follows R(k+1)=R(k)+f⊙exp⁡(−R(k))−.R^(k+1)=R^(k)+f (-R^(k) )-I. (15) We start from R(0)=0R^(0)=0. This ensures that ln⁡(f⁡(x))=R(∞)(x) (f(x))=R^(∞)(x). We then compute the exponentiated T, exp⁡(−R(k)) (-R^(k) ), using the ITNT algorithm defined in Eq. (8). Convergence is guaranteed when f⁡(x)>0f(x)>0 for all x. V.4 Accelerating element-wise sign function evaluation It was shown for the sign transformation (36) that applying the transformation G⁡(R(k))=−1210(1359((R(k))⊙3)−2126R(k)) G(R^(k))=- 12^10 (1359((R^(k)) 3)-2126R^(k) ) (16) before iterating with F (Eq. (1)) further accelerates convergence. In practice, applying the auxiliary transformation G twice serves as an optimal preconditioning step, allowing the T f to converge more rapidly to SignfSign\f\. The numerical performance of this iterative composition is illustrated in Fig. 4 (b). Starting from the identity function f⁡(x)=xf(x)=x, the composition F(3)∘G(2)F^(3) G^(2) rapidly develops a sharp, step-like profile near the origin, effectively approximating the sign function within the compressed T 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 composition F∘GF G is 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.5 Extremum-finding ITNT with a limited number of self-multiplications Figure 4: ITNTs for element-wise nonlinear operations in T representation. (a) Extrema detection with T. i Starting from a T representation of Ψ (black curve), repeated element-wise self-multiplication sharpens the state into a multi-peak profile Φ , where its integral is dominated by extremal configurations (red curve). i 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 between j2=0j_2=0 (orange) and j2=1j_2=1 (cyan) is performed, identifying j2=1j_2=1 as a candidate for the optimal configuration. In the subsequent next panel, j1=0j_1=0 (orange) and j1=1j_1=1 (cyan) are compared and j1=0j_1=0 is selected, while the integrated out region from the fixation of j2=0j_2=0 is indicated in gray. Continuing these steps with a random sequence of local tensors identifies the index I of the maximum (orange, right most panel). i 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 functions F (Eq. (1)) and G (Eq. (16)) that produce an element-wise sign function Signf[j]=Sign⁡f[j]Sign\f\^[j]=Sign\f^[j]\. As an illustrative example, the blue curve shows the output of ITNT for the input field f⁡(x)=xf(x)=x (black dashed line) after applying G twice followed by repeated iterations of F. The red shaded region marks the regime of input values where the composed iteration F⁡(…F(G⁡(G⁡(f))))F(… 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 (a), in this subsection, we describe the ITNT algorithm for extrema location in more detail. Given an n-tensor T f encoding a real-valued field, we amplify large elements by iterative element-wise self-multiplication, R(k+1) R^(k+1) =R(k)⊙R(k) =R^(k) R^(k) (R(0)=f). (R^(0)=f). (17) We then apply normalization after each step, R(k)←R(k)‖R(k)‖.R^(k)← R^(k)\|R^(k)\|. (18) After k iterations, R(k)R^(k) represents f⊙2kf 2^k up to normalization and typically develops a small set of dominant peaks (cf. Fig. 4 (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 index I from a multi-peak T R, we randomly choose a site ordering (p0,…,pn−1)(p_0,…,p_n-1). At a binary fixation step t, with sites p0,…,pt−1\p_0,…,p_t-1\ already fixed, we compare the magnitudes of the two partial integrations for pt=0p_t=0 and pt=1p_t=1, after all other free tensors have been summed out. We then set ptp_t to the branch with the larger magnitude, and continue this process until all sites are fixed, yielding a candidate index I (Fig. 4 (a)i). Deflation Once a candidate index I yielding a maximum element of R is obtained, we remove the corresponding peak from R by subtracting its contribution (Fig. 4 (a)i) R←R−R[I]δI,R← R-R^[I] _I, (19) where δI _I is a χ=1χ=1 T whose local tensors are given by δ[q]0,0[i]=δiq,Iqδ[q]^[i]_0,0= _i_q,I_q. Here δa,b _a,b is a Kronecker delta, defined as δa,b=0 _a,b=0 (a≠ba≠ b) and δa,b=1 _a,b=1 (a=ba=b). Iterating binary fixation → deflation 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,…,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 for an application in energy matching in the Ising model). V.6 Encoding 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 T. We define the classical Ising model on n sites with the Hamiltonian H=∑i=1n∑j=1iJijσiσj,H= _i=1^n _j=1^iJ_ij _i _j, (20) where σi∈0,1 _i∈\0,1\ are the spin variables. Here, we consider Hamiltonians with only two-body interactions; however, the result trivially generalizes to Hamiltonians with arbitrary m-body interactions. For a given configuration =(σ0,σ1,…,σn−1) σ=( _0, _1,…, _n-1), the energy is given by E⁡()=∑i=1n∑j=1iJijσiσj.E( σ)= _i=1^n _j=1^iJ_ij _i _j. (21) Here, JijJ_ij denotes the Ising couplings for i≠ji≠ j, and JiiJ_i represents the on-site terms, since σi2=σi _i^2= _i. To find the configuration that maximizes E⁡()E( σ) by performing the extremum-finding ITNT, we encode the energy landscape into the T E by setting the values of E at index σ to the energy E⁡()E( σ), i.e., E[σ]=E⁡()E^[σ]=E( σ), where σ is the index obtained by interpreting the configuration as a bit string. The T representation of the energy landscape E⁡()E( σ) is sparse. The bond dimension χ of E is bounded by the number of non-zero coupling terms in JijJ_ij, because we can decompose E as E=∑i=1n∑j=1iE⁡[i,j]E= _i=1^n _j=1^iE[i,j], where E⁡[i,j]E[i,j] is a T encoding the coupling energy corresponding to the term JijJ_ij across all configurations. Crucially, E⁡[i,j]E[i,j] can be recast into a χ=1χ=1 T as follows, for i≠ji≠ j, E[i,j][k]=Jij∑ E[i,j]^[k]=J_ij _ α 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]⋯ A[0]_ _0, _1^[k_0]·s A[i-1]_ _i-1, _i^[k_i-1]B[i]_ _i, _i+1^[k_i]A[i+1]_ _i+1, _i+2^[k_i+1]·s A[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], A[j-1]_ _j-1, _j^[k_j-1]B[j]_ _j, _j+1^[k_j]A[j+1]_ _j+1, _j+2^[k_j+1]·s A[n-1]_ _n-1, _n^[k_n-1], (22) where A⁡[q]A[q] is the unit tensor with all elements equal to 11, and B⁡[q]B[q] is a χ=1χ=1 local tensor which takes the form B[q]0,0[kq]=kqB[q]^[k_q]_0,0=k_q, for a given site q. To encode the energy landscape of any m-body interaction term, E⁡[i1,i2,i3,⋯,im]E[i_1,i_2,i_3,·s,i_m], one creates the unit T I and then replaces the local tensors corresponding to the sites that appear in the interaction term with the B tensor. According to Ref. (2), the sum of two TTs yields another T with a summed bond dimension. Thus, summing over E⁡[i,j]E[i,j] results in the T E having a bond dimension equivalent to the number of terms in E⁡()E( σ), including the constant offset. Because of this sparsity, we perform ITNTs directly on E to 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 T representation, which acts simultaneously across all spin configurations. This approach differs from tensor network simulations of quantum annealing (55), where the T encodes the quantum state of the annealer rather than the energy landscape considered here. V.7 Truncation error budget for certified optimality of Max-3-SAT solutions Finally, we derive a truncation-error bound extending Eq. (2) that certifies the maximum found by the maximum-finding ITNT. First we derive the sufficient number P of self-multiplication steps for the field E that encodes the energy landscape of the M-clause Max-3-SAT problem as defined in Eq. (5). We then define the normalized field Ψ=E/A =E/A, where A=max⁡EA= \E\ is the maximum number of satisfied clauses. A value of P sufficient for BF to retain an optimal branch is obtained by considering the worst case in which every nonoptimal element of Ψ takes its largest possible value, 1−1/A1-1/A. Let J be the index with Ψ[J]=1 ^[J]=1, assumed unique for simplicity; degeneracy of the optimum only increases with the margin below. The sum of all elements of Ψ⊙2P 2^P other than J then obeys ∑j≠J(Ψ⊙2P)[j]≤2n(1−1/A)2P≤2nexp(−2P/M), _j≠ J( 2^P)^[j]≤ 2^n(1-1/A)^2^P≤ 2^n (-2^P/M ), (23) where the right-hand inequality uses A≤MA≤ M, valid since at most M clauses can be satisfied. It is therefore sufficient to choose P such that 2nexp(−2P/M)<12^n (-2^P/M )<1, and hence we obtain P>log2⁡(nMln⁡2)P> _2(nM 2). For the result in Fig. 3(c) we set P=⌈log2⁡(nM)⌉+1P= _2(nM) +1, which obeys the bound. Given P satisfies the bound, we now derive the truncation budget for this dominance condition to be met at the end of P self-multiplication steps. Let R(p)R^(p) denote the computed field after p≤Pp≤ P steps, with R(0)=ΨR^(0)= and R(p)=(R(p−1))⊙2+τp R^(p)=(R^(p-1)) 2+ _p (24) where τp _p is the total error injected by the truncation at step p. Let tp=ε⁡((R(p−1))⊙2,R(p))=‖τp‖2t_p= ((R^(p-1)) 2,R^(p))= \| _p \|_2 be the ℓ2 _2 norm of the error as defined in Eq. (2) and let Δp=R(p)−Ψ⊙2p _p=R^(p)- 2^p be the error tensor, the accumulated deviation from the true representation of the self-multiplied field. After P iterations, 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 of J correctly. The difference between the partial integral at the qqth site with the correct bit (jq=Jqj_q=J_q) and with the incorrect bit (jq≠Jqj_q≠ J_q) is then bounded from below by Λ(Ψ⊙2P,q)≥1−∑i≠J(Ψ⊙2P)[i]≥1−2nexp(−2P/M), ( 2^P,q)≥ 1- _i≠ J( 2^P)^[i]≥ 1-2^n (-2^P/M), (25) where Λ⁡(X,q)=ℐ(X,q)jq=Jq−ℐ(X,q)jq≠Jq (X,q)=I(X,q)^j_q=J_q-I(X,q)^j_q≠ J_q for a given field X is the difference between the branches, one with the solution and the other without, and the second inequality holds since J contributes 11 to the first integral and every remaining configuration contributes non-negatively to the second. Invoking the triangular inequality, for R(p)R^(p) to identify J correctly it therefore suffices that Λ(R(P),q)≥1−2nexp(−2P/M)−|Λ(ΔP,q)|>0. (R^(P),q)≥ 1-2^n (-2^P/M)- | ( _P,q) |>0. (26) From the definition, Δp _p obeys the recursion Δp=2Ψ⊙2p−1⊙Δp−1+Δp−1⊙2+τp. _p=2 2^p-1 _p-1+ 2_p-1+ _p. (27) From this, provided that the initial truncation is 00, we derive ΔP _P =2P−pΨ⊙(2P−2p)⊙Δp+∑m=p+1P2P−mΨ⊙(2P−2m)⊙Δm−1⊙2+∑m=p+1P−12P−mΨ⊙(2P−2m)⊙τm =2^P-p (2^P-2^p) _p+ _m=p+1^P2^P-m (2^P-2^m) _m-1 2+ _m=p+1^P-12^P-m (2^P-2^m) _m =∑m=1P2P−mΨ⊙(2P−2m)⊙Δm−1⊙2+∑m=1P−12P−mΨ⊙(2P−2m)⊙τm, = _m=1^P2^P-m (2^P-2^m) _m-1 2+ _m=1^P-12^P-m (2^P-2^m) _m, (28) which is merely a telescoping error propagation. Finally, let up=2up−1+up−12+tpu_p=2u_p-1+u^2_p-1+t_p with u0=0u_0=0, combining the linearity of the partial integration with inequality Eq. (28) and applying the bounds |Λ⁡(Ψ⊙(2P−2m)⊙τm,q)|≤‖Ψ⊙(2P−2m)‖2‖τm‖2≤2tm | ( (2^P-2^m) _m,q)|≤\| (2^P-2^m)\|_2\| _m\|_2≤ 2t_m (29) and |Λ⁡(Ψ⊙(2P−2m)⊙Δm−12,q)|≤‖Ψ⊙(2P−2m)‖∞‖Δm−1‖2≤um−12, | ( (2^P-2^m) _m-1^2,q)|≤\| (2^P-2^m)\|_∞\| _m-1\|_2≤ u_m-1^2, (30) we obtain |Λ⁡(ΔP,q)|≤∑m=1P−12P−m2tm+∑m=1P2P−mum−12. | \! ( _P,q ) |≤ _m=1^P-12^P-m 2\,t_m+ _m=1^P2^P-mu_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), the optimality of the solution is certified provided ∑m=1P−12P−m2tm+∑m=1P2P−mum−12<1−2nexp⁡(−2PM). _m=1^P-12^P-m 2\,t_m+ _m=1^P2^P-mu_m-1^2<1-2^n (- 2^PM ). (32) This is exactly the bound that the truncation errors tmt_m must obey over the full P steps 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 Experiment 2007, P10014 (2007). Schollwöck (2011) U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011). Orús (2014) R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals of Physics 349, 117 (2014), arXiv: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 Physics 411, 167998 (2019), arXiv: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 Physics 5, 273 (2023). Lubasch et al. (2018) M. Lubasch, P. Moinier, and D. Jaksch, Multigrid renormalization, Journal of Computational Physics 372, 587 (2018). 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 Science 2, 30 (2022), arXiv:2106.05782 [physics, physics:quant-ph] . Kiffner and Jaksch (2023) M. Kiffner and D. Jaksch, Tensor network reduced order models for wall-bounded flows, Physical Review Fluids 8, 124101 (2023). Gourianov et al. (2025) N. Gourianov, P. Givi, D. Jaksch, and S. B. Pope, Tensor networks enable the calculation of turbulence probability distributions, Science Advances 11, eads5990 (2025). 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 (2025), arXiv:2507.05222 [physics.flu-dyn] . 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 (2026), arXiv:2604.16179 [physics] . Huynh et al. (2023) L. Huynh, J. Hong, A. Mian, H. Suzuki, Y. Wu, and S. Camtepe, Quantum-Inspired Machine Learning: A Survey (2023), arXiv:2308.11269 [cs] . 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 X 8, 031012 (2018). 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 Research 7, 013112 (2025). 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 Research 8, 023006 (2026). Michailidis et al. (2024) A. A. Michailidis, C. Fenton, and M. Kiffner, Tensor Train Multiplication (2024), arXiv:2410.19747 [physics] . 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 Physics 18, 104 (2025), arXiv:2407.02454 [physics] . 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 Physics 529, 113891 (2025). 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 Monthly 100, 768 (1993), 2324783 . Bullen (2003) P. S. Bullen, Handbook of Means and Their Inequalities (Springer Netherlands, Dordrecht, 2003). Jarzynski (1997) C. Jarzynski, Nonequilibrium equality for free energy differences, Physical Review Letters 78, 2690 (1997). Ruel and Ayres (1999) J. J. Ruel and M. P. Ayres, Jensen’s inequality predicts effects of environmental variation, Trends in Ecology & Evolution 14, 361 (1999). 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, Nature 596, 583 (2021). 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 E 110, 045308 (2024). 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 Advances 1, 042172 (2011). Stoudenmire and White (2010) E. M. Stoudenmire and S. R. White, Minimally Entangled Typical Thermal State Algorithms, New Journal of Physics 12, 055026 (2010), arXiv:1002.1305 [cond-mat] . 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 Mathematics 135, 63 (2015). 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 Journal 74, 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 Plus 133, 28 (2018). 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, edited by S. Moriai and H. Wang (Springer International Publishing, Cham, 2020) p. 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) p. 221–256. Oseledets (2011) I. V. Oseledets, Tensor-Train decomposition, SIAM Journal on Scientific Computing 33, 2295 (2011). 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 Access 8, 13937 (2020). 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 (2016). Michailidis et al. (2025) A. A. Michailidis, C. Fenton, and M. Kiffner, Element-wise multiplication of tensor trains, SIAM Journal on Scientific Computing 47, B1158 (2025). 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 Quantum 4, 040318 (2023). 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 . Jiang et al. (2018) S. Jiang, K. A. Britt, A. J. McCaskey, T. S. Humble, and S. Kais, Quantum Annealing for Prime Factorization, Scientific Reports 8, 17667 (2018). 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 (2026), arXiv: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 Research 114, 219 (1999). 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 Letters 132, 040404 (2024), arXiv:2307.01696 [quant-ph] . 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 (2025), arXiv:2502.04425 [quant-ph] . McCulloch and Gulácsi (2002) I. P. McCulloch and M. Gulácsi, The non-Abelian density matrix renormalization group algorithm, Europhysics Letters 57, 852 (2002). 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 (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 1 Clarendon 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 I Application in nonlinear PDE In this section, we demonstrate the application of the Newton-Raphson-based ITNT method in subsection V.3 to solve nonlinear PDEs in the T representation, with particular focus on the Kidder equation, which models flow pressure in porous tubes. I.A Using the Newton-Raphson Method to realize the nonlinearity in Kidder equation To solve the Kidder equation, we need to realize the element-wise transformation w⁡(x)→(1−κw(x))−12w(x)→ (1-κ w(x) )^- 12 in T representation. Following subsection V.3, we obtain the solution of g⁡(R⁡(x))=R(x)−2−1−κ−w⁡(x)=0g(R(x))= R(x)^-2-1-κ-w(x)=0, via the iteration: R(k+1)(x)=R(k)(x)−R(k)(x)−2−(1−κw(x))−2R(k)(x)−3,R^(k+1)(x)=R^(k)(x)- R^(k)(x)^-2- (1-κ w(x) )-2R^(k)(x)^-3, Simplifying this expression, we obtain the following iteration in T representation: R(k+1)=12R(k)⊙(3−(1−κw)⊙R(k)⊙R(k)).R^(k+1)= 12R^(k) (3I- (1-κ w ) R^(k) R^(k) ). (S1) Here, R(k)R^(k) and w are T encoding the function R(k)R^(k) and w⁡(x)w(x) on the equally discretized x-grid. The term I is a χ=1χ=1 unit T where all the element of its local tensors are 11 as explained in the main text. For |κ|<1|κ|<1 and |w⁡(x)|<1|w(x)|<1, the iteration converges when starting from a unit-element T R(0)=R^(0)=I which represents R(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 T, limk→∞|R(k)⟩ _k→∞ R^(k), accurately represents (1−κw(x))−12(1-κ w(x))^- 12. I.B Conducting the self-consistent Picard iteration The Kidder equation is given by: w′(x)+ 2xw′(x)(1−κw(x))−12= 0,w (x)\;+\;2\,x\,w (x)\, (1-κ\,w(x) )^- 12\;=\;0, subject to the boundary conditions w⁡(0)= 1andlimx→∞w⁡(x)= 0.w(0)\;=\;1 _x→∞w(x)\;=\;0. To solve it numerically, we discretize the interval [0,xmax][0,x_ ] into 2n2^n grid points xix_i, with inter-grid point spacing Δx=xmax/2n x=x_ /2^n. Here, n represents the number of tensors required to store the discretized field in the T format. For the Kidder equation with κ=0.5κ=0.5, it is sufficient to set xmax=16x_ =16 to approximate xmax≈∞x_ ≈∞ (33), leading to the grid resolution of Δx=2−12 x=2^-12. We label the approximate solution by wi≈w⁡(xi)w_i≈ w(x_i) and use finite differences to approximate w′(x)w (x) and w′(x)w (x). For the interior points i=1,…,2ni=1,…,2^n, the discrete system typically reads as follows: Fi()=wi+1−2wi+wi−1(Δx)2+ 2xiwi+1−wi−12Δx(1−κwi)−12= 0. splitF_i(w)&\;=\; w_i+1-2\,w_i+w_i-1( x)^2\\ &~~\;+\;2\,x_i\, w_i+1-w_i-12\, x\, (1-κ\,w_i )^- 12\;=\;0. split (S2) We keep w0=1w_0=1 and w2n=0w_2^n=0 fixed to enforce the boundary conditions. A Picard fixed-point iteration solves this nonlinear system by rewriting each equation to isolate wiw_i on one side. Concretely, one arranges: wi=Gi(wi−1,wi,wi+1),w_i\;=\;G_i (w_i-1,\,w_i,\,w_i+1 ), so that each interior variable wiw_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, splitG_i (w^(k) )&= w_i+1^(k)+w_i-1^(k)2\\ &+ x_i\, x2\, [w_i+1^(k)-w_i-1^(k) ] (1-κ\,w_i^(k) )^- 12, split (S3) where wi−1(k),wi+1(k)w_i-1^(k),\,w_i+1^(k), and wi(k)w_i^(k) come from the previous iteration. In T representation, this Picard iteration becomes G⁡(w(k))=12(L^+R^)⋅w(k)+Δx2X⊙(1−κw(k))−12⊙(L^−R^)⋅w(k).. split&G (w^(k) )= 12 ( L+ R )· w^(k)\\ &+ x2X (1-κ\,w^(k) )^- 12 ( L- R )· w^(k). split. (S4) Here, X is a T that encodes the coordinates of the discretized grid (X[i]=iΔxX^[i]=i x)k, L L and R R are linear operators that shift the index by one, i.e. Ψ′=L^Ψ = L s.t. Ψ′[j]=Ψ[j+1] [j]= ^[j+1] and Ψ′=R^Ψ = R s.t. Ψ′[j]=Ψ[j−1] [j]= ^[j-1] (6). The T representation of the term (1−κwi(k))−12(1-κ w_i^(k))^- 12, denoted by (1−κwOPEN(k))−12CLOSE (1-κ\,w^(k) )^- 12, is given by ITNTs in Eq. (S1). We start from an initial guess w(0)w^(0) and perform the iteration w(k+1)=G(w(k)),k=0,1,2,…w^(k+1)=G (w^(k) ), k=0,1,2,… until successive solutions are sufficiently close. This fixed-point iteration converges to the unique T w satisfying the discretized Kidder equation and all of its boundary conditions. I.C Anderson 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 iteration n: w(k+1)=G⁡(w(k)).w^(k+1)=G (w^(k) ). Define the residual as f(k)=G⁡(w(k))−w(k).f^(k)=G (w^(k) )-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 _i that minimize the new residual. A commonly used variant enforces ∑iαi=1 _i _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)\;=\; _i _i\, (w^(k-i)+f^(k-i) ),~with~ _i _i=1. The small system to find a set of scalars αi\ _i\ by computing dot products among the stored residuals. Once those αi _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 produce w(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., storing m previous iterates) to keep the least-squares problem small. In a T-based Kidder solver, the main computational effort lies in evaluating inner products of residuals (k)f^(k) stored in tensor form and carefully truncating the T 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.D Results We set the SVD tolerance to 10−1210^-12 with χ=400χ=400, and conduct the Picard iteration for 10510^5 times with an Anderson acceleration with a memory of m=10m=10. After this, we prolong the T to a finer grid (see Ref. (6)) on which we re-run the above iteration. We start from a grid finesse of Nres=4N_res=4 and end at Nres=17N_res=17, the final result provides an initial slope of -1.1917912520439131, which shows an error of ∼6×10−7 6× 10^-7 compared to the standard result given in Ref. (33). 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 10^-14. In this case, each iteration step becomes slower as the bond dimension increases. Besides, as the grid finesse NresN_res reaches ∼18 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-κ w(x))^-0.5 via the Newton-Raphson method in the T representation, and this approach is highly effective. I Truncation error bounding For a positive field Ψ , with (⋯)SVD(·s)_SVD denoting the field after SVD truncation, the maximum and the sub-maximum of Ψ⊙2p 2^p retain their order provided ε⁡(Ψ⊙2p,(Ψ⊙2p)SVD)≤(max⁡(Ψ⊙2p)[j]−maxsub⁡(Ψ⊙2p)[j])/2 ( 2^p,( 2^p)_SVD)≤( $ 2^p)^[j]\- _sub\( 2^p)^[j]$/ 2 (S5) as presented in Eq. (2). In this appendix, we derive the right-hand side of this inequality by showing that 1/21/ 2 is the largest factor for which order preservation is guaranteed. Let Φ=Ψ⊙2p = 2^p and Φ′=(Ψ⊙2p)SVD =( 2^p)_SVD for the exact and truncated arrays. For simplicity we write ε=‖Φ−Φ′‖2 = \| - \|_2 for the truncation error and Δ=max⁡Φ−maxsub⁡Φ>0 = \ \- _sub\ \>0 for 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⁡Φ _1= \ \- \ \, _2= _sub\ \- _sub\ \ (S6) which measure the depression of the maximum (δ1 _1) and the elevation of the second maximum (δ2 _2). In order for the order of the maximum and second maximum to be conserved, δ1 _1 and δ2 _2 must meet δ1+δ2≤Δ _1+ _2≤ , while the truncation error irself constraints them by δ12+δ22≤ε2 _1^2+ _2^2≤ ^2. Invoking Cauchy-Schwarz inequality, we obtain δ1+δ2≤2δ12+δ22≤2ε≤Δ _1+ _2≤ 2 _1^2+ _2^2≤ 2 ≤ (S7) where the first inequality is saturated when δ1=δ2 _1= _2, and the second when δ1=δ2=ε/2 _1= _2= / 2. With this, Eq. (2) is shown. I Application 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 f1(x)f_1(x) ITNT(f1(x))ITNT(f_1(x)) Figure S1: Result of the Sign-function transformation via ITNT of the function defined in Eq. S8 for different number of roots NrootN_root and different qubit numbers n. For all examples we consider Niter=20N_iter=20 and an SVD tolerance of 10−810^-8. Sign transformation of the function with the following properties: (a) n=14n=14, Nroot=7N_root=7, χmax=20 _max=20. (b) n=14n=14, Nroot=7N_root=7, χmax=25 _max=25. (c) n=18n=18, Nroot=7N_root=7, χmax=25 _max=25. (d) n=18n=18, Nroot=23N_root=23, χmax=75 _max=75. For Nroot=23N_root=23, a noiseless sign transformation requires χmax=75 _max=75. Figure S2: Scaling of the element-wise evaluation of sign function via ITNT (a) Element-wise evaluation of sign function y=Signfω(x)y=Sign\f_ω(x)\ (red) using ITNT of a function y=fω(x)y=f_ω(x) (black) ω=1.0ω=1.0. (b) Required bond dimension in dependence on the number of tensors n for ω=0.5ω=0.5. Considered are the maximal intermediate bond-dimension χpeak _peak (crosses) used by the algorithm if the bond dimension is not restricted and the SVD tolerance is fixed at 10−810^-8. Furthermore, we consider the final bond dimension χfinal _final (circles) required to represent y accurate up to the order of the SVD tolerance of 10−810^-8 with respect to the L2L^2-error. Finally, there is the cutoff bond dimension χmin _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 _final within 200200 iteration steps. (c) Required bond dimensions for different numbers of roots NrootsN_roots by using ω=[0.5,1.0,1.5,2.0]ω=[0.5,1.0,1.5,2.0] corresponding to Nroots=[7,13,19,25]N_roots=[7,13,19,25] and n=18n=18 to converge within a maximum of 200200 iteration steps. The bond dimensions are defined as for (b). (d) The number of iteration steps NiterN_iter needed to reach the final compressed result with respect to the chosen maximal cutoff bond dimension χmax _ max for ω=0.5ω=0.5 for n=14,18,22n=14,18,22 (light to dark). All simulations were implemented using the zip-up algorithm for point-wise multiplication. Figure S3: The bond dimension of the sign-transformed T after NiterN_iter steps of iteration. The qubit number is n=18n=18 with root number Nroot=7N_root=7, and χmax _max denotes the maximum bond dimension allowed during T operations (e.g., element-wise products), with the SVD tolerance fixed at 10−810^-8. For χmax=28 _max=28, the convergence is reached at Niter≈64N_iter≈ 64. For χmax≳35 _max 35, the required iteration count reduces to Niter≈34N_iter≈ 34. fω(x)=cos⁡(ωx)−0.000035(x−2)3+0.001x2−0.02f_ω(x)= (ω x)-0.000035(x-2)^3+0.001x^2-0.02 (S8) encoded across n tensors, where x∈[0,64)x∈[0,64) and the frequency ω determines the number of roots on this interval. Varying the number of roots, we study the impact of the chosen cutoff bond dimension χmax _ max on the convergence of the algorithm. The result of the ITNT-based sign-transformation of fω(x)f_ω(x) is depicted in Fig. S1 for different values of ω, and hence different number of roots NrootN_ root, and with a changing bond-dimension. For the simulation shown in Fig. S1 (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 (b)-(d), the chosen χmax _ max is sufficient for convergence. Hereby, the necessary χmax _ 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 (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 _ final up to an accuracy of order 10−810^-8, saturates quickly with the number of tensors while showing a linear dependence on NrootN_ root. Furthermore, we consider the minimum cutoff bond dimension χmin _ min required during ITNT to allow convergence within 200200 steps. This χmin _ min is significantly lower than the bond dimension χpeak _ peak, the bond dimension used by the algorithm when no cutoff bond dimension is set, and only the SVD tolerance of 10−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 _ min to allow for a quicker convergence (c.f. S2 (c)). In our simulations, we observe that selecting a sufficiently large bond dimension χmax≳5Nroot _max 5N_root is crucial to minimize the convergent iteration step, NiterN_iter. Furthermore, we find that when the maximum bond dimension is large enough (i.e., χmax≳5Nroot _max 5N_root), the sign transformation ITNT fully converges after Niter≈2N_iter≈ 2N steps of iteration. For such a large enough χmax _max, after the ITNT reaches convergence, the bond dimension of the T shrinks back roughly to the root number NrootN_root. We can truncate the bond dimension to χmax _max up to χmax∼4Nroot _max 4N_root (c.f. Fig. S3). In this case, the convergence can still be reached, although it is reached at a larger iteration step NiterN_iter. As we truncate to an even smaller χmax _max, the iteration converges very slowly, and the whole iteration may eventually diverge due to numerical errors. The runtime scalability with qubit number N is evident: as we increase N, the bond dimension and runtime remain manageable. The runtimes required to perform the sign transformation with ITNT are summarized in Tab. 1 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. n NrootN_root χmax _max time (sec) of ITNT-based sign tranforms time (sec) of standard sign transform. 22 7 45 3 1.6 22 19 100 7 1.6 22 51 250 65 1.6 24 7 45 3.5 6.5 24 19 100 8.5 6.5 24 51 250 84 6.5 26 7 45 5 28 26 19 100 11 28 26 51 250 110 28 Table 1: Comparison of runtime performance across different qubit numbers n, root numbers NrootN_ root, and chosen cutoff bond dimensions χmax _ 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 from 10−810^-8 to 10−610^-6 will only slightly increase the error of the ITNT result (i.e., the L2L^2-norm of the difference compared with the element-wise transformed array), from ∼10−8 10^-8 to ∼10−6 10^-6. At the core of the sign-transformation lays the the homomorphic comparison iteration, F⁡(f)=0.5f(3−f2)F(f)=0.5f(3-f^2). In the following we will consider its iteration speed for different distances to the root. Near f=±1f=± 1: Double-exponential Convergence The stable fixed points of F are at f∗=±1f^*=± 1. Near f=1f=1, let ϵ(k)=1−f(k)ε^(k)=1-f^(k) represent the small deviation from the fixed point. The recurrence is approximately ϵ(k+1)≈1.5(ϵ(k))2ε^(k+1)≈ 1.5\,(ε^(k))^2. This leads to double-exponential convergence ϵ(k)=(1.5ϵ(0))2k/1.5ε^(k)=(1.5\,ε^(0))^2^k/1.5. If 1.5ϵ(0)<11.5ε^(0)<1, ϵ(k)ε^(k) decays extremely rapidly, ensuring that the ITNT quickly approaches f∗=1f^*=1. A similar argument applies near f=−1f=-1. Near f=0f=0: Exponential Departure and Resolution Threshold In contrast to the stable fixed points at f=±1f=± 1, the point f=0f=0 is an unstable fixed point. To see this, we linearize F around f=0f=0: F⁡(f)=0.5f(3−f2)≈1.5ffor small f.F(f)=0.5f(3-f^2)≈ 1.5f small f. Thus, any small nonzero value of f near zero grows by approximately 1.5 per iteration, f(k)≈f(0)×(1.5)kf^(k)≈ f^(0)×(1.5)^k. This results in exponential growth away from zero, repelling from f=0f=0. Suppose that after S steps of iteration, the iterative process realizes a sign-like function, i.e., f(S)≈1f^(S)≈ 1, we can estimate the “resolution threshold” for the initial value, f(0)≈(1.5)−Sf^(0)≈(1.5)^-S. This threshold, (1.5)−S(1.5)^-S, sets up a resolution limit in approximating the sign function. After S steps, any initial values f(0)f^(0) (and hence portions of the function f⁡(x)f(x)) larger than (1.5)−S(1.5)^-S are effectively sign-transformed and pushed toward 11, while those smaller than (1.5)−S(1.5)^-S remain evidently far from 1. In other words, after sufficiently S steps 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 function f⁡(x)f(x), the intermediate value theorem implies that whenever there is a change in the sign of f⁡(x)f(x), there exists at least one root of a function. When a curve is represented as an T, roots can be found by first applying the sign-transformation ITNT, then, we apply the differential differential linear operator R^−L R- L for the linear operators R R and L L introduced in Eq. (S4) to obtain the T that marks the root positions as delta peaks. From the resulting multi-delta-peak T, we use the deflation method in Sec. V.5 to iteratively locate each peak position in the T representation. This procedure outputs the set of root positions on the discretized grid. IV Benchmark: performance of extremum finding on discrete problem Here, we focus on a classical Ising model on a square lattice of size 12×1212× 12 with 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 for details of the model). The energy landscape of the classical Ising model can be mapped onto energy-landscape T E with 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χ=409 achieves the exact encoding of the energy landscape. To identify the maximum energy configuration, we fix the Hamiltonian and perform self-multiplication p=8p=8 times, then perform 800 independent deflation runs using 100 different stochastic seeds. In Fig. S4 (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 (b), illustrating how ITNT identifies the extremum within T framework. This task involves a 144-site system spanning a configuration space of 2144∼10432^144 10^43. Performance-wise, a single self-multiplication of a χ=409χ=409 T takes ∼2000 2000 seconds 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 the E=0E=0 configurations of the Ising model and find a configuration ∈0,1n σ∈\0,1\^n, with the corresponding index σ, that satisfies |E[σ]|<10−4|E^[σ]|<10^-4. Figure S4: Accuracy and scaling benchmarks. (a) Binary-discard deflation for locating extrema of a 12×1212× 12 classical Ising model. The algorithm is performed for 100 runs (each including 800 extremum-finding ITNT steps i), 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χ=409 and iteratively self-multiply E that encodes the energy of all the configurations 8 times. In the resulting T, we identified 800800 sub-leading peaks, and found the extremum |E|max∼33.35|E|_max 33.35 after the i=324i=324th step. Here, the configuration returned by the iith deflation step is denoted as ⁡(i) σ(i) with corresponding index σ⁡(i)σ(i), and E[σ⁡(i)]E^[σ(i)] is the corresponding energy of the found configuration. IV.A Importance of self-multiplication steps 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. S4b., 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, 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.B Application in energy matching Figure S6: Energy matching. Energy matching E⁡()=0E( σ)=0 for the 144144-qubit model is conducted and 10410^4 candidate are identified. The best match, |E⁡(best)|∼1.38×10−4|E( σ_best)| 1.38× 10^-4, is found after the i=784i=784th step. The second-best match, |E⁡(next-best)|∼1.81×10−4|E( σ_next-best)| 1.81× 10^-4, is found after the i=5813i=5813rd deflation. During the self-multiplication, the bond dimension is fixed at 818818 for 1414 iterations (each taking ∼104 10^4 seconds). During the following deflation, the bond dimension is truncated to 19001900 whenever it reaches 20002000, yielding a per-extremum-finding ITNT runtime of ∼60 60 seconds. In this section we detail the energy-matching protocol referenced in the main text (Results), whose aim is to find configurations σ with E⁡()≈0E( σ)≈ 0 for the same 12×1212× 12 random-field Ising instance. Following the notation from Methods V.6), We encode the energy landscape as Eσ=E⁡()E^σ=E( σ), and form the transformed landscape E′()=E()2−|E|max2E ( σ)=E( σ)^2- E _ ^2 with |E|max=max⁡|E⁡()| E _ = _ σ E( σ) estimated via the extrema-locating ITNT. Since E′()≤0E ( σ)≤ 0, its extremum (most negative) value identifies the configuration with E⁡()E( σ) closest to zero, which we then locate using the same extremum-finding ITNT steps. The landscape of E′()E ( σ) can be highly degenerate and rugged near its extrema. Thus, to locate the extremum of E′()E ( σ), the deflation method is again very useful to mitigate the error arising from the bond-dimension truncation, allowing us to find a string matching E⁡()E( σ) closer to 00. We find that after performing the iterative self-multiplication for 14 times (with bond dimension fixed at a larger value χ=818χ=818 to ensure a better accuracy of the ITNT algorithms), the random binary fixation finds matching accuracy at E⁡(first-try)∼10−1E( σ_first-try) 10^-1, while the subsequent iterative deflation (locating 10410^4 sub-leading peaks) improves the matching accuracy to E⁡(best)∼1.38×10−4E( σ_best) 1.38× 10^-4. Notably, in Fig. S6, the best matching configuration separates from the next-best configuration, E⁡(next-best)∼1.81×10−4E( σ_next-best) 1.81× 10^-4, by a Hamming distance of 4545. 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 T deflation method inherently bypasses.

Instant research discovery