Tensor decomposition sampling algorithms for Bayesian inverse problems

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.

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.