# Tensor decomposition sampling algorithms for Bayesian inverse problems

Lead Research Organisation:
University of Bath

Department Name: Mathematical Sciences

### Abstract

Understanding the cause of events we observe is a vital task in both everyday life and most pressing global challenges.

Was the rainy summer really due to the polar ice meltdown?

Will an aircraft wing made from lightweight composite materials break?

Or will a nuclear waste leak from the storage and reach a drinking water well?

Often, a simple yes/no answer is impossible.

Virtually all statistical forecasts operate with the probability of an event to happen,

that is a quantitative characteristics of the uncertainty in the prediction.

However, such predictions can be very poor if only some observations are given, and nothing is known about the underlying natural processes.

On the other extreme, even a super-accurate model would be meaningless if there is no data to initialise it.

How can we predict weather for tomorrow if we have collected no measurements today?

In practice we usually have both model and data, but of limited quality:

a partially inaccurate model, and a partially incomplete and noisy observation data.

How can we produce a forecast that is best in some sense, together with its uncertainty?

A mathematically rigorous answer to this question is known for centuries: the Bayes theorem.

However, it might be extremely challenging to employ it in practice due to the so-called curse of dimensionality.

The Bayes theorem describes the answer in the form of a joint probability distribution function that depends on all tunable parameters of the model.

Although a forecast of interest can be just a single number,

computing this number requires numerical integration of the probability function.

Straightforward attempt to do so involves computing probability values for all possible combinations of the parameters.

This renders the amount of computations growing exponentially with the dimensionality, that is the number of parameters, in the problem.

While some simple case with only one parameter might be calculable in milliseconds,

for high-dimensional problems with tens of parameters even the lifetime of the Universe could be not enough to solve them straightforwardly.

However, many probability functions arising in the Bayesian approach contain hidden structure that may aid computational methods significantly.

This project aims to reveal and exploit this structure to make Bayesian statistical predictions computationally tractable.

I will approach this by developing new algorithms that combine advantages of several classical mathematical methods.

The core of the project is the tensor product decompositions.

This is a powerful family of methods for data compression that originate from the simple separation of variables.

The efficiency of tensor decompositions relies on assumption that the model parameters are weakly dependent in a certain sense (for example, the first parameter has little influence on the last one).

Another classical method from probability, the Rosenblatt transformation, will be exploited

to develop an adaptive procedure to compute a change of coordinates that fulfils the assumption of weak dependence for the transformed parameters.

The new methods will enable better predictions driven by Bayes-optimal statistical analysis in complicated inverse problems such as those arising in testing and certification of new composite materials for aerospace industry.

Moreover, embodying the new algorithms in open-source software in collaboration with academic groups in statistics and engineering will pave the way to even wider uptake of the proposed methodology for treatment of uncertainty.

Was the rainy summer really due to the polar ice meltdown?

Will an aircraft wing made from lightweight composite materials break?

Or will a nuclear waste leak from the storage and reach a drinking water well?

Often, a simple yes/no answer is impossible.

Virtually all statistical forecasts operate with the probability of an event to happen,

that is a quantitative characteristics of the uncertainty in the prediction.

However, such predictions can be very poor if only some observations are given, and nothing is known about the underlying natural processes.

On the other extreme, even a super-accurate model would be meaningless if there is no data to initialise it.

How can we predict weather for tomorrow if we have collected no measurements today?

In practice we usually have both model and data, but of limited quality:

a partially inaccurate model, and a partially incomplete and noisy observation data.

How can we produce a forecast that is best in some sense, together with its uncertainty?

A mathematically rigorous answer to this question is known for centuries: the Bayes theorem.

However, it might be extremely challenging to employ it in practice due to the so-called curse of dimensionality.

The Bayes theorem describes the answer in the form of a joint probability distribution function that depends on all tunable parameters of the model.

Although a forecast of interest can be just a single number,

computing this number requires numerical integration of the probability function.

Straightforward attempt to do so involves computing probability values for all possible combinations of the parameters.

This renders the amount of computations growing exponentially with the dimensionality, that is the number of parameters, in the problem.

While some simple case with only one parameter might be calculable in milliseconds,

for high-dimensional problems with tens of parameters even the lifetime of the Universe could be not enough to solve them straightforwardly.

However, many probability functions arising in the Bayesian approach contain hidden structure that may aid computational methods significantly.

This project aims to reveal and exploit this structure to make Bayesian statistical predictions computationally tractable.

I will approach this by developing new algorithms that combine advantages of several classical mathematical methods.

The core of the project is the tensor product decompositions.

This is a powerful family of methods for data compression that originate from the simple separation of variables.

The efficiency of tensor decompositions relies on assumption that the model parameters are weakly dependent in a certain sense (for example, the first parameter has little influence on the last one).

Another classical method from probability, the Rosenblatt transformation, will be exploited

to develop an adaptive procedure to compute a change of coordinates that fulfils the assumption of weak dependence for the transformed parameters.

The new methods will enable better predictions driven by Bayes-optimal statistical analysis in complicated inverse problems such as those arising in testing and certification of new composite materials for aerospace industry.

Moreover, embodying the new algorithms in open-source software in collaboration with academic groups in statistics and engineering will pave the way to even wider uptake of the proposed methodology for treatment of uncertainty.

### Planned Impact

The most immediate impact will be new computational algorithms for function approximations, equipped with numerical analysis and efficient open-source implementation.

Tensor product methods have already reached a considerable state of maturity, although their uptake in probability and statistics was somehow limited.

A combination with the adaptive transformation of variables, provided by the Rosenblatt transformation,

will significantly expand the class of tractable high-dimensional functions.

This knowledge will be of interest to the international community of researchers in mathematics and statistics, as well as more applied areas such as physics, chemistry and engineering, where high-dimensional problems become more and more ubiquitous.

The mathematical theory of the algorithms established in the course of the project will be the starting point for further developments of function approximation methods in applied mathematics, while the publicly available implementation will simplify the practical computations using the new techniques in applications.

Rigorous quantification of uncertainty in computer modelling is often dismissed in existing engineering practice (for example in materials science)

due to the lack of robust numerical methods.

Traditional Monte Carlo-based methods may involve solution of tens-hundreds thousands high-resolution deterministic problems, and take days or months to complete.

The project strives at producing new efficient computational tools for real-life UQ problems,

with the strong potential of replacing currently used techniques (such as random walk MCMC).

In particularly challenging cases, such as the inverse stochastic elasticity problem,

these tools will essentially become an enabling technology, allowing optimal and certified statistical modelling of e.g. manufacturing defects in new composite materials and consequences of those.

A major reduction in testing requirements and hence in the development time and costs, will encourage the development of new materials and accelerate their commercialisation and adoption in industry.

Moreover, faster numerical methods can reduce the amount of high-load scientific computations, and by that cut down the energy consumption of supercomputers and the associated carbon footprint.

Tensor product methods have already reached a considerable state of maturity, although their uptake in probability and statistics was somehow limited.

A combination with the adaptive transformation of variables, provided by the Rosenblatt transformation,

will significantly expand the class of tractable high-dimensional functions.

This knowledge will be of interest to the international community of researchers in mathematics and statistics, as well as more applied areas such as physics, chemistry and engineering, where high-dimensional problems become more and more ubiquitous.

The mathematical theory of the algorithms established in the course of the project will be the starting point for further developments of function approximation methods in applied mathematics, while the publicly available implementation will simplify the practical computations using the new techniques in applications.

Rigorous quantification of uncertainty in computer modelling is often dismissed in existing engineering practice (for example in materials science)

due to the lack of robust numerical methods.

Traditional Monte Carlo-based methods may involve solution of tens-hundreds thousands high-resolution deterministic problems, and take days or months to complete.

The project strives at producing new efficient computational tools for real-life UQ problems,

with the strong potential of replacing currently used techniques (such as random walk MCMC).

In particularly challenging cases, such as the inverse stochastic elasticity problem,

these tools will essentially become an enabling technology, allowing optimal and certified statistical modelling of e.g. manufacturing defects in new composite materials and consequences of those.

A major reduction in testing requirements and hence in the development time and costs, will encourage the development of new materials and accelerate their commercialisation and adoption in industry.

Moreover, faster numerical methods can reduce the amount of high-load scientific computations, and by that cut down the energy consumption of supercomputers and the associated carbon footprint.

### Organisations

## People |
## ORCID iD |

Sergey Dolgov (Principal Investigator) |

### Publications

Antil H
(2022)

*TTRISK: Tensor train decomposition algorithm for risk averse optimization*in Numerical Linear Algebra with Applications
Cui T
(2024)

*Deep Importance Sampling Using Tensor Trains with Application to a Priori and a Posteriori Rare Events*in SIAM Journal on Scientific Computing
Cui T
(2021)

*Deep Composition of Tensor-Trains Using Squared Inverse Rosenblatt Transports*in Foundations of Computational Mathematics
Cui T
(2023)

*Scalable conditional deep inverse Rosenblatt transports using tensor trains and gradient-based dimension reduction*in Journal of Computational Physics
Dolgov S
(2024)

*Tensor product approach to modelling epidemics on networks*in Applied Mathematics and Computation
Gourianov N
(2021)

*A Quantum Inspired Approach to Exploit Turbulence Structures*
Rohrbach P
(2022)

*Rank Bounds for Approximating Gaussian Densities in the Tensor-Train Format*in SIAM/ASA Journal on Uncertainty QuantificationDescription | This project concerns the development of fast computational algorithms for solving inverse problems in Bayesian way. The ultimate solution to those is given by a multivariate posterior probability density function, which may concentrate along complicated manifolds in a high-dimensional space if the data is not fully informative. We developed a method to approximate this function by a sequence of coordinate transformations driven by simpler density functions, such as different temperings of the final function. Each function is approximated by a tensor train decomposition, which allows both fast construction and evaluation of the coordinate transformation needed in the next step. The new method has demonstrated higher accuracy-complexity efficiency than a neural network coordinate transformation sampler in some examples of inverse and calibration problems constrained by differential equations, such as epidemiological models, subsurface flows and elasticity. Further strands started are risk-averse optimisation and estimation of rare event probabilities. |

Exploitation Route | Two software packages implementing the main coordinate transformation method (https://github.com/dolgov/TT-IRT and https://github.com/fastfins/ftt.m) are publicly available for use by a wide range of researchers and engineers in statistics, inverse problems and data science. |

Sectors | Digital/Communication/Information Technologies (including Software) |

URL | https://github.com/dolgov/TT-IRT |

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 |

Title | TTRISK |

Description | TTRISK: Tensor Train Algorithm for Risk Averse Optimization. Approximates the Conditional Value at Risk (CVaR) by using Tensor Train approximations and adaptively smoothed indicator functions to compute expectations over random variables in CVaR. |

Type Of Technology | Software |

Year Produced | 2022 |

Open Source License? | Yes |

Impact | not yet |

URL | https://github.com/dolgov/TTRISK |