Source-linked AI summary
A Polynomial-Scaling PDE Solver with Entanglement-Basis Tensor Networks
Abhijatmedhi Chotrattanapituk, Michael J. Landry, Chu-Liang Fu, Mingda Li
TL;DR
Conventional FEM realizations remain problem-dependent, while nonlinear analytic operators require additional variable substitution or iterative root finding in a direct coefficient representation. The method treats local element representations as independent degrees of freedom and imposes PDE, continuity, and boundary conditions through residuals optimized in an augmented coefficient space represented by tensor networks. With J, S, and network degree M bounded, storage remains polynomial in N; TN variational method (N = 10, B = 4, S = 1) differed by less than five percent of original sinusoidal amplitude.
Problem
Conventional FEM realizations remain problem-dependent, while nonlinear analytic operators require additional variable substitution or iterative root finding in a direct coefficient representation.
Source: §I · p. 1, B. Augmented coefficient space · p. 2Method
The method treats local element representations as independent degrees of freedom and imposes PDE, continuity, and boundary conditions through residuals optimized in an augmented coefficient space represented by tensor networks.
Source: §I · p. 1, §IX · p. 13Results
With J, S, and network degree M bounded, storage remains polynomial in N; TN variational method (N = 10, B = 4, S = 1) differed by less than five percent of original sinusoidal amplitude.
Source: C. Tensor Network Compression · p. 3, §VIII · p. 12Takeaways & Limitations
Reusable local tensor operators and explicit residuals let the framework accommodate different local approximation spaces and boundary relations within the same optimization framework.
Source: §I · p. 2, §IX · p. 13
Abstract
from arXiv · showhide
We develop a finite element method (FEM) for partial differential equation (PDE) solver based on the entanglement-basis representation introduced in our companion work. By lifting non-linear finite-element equations into an augmented coefficient space, the governing PDE together with boundary, initial, and inter-element constraints can be expressed through a unified quadratic residual minimization. Although this augmented space grows exponentially with the number of elements, its tensor-product structure allows it to be represented efficiently using tensor networks. Using the matrix product state (MPS) as a concrete example, we show that density matrix renormalization group (DMRG) sweeps enable element-by-element optimization without explicitly constructing the full augmented space. For bounded bond dimension, the resulting computational cost scales polynomially with the number of finite elements. We extend the framework to time-dependent problems through implicit temporal discretization and demonstrate convergence under both mesh and polynomial refinement using diffusion equations.
I. INTRODUCTION
The introduction positions this work as a practical PDE implementation of the companion TN-FE framework, reorganizing PDE, continuity, and boundary requirements as residuals while retaining local finite-element approximations. It also establishes notation and outlines how finite-element discretization turns analytic operator equations into algebraic coefficient problems.
I. INTRODUCTION
- The formulation imposes PDE, inter-element continuity, and boundary conditions through separate residual operators, allowing local approximation choices to be separated from element-coupling constraints.The same solver can accommodate different elemental representations and boundary relations, including analytic nonlinearities.Source: §I · pp. 1–2
II. NOTATIONS
- The notation uses boldface Roman symbols for non-scalar tensors, Greek characters for tensor indices, upright Roman labels, and brackets for iterated arrays.Repeated indices use Einstein summation; the text also assigns meanings to tensor product, contraction, direct sum, and block direct sum.Source: §II · p. 2
A. Finite Element Representation of Analytic Operator Equations
- The finite-element representation partitions the domain into elements with local basis functions, then Galerkin projection reduces the operator equation to an algebraic problem for the coefficient tensor.The test functions belong to the finite-element test space.Source: A. Finite Element Representation of… · p. 2
- An analytic operator can be expanded as a converging Fréchet-Taylor series, yielding a hierarchy of multilinear algebraic contributions in the weak form.The interaction kernels contract with the corresponding derivatives of the solution.Source: A. Finite Element Representation of… · p. 2
B. Augmented coefficient space
Nonlinear analytic operator equations motivate an augmented coefficient representation, while tensor-network compression addresses its exponential storage growth. When correlation order, bond dimension, and network degree remain bounded, storage scales polynomially with the number of elements.
B. Augmented coefficient space
- For nonlinear analytic operators, lifting the coefficient representation to Fock space represents correlations linearly, avoiding the need for variable substitution or iterative root finding.In the conventional representation, nonlinear operators require additional variable substitution and/or iterative root finding.Source: B. Augmented coefficient space · p. 2
- The augmented construction allocates J slots per element, each with B + 1 basis functions including an empty basis function, to represent correlations up to order J.The multilinear operator permits product terms to be reordered by element label.Source: B. Augmented coefficient space · p. 3
C. Tensor Network Compression
- Direct storage of the Fock-space coefficient tensor grows exponentially with the number of elements and correlation order, so tensor-network compression is needed.Source: C. Tensor Network Compression · p. 3
- Locality-based tensor-network choices include MPS for a one-dimensional chain, PEPS for a two-dimensional grid, and higher-order networks for correlations beyond pairs.Source: C. Tensor Network Compression · p. 3, Figure 2 · p. 3
- With bounded local correlation order, bond dimension, and network degree, tensor-network storage remains polynomial in the number of elements.This avoids the exponential growth of direct Fock-space storage.Source: C. Tensor Network Compression · p. 3
D. Variational solution of the discretized equation
The discretized problem is solved by minimizing a weighted residual over a tensor-network coefficient state, with single-site updates enabling convex subproblems. The formulation specializes analytic PDEs to operator equations and treats temporal dependence separately for the IVP demonstration.
D. Variational solution of the discretized equation
- The TN-FE solution minimizes a weighted residual over the tensor-network ansatz, whose global optimization is generally non-convex but whose single-site updates are convex weighted least-squares problems.The residual-space metric is positive-definite.Source: D. Variational solution of the… · pp. 3–4
E. Specialization to PDE
- Analytic PDEs can be reduced to operator-equation problems, and this work emphasizes explicitly constructing and implementing the PDE solver.When the nonlinearity is analytic, it can be expanded as a convergent series.Source: E. Specialization to PDE · p. 4
A. Target Equation
- The temporal variable is singled out so an MPS ansatz and DMRG can be used for the one-dimensional initial-value problem demonstration.The authors state that time could instead be treated as a dimension of q.Source: A. Target Equation · p. 4
B. Temporal discretization
The time-dependent problem is advanced iteratively with a finite-difference scheme, while reusable weak-form operators address repeated assembly across time steps. Spatial derivatives are approximated by elementwise projection tensors.
B. Temporal discretization
- The IVP is advanced across time steps using a finite-difference scheme, with implicit Euler chosen for simplicity; the rearranged formulation can be applied to any PDE.The method seeks the solution and its time derivatives at the next time step.Source: B. Temporal discretization · p. 4
C. Building Block Operators
- Recomputing G from scratch at every time step motivates precomputing the time-independent PDE terms into a weak form for connection to time-dependent terms.The interaction kernel may depend on time, while the other PDE terms remain the same.Source: C. Building Block Operators · p. 4
1. Derivative
- A spatial derivative is approximated on each element by a projection tensor obtained through the weak formulation.The derivative tensor carries components for different differentiation axes.Source: §1 · p. 5
2. Multiplication
Elementwise multiplication, constant factors, and boundary projections are represented as tensor-network projection operators. These building blocks support spatial operations and boundary coupling in the TN formulation.
2. Multiplication
- Spatial multiplication is approximated with a projection tensor.The tensor-network representation is shown in Fig. 3b.Source: §2 · p. 5
3. Constant Factor
- A constant factor acts on non-basis indices while leaving the basis action as an identity projection, supporting tensor-valued PDE functions.The constant-factor projection uses a Kronecker delta on the basis indices.Source: §3 · p. 5
4. Boundary Projection
- Projecting an element’s solution to a shared edge provides the boundary representation used to connect neighboring elements.The boundary is defined as the intersection of the elements, with its own basis functions.Source: §4 · p. 5, Figure 3 · p. 5
5. Summation
Tensor-network summation uses block direct sums when summands share a structure and an augmentation-and-filter construction when they do not. Shared components can be factored out so summation is localized to differing nodes.
5. Summation
- The formalism distinguishes summations whose summands share a TN structure from those whose structures differ.These two cases call for direct nodewise block sums and the augmentation-and-filter construction, respectively.Source: §5 · pp. 5–6
- For summands with the same TN structure, each node can be replaced by the block direct sum of corresponding summand nodes.For multiple bonding indices, the block form is required to contract the matrices to a scalar.Source: §5 · pp. 5–6, Figure 4 · p. 6
- When TN structures differ, augmenting each tensor with 1 before taking their tensor product provides a construction for extracting and adding their components.A filter extracts the components linear to A and B.Source: §5 · p. 6, Figure 5 · p. 7
- Factoring out shared TN components lets the summation be performed only at nodes that differ.When summands share a Fock-space coefficient tensor, compatible node tensors can be added before contracting with the common part.Source: §5 · pp. 6–7
6. Compound Network
The formulation keeps element coefficient tensors independent by enforcing compatibility through residual terms. It supplies compound-network operations and Fourier bases with closed-form tensors for element calculations.
6. Compound Network
- For compound operators, directly applying the weak formalism to the whole operator L can give a better approximation, at the cost of additional per-time-step calculations.The framework can construct higher-order derivatives, tensor-product variations, and tensor-calculus operators.Source: §6 · p. 7
1. Independent Local Approximation Spaces
- Additional residual terms enforce boundary conditions, continuity, and derivative compatibility, so coefficient tensors can remain independent rather than forming a globally conforming basis.The section explores bases with closed-form derivative, multiplication, and boundary-projection tensors.Source: §1 · p. 7
2. Fourier Basis
- Fourier basis derivatives can be matched to basis functions to construct differentiation matrices, while boundary projection reduces to evaluation at the element edge.The one-dimensional basis uses sine and cosine functions selected by the sign of the basis index; the boundary is represented by a scalar basis function.Source: §2 · p. 8
- Fourier multiplication tensors have a closed form; direct multiplication followed by basis-order truncation can be unstable because the basis functions are nonorthogonal.Element size and wavelength are independent, so terms above the cutoff can still contribute.Source: §2 · p. 8
3. Lagrange Polynomial Basis
The Lagrange polynomial basis supports derivative and multiplication tensors through polynomial operations, with Gauss-Lobatto-Legendre points offering quadrature and interpolation properties. Other bases can also be used with suitable projection and stability treatment.
3. Lagrange Polynomial Basis
- Lagrange basis derivatives can be rearranged as linear combinations of basis functions, with term matching guaranteed because the distinct-zero basis functions are independent and differentiation lowers polynomial order.Gauss-Lobatto-Legendre points support Gaussian quadrature exact through polynomial order 2B −3.Source: §3 · p. 8
- For Lagrange polynomials, multiplication tensors are obtained by evaluating products at basis zeroes, and boundary projection has the same form as for the Fourier basis.With GLL points, the basis also acts as an interpolation basis for functions up to polynomial order 2B −3.Source: §3 · pp. 8–9
4. Other Bases
- Other bases can be used if their derivative, multiplication, and boundary-projection tensors are calculated by projection to the original basis and any projection instability is addressed.Source: §4 · p. 9
V. UNIFIED RESIDUAL CONSTRUCTION
- Additional solution requirements can be added to the variational optimization target by framing each requirement as another variational problem.Source: §V · p. 9
A. Inter-element Continuity Residual
The formulation expresses inter-element compatibility and boundary conditions as residuals based on differences between boundary projections. It combines these with the PDE residual in a weighted variational objective, then uses an MPS ansatz for optimization.
A. Inter-element Continuity Residual
- Inter-element continuity residuals measure the difference between neighboring elements’ boundary projections, and can also enforce continuity of their pth derivatives.Source: A. Inter-element Continuity Residual · p. 9
B. Boundary Condition Residuals
- A boundary condition is treated as continuity between an element and an environment element whose boundary value is supplied by the problem.The residual substitutes the environment boundary value for one element’s boundary projection.Source: B. Boundary Condition Residuals · p. 9, Figure 7 · p. 9
C. Combined Variational Objective
- The variational target sums the individually applied weighted residual norms for the PDE, inter-element continuities, and boundary conditions.Minimizing this target gives a solution when one satisfying all continuity and boundary conditions exists.Source: C. Combined Variational Objective · p. 9, §VI · p. 10
VI. TARGET OPTIMIZATION
- Direct optimization over the Fock space is not viable because of its size; a tensor-network ansatz handles the structure but makes the optimization non-convex.The discussion focuses on one-dimensional problems with an MPS ansatz and DMRG.Source: §VI · p. 10
A. MPS Representation
- In the MPS representation, physical indices label local element coefficients, while bond dimension S controls retained correlations between adjacent elements.Source: A. MPS Representation · p. 10
B. Effective Single-site Problem
DMRG optimizes one MPS node at a time using an effective quadratic objective formed by contracting the fixed network. Recursively updated left and right environment contractions support sequential sweeps.
B. Effective Single-site Problem
- With all but one coefficient node fixed, the single-site update uses an effective quadratic objective and reduces to a convex weighted least-squares problem.The local update can be solved as a least-squares problem or an explicit matrix equation.Source: B. Effective Single-site Problem · pp. 10–11
C. Environment contractions
- The effective operator at an active site is formed by contracting the fixed tensors to its left and right, with contractions updated dynamically to avoid rebuilding them at every update.Source: C. Environment contractions · p. 11
D. DMRG algorithm
- DMRG optimizes MPS tensors sequentially back and forth along the chain, using prepared left and right contractions and repeated single-site updates.The sweeps can use a fixed number of rounds or be replaced by a stopping criterion.Source: D. DMRG algorithm · p. 11
VII. COMPLEXITY ANALYSIS
The complexity analysis tracks storage and runtime through element count, local basis size, bond dimension, and correlation order. The diffusion demonstration uses implicit Euler with DMRG sweeps and reports agreement with explicit Runge–Kutta alongside constraint preservation.
A. Parameter Count
- The storage analysis parameterizes the DMRG representation by N elements, basis size B, bond dimension S, and correlation order J.The passage also describes preparing H by contracting the PDE, continuity, and boundary terms through summation filters, using the same storage as H.Source: A. Parameter Count · p. 11
B. Runtime Complexity
- For PDEs whose residual terms are local to each element, the PDE contribution to runtime involves NJ possible terms; local continuity and boundary terms can reduce the other contribution.The text states that known structure in L can improve the bound and that locality allows terms to be added in local space.Source: B. Runtime Complexity · p. 12
VIII. DEMONSTRATION
- With implicit Euler time steps, DMRG sweeps minimize the weighted residual over an MPS representation for N uniform elements with B Lagrange basis functions and bond dimension S.The example uses a space-dependent diffusion equation and represents its IVP as a tensor circuit.Source: §VIII · p. 12, Figure 10 · p. 12
- For N = 10, B = 4, and S = 1, the TN variational solution differs from the explicit Runge–Kutta simulation by less than five percent of the sinusoid’s original amplitude, while continuity and smoothness hold within numerical precision.The authors interpret the constraint behavior as confirming that the variational target encodes the boundary constraints.Source: §VIII · p. 12
A. Convergence and Accuracy
The diffusion example shows convergence with spatial refinement up to a level set by the time step, and the conclusion reports controlled-error evolution with preserved constraints. The framework’s polynomial scaling depends on bounded tensor-network parameters and does not guarantee an advantage over conventional PDE solvers.
A. Convergence and Accuracy
- At fixed ∆t, error converges as basis order B or element count N increases, then plateaus until temporal resolution improves.The reported trends are shown for different time-step sizes in Fig. 13.Source: A. Convergence and Accuracy · p. 12, Figure 13 · p. 13
IX. CONCLUSION
- The variable-coefficient diffusion example reproduces expected time evolution with controlled error, preserves inter-element and boundary constraints, and converges under element and basis refinement.The convergence plot indicates that fixed temporal resolution limits further error reduction.Source: A. Convergence and Accuracy · p. 12, §IX · p. 13
- For fixed correlation order, bond dimension, and local basis size, storage and computational cost scale polynomially with the number of finite elements.The conclusion attributes access to the exponentially large augmented coefficient space to restricting the coefficient tensor to a tensor-network manifold.Source: §IX · p. 13
- Polynomial scaling is not a universal computational advantage: practical efficiency requires compact tensor-network representations, and the present implementation is limited to one-dimensional MPS/DMRG.Polynomial prefactors can grow substantially with local basis size, nonlinear correlation order, or required bond dimension.Source: §IX · pp. 13–14