Enhancing Tensor Network Algorithms for Many-Body Quantum Systems with MATLAB
By Dr. Anna Francuz, University of Vienna
“Using MATLAB, my team has developed a more stable and efficient approach for the differentiation of tensor network algorithms.… We are actively applying this workflow to deepen our understanding of strongly entangled quantum systems.”
A lodestar of quantum computing research has been achieving quantum advantage: demonstrating quantum computers that can outperform classical systems on practical tasks such as optimization, cryptography, or drug discovery. Before building large-scale quantum computers, however, it is necessary to run smaller scale tests and compare the results of those tests with theory. That work—with an emphasis on strongly entangled, many-body quantum systems—is the primary focus of my research at the University of Vienna.
Over the past 20 years or so, tensor networks have been shown to be among the best available tools for studying strongly entangled systems. This is due, in large part, to the way they reduce computational complexity (compared to quantum Monte Carlo methods, for example) by exploiting entanglement structure. Working with tensor networks, researchers can employ optimization techniques to find the ground state of a quantum system, from which they can calculate magnetization and other observables to understand the system’s properties and directly compare the results with experimental data.
Until recently, performing this optimization had been problematic because previous approaches relied on approximate calculations for the gradient, a crucial element in the differentiation process that drives optimization. In addition to significant errors, this drawback also frequently resulted in optimizations failing to find the global minimum. Researchers addressed some of these shortcomings by shifting to automatic differentiation (AD), but other issues emerged, such as divergences when differentiating a singular value decomposition (SVD) and fundamental problems with the way the AD was implemented in code.
Using MATLAB®, my team has developed a more stable and efficient approach for the differentiation of tensor network algorithms. We combine this novel AD implementation with two other key elements: tensor network contraction using a MATLAB community toolbox and high-performance computing using MATLAB Parallel Server™. We are actively applying this workflow to deepen our understanding of strongly entangled quantum systems.
MATLAB and Tensor Networks in Quantum Physics
By Temo Vekua, MathWorks
The application of tensor networks in quantum physics was pioneered by a preeminent physicist and leader in the field of quantum information theory, Dr. Ignacio Cirac from the Max Planck Institute of Quantum Optics in Munich. Two decades ago, Dr. Cirac and Dr. Frank Verstraete developed a theory of tensor networks to describe quantum many-body systems. Dr. Francuz traces her academic lineage back to Dr. Cirac: Her mentor, Dr. Norbert Schuch, was a postdoctoral researcher at Max Planck Institute, lecturer at the Technical University of Munich with Dr. Cirac, and later a tenured research group leader at the Max Planck Institute of Quantum Optics in Garching before becoming a professor at the University of Vienna.
Dr. Cirac relied on MATLAB for his earliest tensor network implementations. He and his team of approximately 40 researchers continue to use it today. “I mainly use MATLAB for scientific programming,” he explains. “MATLAB works well for me, and it can be as efficient as compiled languages when properly used.” The recognition of MATLAB as a vital tool in this field echoes the sentiments of Dr. Francuz, who notes, “I don’t think of myself as an exceptional programmer. Rather, I’m a physicist using MATLAB to do work with tensor networks that hasn’t been done before.”
Having shown the value of tensor networks for studying low-dimensional (1D and 2D), strongly correlated quantum systems, Dr. Cirac and his team are currently working to extend the method to 3D systems.
Tensor Network Contraction, Automatic Differentiation, and Optimization
When researchers in quantum physics employ tensor networks, we often work with translationally invariant infinite systems. Such systems extend indefinitely and have properties that are the same at every point when shifted (translated) by a certain distance This means that the Hamiltonian—an operator representing the total energy of the quantum system—exhibits either periodic repetition or uniformity throughout the system. So, whenever we calculate one of these observables, we must first contract the tensor network to something tractable (Figure 1). To do this in our workflow, we use ncon, a MATLAB function developed and shared with the community by Robert N. C. Pfeifer, Glen Evenbly, and their colleagues. The ncon function enables us to efficiently contract tensor networks and approximate the infinite environment of a local observable with a finite value.
Figure 1. A state |Ψ(A)⟩ in the square lattice shown in B (on the right) is generated by contracting the virtual (black) indices of a tensor in A (on the left). (Image credit: Francuz, A., Schuch, N., & Vanhecke, B. (2025). Stable and efficient differentiation of tensor network algorithms. Phys. Rev. Res., 7(1), 013237. American Physical Society. https://doi.org/10.1103/PhysRevResearch.7.013237. Licensed under CC BY 4.0.)
Once we have contracted the tensor network, which is represented as dlarray objects in MATLAB, the next step is to use AD to calculate the gradient of the energy, which is to be minimized when calculating the ground state. In implementing our AD approach in MATLAB, we addressed three critical issues with existing AD implementations: high memory usage, instability in gradient computations, and inaccuracies in backpropagating gradients.
The first issue was the high memory usage associated with AD. Generic AD approaches require substantial memory to store the intermediate objects generated during the iterative corner transfer matrix (CTM) procedure to compute local observables efficiently. (CTM is an iterative method used to approximate the contraction of two-dimensional tensor networks, particularly in projected entangled pair states, or PEPS.) While an alternative approach involving differentiation of the fixed-point equation can alleviate this memory issue, it introduces complications related to gauge fixing. Gauge fixing ensures consistency in tensor representation by resolving ambiguities, but its reliable implementation has been elusive. We identified the root cause of these difficulties and proposed two practical methods to reliably fix the gauge, thereby resolving the memory problem.
The second significant issue was the instability in the gradient computation of the SVD, a crucial component of the CTM algorithm. The gradient is poorly conditioned and becomes undefined in cases of degeneracies, where multiple singular values are equal. These degeneracies can destabilize gradient computations, as even small changes in the input can cause disproportionately large variations in the singular vectors. To address this, we uncovered and exploited a hidden gauge symmetry within the CTM tensor network. While this approach extended beyond the conventional scope of AD, it enabled us to develop a stable solution that remains compatible with AD frameworks.
Finally, we tackled inaccuracies in the previously used equations for the backpropagation of gradients through the truncated SVD. These equations incorrectly assumed that the truncated spectrum was zero, omitting an essential term and leading to significant errors. Backpropagating gradients, which involves computing how parameter changes propagate through a sequence of operations, requires precise gradient calculations for accurate optimization. We derived a corrected gradient formulation for the truncated SVD, accounting for the truncated spectrum. This correction greatly improved accuracy while requiring only minor modifications to existing AD code.
After implementing the AD improvements in MATLAB, we were ready to use the gradients generated by our AD implementation in optimization tasks, such as minimizing energy to find the ground state of quantum systems. For this part of the workflow, we use a limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm based on the fminunc function from Optimization Toolbox™, passing in the custom function we wrote to calculate the gradient (Figure 2).
Figure 2. Optimization plot showing two different gradients. The optimization quickly gets stuck by using the gradient from conventional AD (orange squares). The optimization converges to a global minimum by restarting the optimization using the gradient from the improved AD implementation (green circles). (Image credit: Image credit: Francuz, A., Schuch, N., & Vanhecke, B. (2025). Stable and efficient differentiation of tensor network algorithms. Phys. Rev. Res., 7(1), 013237. American Physical Society. https://doi.org/10.1103/PhysRevResearch.7.013237. Licensed under CC BY 4.0.)
Parallel Computing with MATLAB on the Austrian Scientific Computing
Phase diagrams are a crucial tool in the study of quantum many-body systems because they provide a visual representation of the different phases (states) of a system as a function of various parameters, including time. To create phase diagrams, we need to vary a parameter in the Hamiltonian over a certain range and calculate the ground state for each parameter value. This allows us to identify phase transitions when the system abruptly changes its behavior: for example, changing from a superconducting to an insulating state.
Because calculating the ground state for numerous parameter values sequentially on a single laptop can be a painfully slow process—and because the memory requirements are often beyond what standard laptops provide—we plan to perform this processing on the Austrian Scientific Computing (ASC) using MATLAB Parallel Server™. ASC is home to the most powerful supercomputers in Austria, including VSC-4, which has 790 nodes totaling 37,920 cores, and VSC-5, which has 710 nodes totaling 98,560 cores.
To get started, I attended a workshop conducted by MathWorks engineers covering parallel computing with MATLAB and using MATLAB on the HPC systems of the ASC. Running parallel ground state calculations on the ASC will enable our team to rapidly calculate properties of virtual quantum systems that we can directly compare with properties of real systems measured experimentally.
In fact, we are currently planning to use MATLAB Parallel Server on ASC to study phase diagrams for a programmable quantum simulator that can operate with 256 quantum bits (qubits). Exploring this system and others like it with 2D tensor networks is a focus of our ongoing research.
Published 2025
