Tensor product numerical methods for high-dimensional problems in probability and quantum calculations
Lead Research Organisation:
University of Bath
Department Name: Mathematical Sciences
Abstract
The tremendous complexity of contemporary technological processes induces a extremely high complexity in mathematical models that describe the physical laws of nature, but nonetheless computer simulations are indispensable for accurate quantitative predictions. Innovative models accounting for quantum effects or uncertain information are high-dimensional and hugely expensive - the lifetime of the Universe would not be enough to solve them by classical means. However, many such problems exhibit structure that if exploited can significantly reduce the computational effort. This project aims to break the complexity of high-dimensional models, and open them for routine use in computer simulation. This will be accomplished by the development of new methods that reveal hidden low-dimensional structures using classic mathematical methods, such as separation of variables and singular value decomposition. I will substantially extend the power of these methods, and apply them across a variety of important physical problems.
We know that our world is three-dimensional, so how do high-dimensional models arise? Imagine a bacterium in a lake. At any point in time it makes a move in an arbitrary direction. We cannot say definitely that the bacterium will reach a certain region in the lake and infect a plant growing there, but we can calculate the probability that this will happen. If we would like to consider all plants growing in a lake together we will have to store in computer memory the probability values for all plants. If we describe the behaviour of two bacteria at the same time, we will square the consumption of computer memory, since independent of the position of the first bacterium the second one still has the freedom to go anywhere. With an increasing number of bacteria, the amount of data explodes exponentially and quickly exhausts any reasonable memory. So the term "dimension" refers to the number of bacteria, and is naturally very high.
However, if bacteria move independently it is enough to store probability values just for one of them: the joint probability of the overall situation in the lake is then simply the product of the marginal probabilities. In reality, the bacteria will interact to some extent with near-by bacteria and will be influenced by the ambient flow in the lake. However, it is highly unlikely that a bacterium on one side of the lake will be affected by bacteria at the other. So the volume of data that effectively approximates the whole system will be many orders of magnitude smaller than the total number of degrees of freedom. The concept of high dimension in this example is in fact ubiquitous. An airplane, experiencing a fluctuating load on its wings, quantum effects unravelled by a magnetic resonance spectrometer, circadian rhythms or virus replication - all these phenomena are described by high-dimensional models. I will design methods that set new levels of prediction accuracy in these existing models, in particular such vital problems as subsurface flows of pollutants in an uncertain medium or stochastic virus dynamics, and extend the class of problems that can be tackled.
The new methods I am going to develop combine a range of mathematical techniques. The tensor product concept is a powerful data compression method, but it is in fact nothing more than an immediate generalisation of the separation of variable concept for continuous functions that can be accurately computed via the singular value decomposition. The workhorses in the tensor product framework are alternating optimisation algorithms. Their convergence can be substantially enhanced by employing ideas from classical Krylov iterative methods for matrix equations. I will extend tensor product methods further, and embody them in publicly available software with transparent user interfaces to popular scientific packages in order to encourage other researchers to try the new methodology in real-life problems.
We know that our world is three-dimensional, so how do high-dimensional models arise? Imagine a bacterium in a lake. At any point in time it makes a move in an arbitrary direction. We cannot say definitely that the bacterium will reach a certain region in the lake and infect a plant growing there, but we can calculate the probability that this will happen. If we would like to consider all plants growing in a lake together we will have to store in computer memory the probability values for all plants. If we describe the behaviour of two bacteria at the same time, we will square the consumption of computer memory, since independent of the position of the first bacterium the second one still has the freedom to go anywhere. With an increasing number of bacteria, the amount of data explodes exponentially and quickly exhausts any reasonable memory. So the term "dimension" refers to the number of bacteria, and is naturally very high.
However, if bacteria move independently it is enough to store probability values just for one of them: the joint probability of the overall situation in the lake is then simply the product of the marginal probabilities. In reality, the bacteria will interact to some extent with near-by bacteria and will be influenced by the ambient flow in the lake. However, it is highly unlikely that a bacterium on one side of the lake will be affected by bacteria at the other. So the volume of data that effectively approximates the whole system will be many orders of magnitude smaller than the total number of degrees of freedom. The concept of high dimension in this example is in fact ubiquitous. An airplane, experiencing a fluctuating load on its wings, quantum effects unravelled by a magnetic resonance spectrometer, circadian rhythms or virus replication - all these phenomena are described by high-dimensional models. I will design methods that set new levels of prediction accuracy in these existing models, in particular such vital problems as subsurface flows of pollutants in an uncertain medium or stochastic virus dynamics, and extend the class of problems that can be tackled.
The new methods I am going to develop combine a range of mathematical techniques. The tensor product concept is a powerful data compression method, but it is in fact nothing more than an immediate generalisation of the separation of variable concept for continuous functions that can be accurately computed via the singular value decomposition. The workhorses in the tensor product framework are alternating optimisation algorithms. Their convergence can be substantially enhanced by employing ideas from classical Krylov iterative methods for matrix equations. I will extend tensor product methods further, and embody them in publicly available software with transparent user interfaces to popular scientific packages in order to encourage other researchers to try the new methodology in real-life problems.
Planned Impact
* Academic dissemination
Tensor product methods can be interesting for different beneficiaries in both applied mathematics and in physics, chemistry and biology, which are closer to real problems. Efficient data compression techniques based on separation of variables could become a method of choice for complex stochastic models that otherwise require too much computational effort with existing methods. To name just a few examples:
- Monte Carlo-based methods for uncertainty quantification may require the solution of thousands or even millions of high-resolution deterministic problems, and take days or months to complete.
- Stochastic models of virus replication in vaccine design have been avoided for a long time due to their high complexity, despite the fact that replication of viruses is an intrinsically stochastic process.
- State of the art model reduction methods in magnetic resonance spectroscopy do make it possible to simulate NMR experiments for proteins, but they still consume up to a terabyte of shared memory. To compare: our preliminary tensor product calculations for the pulse-acquire model (one of the simplest experiments) needed a couple of gigabytes, and were tractable on a laptop.
- Leading Density Functional Theory software packages are still prohibitively expensive, and may take weeks to simulate a point defect in a crystal with extremely simplified boundary conditions.
New numerical methods can reduce the amount of high-load scientific computations, and by that cut down the energy consumption of supercomputers. Besides, if new accurate models are available for numerical simulations, we will have more chances to discover some new physical phenomena, and broaden our understanding of nature.
To disseminate my work I will publish mathematical aspects, descriptions of algorithms and their analysis in top rated mathematical journals (SIAM, LAA), complemented by papers in physical journals (JCP, CPC, Phys. Review), to present results of the application of tensor product methods to particular problems in collaboration with experts in these areas. To provide equal opportunity to access to these publications, I will accompany journal articles with up to date preprint versions. I will also continue to give talks at major international and UK conferences to engage other researchers with my project. I plan to attend several established annual conferences in applied mathematics (see Justification of resources for details), as well as some events in application areas, which will be identified together with my beneficiaries (for example, I was invited to the 9th European Conference on Mathematical and Theoretical Biology). I will be also glad to visit local workshops and seminars
organised by collaborators.
* Software development and user support
Numerical algorithms are often only as good as their particular software implementation. The main reason why tensor approximation methods based on the singular value decomposition have quickly become preferred to generic error minimisation techniques (such as the Newton method) is the incredible efficiency and reliability of the SVD algorithm realised in LAPACK and MKL libraries. I am going to take into account my previous experience and the needs in applications to develop a versatile multi-level tensor product software package which will be made openly available.
To encourage beneficiaries to try the new methodology in their applications, I will implement a MATLAB interface to tensor product algorithms, with as transparent overloading as possible. It will be complemented by high-performance low-level codes for computationally demanding problems. The software will be accompanied by an extensive documentation and typical benchmark and tutorial examples. An open-source toolkit with quick start possibilities will make an instant impact on recognition of tensor product methods in a broad community and I see the biggest impact of this project to come from this.
Tensor product methods can be interesting for different beneficiaries in both applied mathematics and in physics, chemistry and biology, which are closer to real problems. Efficient data compression techniques based on separation of variables could become a method of choice for complex stochastic models that otherwise require too much computational effort with existing methods. To name just a few examples:
- Monte Carlo-based methods for uncertainty quantification may require the solution of thousands or even millions of high-resolution deterministic problems, and take days or months to complete.
- Stochastic models of virus replication in vaccine design have been avoided for a long time due to their high complexity, despite the fact that replication of viruses is an intrinsically stochastic process.
- State of the art model reduction methods in magnetic resonance spectroscopy do make it possible to simulate NMR experiments for proteins, but they still consume up to a terabyte of shared memory. To compare: our preliminary tensor product calculations for the pulse-acquire model (one of the simplest experiments) needed a couple of gigabytes, and were tractable on a laptop.
- Leading Density Functional Theory software packages are still prohibitively expensive, and may take weeks to simulate a point defect in a crystal with extremely simplified boundary conditions.
New numerical methods can reduce the amount of high-load scientific computations, and by that cut down the energy consumption of supercomputers. Besides, if new accurate models are available for numerical simulations, we will have more chances to discover some new physical phenomena, and broaden our understanding of nature.
To disseminate my work I will publish mathematical aspects, descriptions of algorithms and their analysis in top rated mathematical journals (SIAM, LAA), complemented by papers in physical journals (JCP, CPC, Phys. Review), to present results of the application of tensor product methods to particular problems in collaboration with experts in these areas. To provide equal opportunity to access to these publications, I will accompany journal articles with up to date preprint versions. I will also continue to give talks at major international and UK conferences to engage other researchers with my project. I plan to attend several established annual conferences in applied mathematics (see Justification of resources for details), as well as some events in application areas, which will be identified together with my beneficiaries (for example, I was invited to the 9th European Conference on Mathematical and Theoretical Biology). I will be also glad to visit local workshops and seminars
organised by collaborators.
* Software development and user support
Numerical algorithms are often only as good as their particular software implementation. The main reason why tensor approximation methods based on the singular value decomposition have quickly become preferred to generic error minimisation techniques (such as the Newton method) is the incredible efficiency and reliability of the SVD algorithm realised in LAPACK and MKL libraries. I am going to take into account my previous experience and the needs in applications to develop a versatile multi-level tensor product software package which will be made openly available.
To encourage beneficiaries to try the new methodology in their applications, I will implement a MATLAB interface to tensor product algorithms, with as transparent overloading as possible. It will be complemented by high-performance low-level codes for computationally demanding problems. The software will be accompanied by an extensive documentation and typical benchmark and tutorial examples. An open-source toolkit with quick start possibilities will make an instant impact on recognition of tensor product methods in a broad community and I see the biggest impact of this project to come from this.
Organisations
People |
ORCID iD |
| Sergey Dolgov (Principal Investigator / Fellow) |
Publications
Benner P
(2016)
Low-rank solvers for unsteady Stokes-Brinkman optimal control problem with random data
in Computer Methods in Applied Mechanics and Engineering
Benner P
(2017)
Fast iterative solution of the Bethe-Salpeter eigenvalue problem using low-rank and QTT tensor approximation
in Journal of Computational Physics
Dolgov S
(2019)
Approximation and sampling of multivariate probability distributions in the tensor train decomposition
in Statistics and Computing
Dolgov S
(2021)
Tensor Decomposition Methods for High-dimensional Hamilton--Jacobi--Bellman Equations
in SIAM Journal on Scientific Computing
Dolgov S
(2021)
Guaranteed a posteriori error bounds for low-rank tensor approximate solutions
in IMA Journal of Numerical Analysis
Dolgov S
(2020)
Parallel cross interpolation for high-precision calculation of high-dimensional integrals
in Computer Physics Communications
Dolgov S
(2019)
Preconditioners and Tensor Product Solvers for Optimal Control Problems from Chemotaxis
in SIAM Journal on Scientific Computing
Dolgov S
(2019)
A Hybrid Alternating Least Squares--TT-Cross Algorithm for Parametric PDEs
in SIAM/ASA Journal on Uncertainty Quantification
Dolgov S
(2017)
Low-Rank Solution to an Optimization Problem Constrained by the Navier--Stokes Equations
in SIAM Journal on Scientific Computing
Dolgov S
(2018)
Direct tensor-product solution of one-dimensional elliptic equations with parameter-dependent coefficients
in Mathematics and Computers in Simulation
| Description | Reduced order modelling is a powerful technique for speeding up simulation of large dynamical systems, such as non-stationary partial differential equations with uncertain data. Tensor product algorithms allow to compute a reduced model without ever solving the full original problem. However, an optimal control of such systems is more challenging. It requires solving optimality equations, which led to instabilities in existing tensor algorithms. In the framework of the project, I developed a new technique that mitigates this issue by reducing forward, adjoint and gradient equations separately. It allowed fast modelling of inverse stochastic problems in fluid dynamics (http://dx.doi.org/10.1137/15M1040414 and http://dx.doi.org/10.1016/j.cma.2016.02.004). Another difficulty stems from the problems with high variance of coefficients, and hence a larger storage needed by the reduced tensor model. Efficient simulation of such scenarios requires more sophisticated algorithms that exploit as much of the problem structure as possible. The new method I developed in the second year addresses this issue by combining the classical orthogonal projection and the optimised sampling of the random variables associated with the uncertainty (http://dx.doi.org/10.1137/17M1138881). This has drastically improved the performance of the tensor product approach in solving stochastic PDE models with highly varying coefficients and made it more efficient than traditional sampling techniques. The tensorised stochastic PDE solution turned out to be an accurate surrogate model for the Bayesian approach to the inverse problems, which I started investigating in the third year. Moreover, I developed a completely new paradigm of combining tensor approximations and Monte Carlo sampling. A tensor product approximation of a posterior density function allows fast high-dimensional integration, and hence uncorrelated sampling from the approximate density via the conditional distribution method (https://doi.org/10.1007/s11222-019-09910-z). This method is quite general, as it can work in both high- and low-accuracy regimes. An accurate approximation allows the almost-exact samples to be used directly as an optimal quadrature for posterior statistics, while even a sketch of a function learned by the tensor interpolation may generate independence proposals for Markov Chain Monte Carlo or Importance Sampling, producing low correlation of samples or small variance of importance weights. The conditional distribution approach is also embarrassingly parallel over samples, in contrast to classical MCMC. After the project was completed, the conditional distribution method was extended to the Deep Inverse Rosenblatt Transport (DIRT), an adaptive mesh generator in high dimensions. This allows approximation of complicated and concentrated distributions by an iterative refinement of the coordinates using a sequence of tempered distributions. Moreover, any hierarchy in the model (such as different resolutions of the discretisation) can be incorporated seamlessly at different DIRT layers to speed up the construction of intermediate layers, while ensuring a high accuracy of the model by switching to the high resolution at the final level. Further high performance developments include implementation of a greedy low-rank solution algorithm for matrix equations on GPU (https://doi.org/10.1186/s13408-019-0077-0) and a hybrid parallel tensor interpolation algorithm (manuscript in preparation). High arithmetic intensity of dense tensor algebra allows efficient shared memory parallelisation. Using the greedy algorithm, we were able to reconstruct the whole brain connectivity map from viral tracing experiments in a few hours on a single GPU card, which was inaccessible previously due to colossal computing costs. Parallel tensor interpolation algorithm (https://doi.org/10.1016/j.cpc.2019.106869), scaled up to 512 cores, enabled calculation of Ising integrals in spontaneous ferromagnetic magnetisation with high precision (up to 100 decimal digits). Another developments in the end of the project was a tensor decomposition algorithm for the stationary Hamilton-Jacobi-Bellman (HJB) equations (https://arxiv.org/pdf/1908.01533.pdf, accepted in SISC). This can break the curse of dimensionality in dynamic programming for stochastic nonlinear PDEs (for example the Allen-Cahn mixture separation model), and synthesize a much more efficient feedback control signal compared to the traditionally used Linear-Quadratic Regulator. This work will be extended to time-dependent (finite horizon) HJB equations, and to higher dimensions by developing tailored Model Order Reduction techniques. New directions emerged in the New Horizons project are optimal linearisation and reduction of unstable systems. |
| Exploitation Route | Fast computation of an optimal control under a lack of knowledge should allow researchers and engineers to investigate more scenarios and designs of a technological process in a timely manner. Moreover, the developed algorithms reduce the computer memory footprint, and allow to carry out on a desktop (or even on a laptop) simulations that were previously only accessible with a supercomputer. The new hybrid algorithm is implemented in a publicly available Matlab package (http://people.bath.ac.uk/sd901/als-cross-algorithm/) that can be used by other researchers to speed up their simulation of models with uncertainty. The conditional distribution package (https://github.com/dolgov/TT-IRT) can be potentially used in any statistical application, significantly improving prediction accuracy or computational costs. Higher efficiency is attained of course for the probability density functions that admit more accurate tensor product approximations. These include smooth functions and weakly correlated systems, where correlations between random variables decay with the distance between indices of these variables (this might require reordering of variables). Hierarchical, cascade and decaying models (e.g. PDEs with Karhunen-Loeve parametrisation of random coefficients, reliability estimation, Bayesian networks) have proved to be particularly suitable for this approach. On the other hand, the sampling nature of the method imposes no restriction on statistics of interest. The Version 2.0 of this package has been extended to the Deep Inverse Rosenblatt Transport algorithms, which can approximate much wider classes of density functions, included concentrated densities emerging in Bayesian inverse problems with low noise and many measurements. The toolbox is equipped with automatic installation scripts to simplify the uptake by end users, and an extensive documentation on both the global structure and individual functions. The code is written in a simple procedural way with bottleneck operations optimized for efficiency. A more object-oriented version of the Functional Tensor Train and Deep Inverse Rosenblatt Transport toolbox (https://github.com/fastfins/ftt.m) that was developed independently by T. Cui has followed. Tensor methods for the feedback control using the Hamilton-Jacobi-Bellman (HJB) formalism have been implemented in another package https://github.com/dolgov/TT-HJB. Again, the construction and solution of the HJB equations is automated in a dedicated procedure, which requires only the dynamical system description from the user. |
| Sectors | Digital/Communication/Information Technologies (including Software),Environment,Healthcare,Manufacturing, including Industrial Biotechology |
| URL | https://people.bath.ac.uk/sd901/research/software/ |
| Description | Research outputs from this project have been taken up by Terra Quantum AG (https://terraquantum.swiss/) who has signed a year-long consultancy agreement with the PI on the adaptation and development of tensor method software for applied problems brought by clients of Terra Quantum (precise details are currently under NDA). This project was extended for another year. |
| First Year Of Impact | 2021 |
| Sector | Digital/Communication/Information Technologies (including Software) |
| Impact Types | Economic |
| Description | Overcoming the curse of dimensionality in dynamic programming by tensor decompositions |
| Amount | £202,430 (GBP) |
| Funding ID | EP/V04771X/1 |
| Organisation | Engineering and Physical Sciences Research Council (EPSRC) |
| Sector | Public |
| Country | United Kingdom |
| Start | 01/2021 |
| End | 01/2023 |
| Description | Tensor decomposition sampling algorithms for Bayesian inverse problems |
| Amount | £150,426 (GBP) |
| Funding ID | EP/T031255/1 |
| Organisation | Engineering and Physical Sciences Research Council (EPSRC) |
| Sector | Public |
| Country | United Kingdom |
| Start | 01/2021 |
| End | 12/2023 |
| Description | University of Sydney Mathematical Research Institute International visitor program |
| Amount | $6,500 (AUD) |
| Organisation | University of Sydney |
| Sector | Academic/University |
| Country | Australia |
| Start | 06/2019 |
| End | 07/2019 |
| Title | ALS-Cross algorithm |
| Description | This is a small Matlab toolkit for Uncertainty Quantification using the tensor decompositions. The main algorithm is a combination of the ALS and TT Cross methods. The package contains high-level drivers for setting and solving high-dimensional linear systems with block-diagonal matrices arising from stochastic and parametric PDEs. Moreover, routines for comparing different approaches (tensor decomposition, sparse collocation, quasi Monte Carlo) are provided. |
| Type Of Technology | Software |
| Year Produced | 2017 |
| Open Source License? | Yes |
| Impact | - |
| URL | http://people.bath.ac.uk/sd901/als-cross-algorithm/ |
| Title | Low Rank Connectome |
| Description | This package implements the greedy low-rank algorithm for general matrix equations (see https://arxiv.org/abs/1808.05510), applied to the completion of a spatial connectome of a brain from a few viral tracing measurements. The software is optimised for CUDA GPU acceleration in Matlab. This enables inference of the whole flatmap connectivity in a couple of hours on a single Nvidia card, in contrast to days of computing with the gradient descent optimisation. |
| Type Of Technology | Software |
| Year Produced | 2018 |
| Open Source License? | Yes |
| Impact | - |
| URL | https://gitlab.mpi-magdeburg.mpg.de/kuerschner/lowrank_connectome |
| Title | TT-HJB |
| Description | Matlab package to form and solve Hamilton-Jacobi-Bellman (HJB) equations for the optimal feedback control synthesis for stochastic dynamical systems. |
| Type Of Technology | Software |
| Year Produced | 2020 |
| Open Source License? | Yes |
| Impact | No **notable** impacts traced yet, but we have received a few requests from PhD students working on similar problems who would like to build on our code. |
| URL | https://github.com/dolgov/TT-HJB |
| Title | TT-IRT |
| Description | This toolkit implements the Tensor Train Cross approximation of a multivariate joint probability density function and the Inverse Rosenblatt Transform (Conditional Distribution) + MCMC sampling of that distribution using its Tensor Train approximation, see https://arxiv.org/abs/1810.01212 for details. The conditional distribution sampling algorithm is implemented in C and linked to high-level interfaces in Matlab and Python. The package is equipped with two examples from reliability estimation (lifetime of a shock absorber) and PDE-constrained Bayesian inverse problems (elliptic equation with random coefficients). |
| Type Of Technology | Software |
| Year Produced | 2018 |
| Open Source License? | Yes |
| Impact | - |
| URL | https://github.com/dolgov/TT-IRT |