Search

The Kullback-Leibler (KL) divergence is frequently used in data science. For discrete distributions on large state spaces, approximations of probability vectors may result in a few small negative entries, rendering the KL divergence undefined. We address this problem by introducing a parameterized family of substitute divergence measures, the shifted KL (sKL) divergence measures. Our approach is generic and does not increase the computational overhead. We show that the sKL divergence shares important theoretical properties with the KL divergence and discuss how its shift parameters should be chosen. If Gaussian noise is added to a probability vector, we prove that the average sKL divergence converges to the KL divergence for small enough noise. We also show that our method solves the problem of negative entries in an application from computational oncology, the optimization of Mutual Hazard Networks for cancer progression using tensor-train approximations., Comment: 11 pages, 6 figures, minor clarifications, published version

The quality of the inverse approach in electroencephalography (EEG) source analysis is - among other things - depending on the accuracy of the forward modeling approach, i.e., the simulation of the electric potential for a known dipole source in the brain. Here, we use multilayer sphere modeling scenarios to investigate the performance of three different finite element method (FEM) based EEG forward approaches - subtraction, Venant and partial integration - in the presence of tissue conductivity anisotropy in the source space. In our studies, the effect of anisotropy on the potential is related to model errors when ignoring anisotropy and to numerical errors, convergence behavior and computational speed of the different FEM approaches. Three different source space anisotropy models that best represent adult, child and premature baby volume conduction scenarios, are used. Major findings of the study include (1) source space conductivity anisotropy has a significant effect on electric potential computation: The effect increases with increasing anisotropy ratio; (2) with numerical errors far below anisotropy effects, all three FEM approaches are able to model source space anisotropy accordingly, with the Venant approach offering the best compromise between accuracy and computational speed; (3) FE meshes have to be fine enough in the subdomain between the source and the sensors that capture its main activity. We conclude that, especially for the analysis of cortical development, but also for more general applications using EEG source analysis techniques, source space conductivity anisotropy should be modeled and the FEM Venant approach is an appropriate method., Comment: 20 pages, 4 figures, 5 tables

Motivation: We consider continuous-time Markov chains that describe the stochastic evolution of a dynamical system by a transition-rate matrix $Q$ which depends on a parameter $\theta$. Computing the probability distribution over states at time $t$ requires the matrix exponential $\exp(tQ)$, and inferring $\theta$ from data requires its derivative $\partial\exp\!(tQ)/\partial\theta$. Both are challenging to compute when the state space and hence the size of $Q$ is huge. This can happen when the state space consists of all combinations of the values of several interacting discrete variables. Often it is even impossible to store $Q$. However, when $Q$ can be written as a sum of tensor products, computing $\exp(tQ)$ becomes feasible by the uniformization method, which does not require explicit storage of $Q$. Results: Here we provide an analogous algorithm for computing $\partial\exp\!(tQ)/\partial\theta$, the differentiated uniformization method. We demonstrate our algorithm for the stochastic SIR model of epidemic spread, for which we show that $Q$ can be written as a sum of tensor products. We estimate monthly infection and recovery rates during the first wave of the COVID-19 pandemic in Austria and quantify their uncertainty in a full Bayesian analysis. Availability: Implementation and data are available at https://github.com/spang-lab/TenSIR.

The reconstruction of the structure of biological tissue using electromyographic data is a non-invasive imaging method with diverse medical applications. Mathematically, this process is an inverse problem. Furthermore, electromyographic data are highly sensitive to changes in the electrical conductivity that describes the structure of the tissue. Modeling the inevitable measurement error as a stochastic quantity leads to a Bayesian approach. Solving the discretized Bayes-inverse problem means drawing samples from the posterior distribution of parameters, e.g., the conductivity, given measurement data. Using, e.g., a Metropolis-Hastings algorithm for this purpose involves solving the forward problem for different parameter combinations which requires a high computational effort. Low-rank tensor formats can reduce this effort by providing a data-sparse representation of all occurring linear systems of equations simultaneously and allow for their efficient solution. The application of Bayes' theorem proves the well-posedness of the Bayes-inverse problem. The derivation and proof of a low-rank representation of the forward problem allow for the precomputation of all solutions of this problem under certain assumptions, resulting in an efficient and theory-based sampling algorithm. Numerical experiments support the theoretical results, but also indicate that a high number of samples is needed to obtain reliable estimates for the parameters. The Metropolis-Hastings sampling algorithm, using the precomputed forward solution in a tensor format, draws this high number of samples and therefore enables solving problems which are infeasible using classical methods.

The solution of parameter-dependent linear systems, by classical methods, leads to an arithmetic effort that grows exponentially in the number of parameters. This renders the multigrid method, which has a well understood convergence theory, infeasible. A parameter-dependent representation, e.g., a low-rank tensor format, can avoid this exponential dependence, but in these it is unknown how to calculate the inverse directly within the representation. The combination of these representations with the multigrid method requires a parameter-dependent version of the classical multigrid theory and a parameter-dependent representation of the linear system, the smoother, the prolongation and the restriction. A derived parameter-dependent version of the smoothing property, fulfilled by parameter-dependent versions of the Richardson and Jacobi methods, together with the approximation property prove the convergence of the multigrid method for arbitrary parameter-dependent representations. For a model problem low-rank tensor formats represent the parameter-dependent linear system, prolongation and restriction. The smoother, a damped Jacobi method, is directly approximated in the low-rank tensor format by using exponential sums. Proving the smoothing property for this approximation guarantees the convergence of the parameter-dependent method. Numerical experiments for the parameter-dependent model problem, with bounded parameter value range, indicate a grid size independent convergence rate.

Continuous-time Markov chains describing interacting processes exhibit a state space that grows exponentially in the number of processes. This state-space explosion renders the computation or storage of the time-marginal distribution, which is defined as the solution of a certain linear system, infeasible using classical methods. We consider Markov chains whose transition rates are separable functions, which allows for an efficient low-rank tensor representation of the operator of this linear system. Typically, the right-hand side also has low-rank structure, and thus we can reduce the cost for computation and storage from exponential to linear. Previously known iterative methods also allow for low-rank approximations of the solution but are unable to guarantee that its entries sum up to one as required for a probability distribution. We derive a convergent iterative method using low-rank formats satisfying this condition. We also perform numerical experiments illustrating that the marginal distribution is well approximated with low rank., Comment: 20 pages, 10 figures

Low-rank tensor approximations have shown great potential for uncertainty quantification in high dimensions, for example, to build surrogate models that can be used to speed up large-scale inference problems (Eigel et al., Inverse Problems 34, 2018; Dolgov et al., Statistics & Computing 30, 2020). The feasibility and efficiency of such approaches depends critically on the rank that is necessary to represent or approximate the underlying distribution. In this paper, a-priori rank bounds for approximations in the functional tensor-train representation for the case of Gaussian models are developed. It is shown that under suitable conditions on the precision matrix, the Gaussian density can be approximated to high accuracy without suffering from an exponential growth of complexity as the dimension increases. These results provide a rigorous justification of the suitability and the limitations of low-rank tensor methods in a simple but important model case. Numerical experiments confirm that the rank bounds capture the qualitative behavior of the rank structure when varying the parameters of the precision matrix and the accuracy of the approximation. Finally, the practical relevance of the theoretical results is demonstrated in the context of a Bayesian filtering problem., Comment: 25 pages, 6 figures

We present an iterative method for the search of extreme entries in low-rank tensors which is based on a power iteration combined with a binary search. In this work we use the HT-format for low-rank tensors but other low-rank formats can be used verbatim. We present two different approaches to accelerate the basic power iteration: an orthogonal projection of Rayleigh-Ritz type, as well as an acceleration of the power iteration itself which can be achieved due to the diagonal structure of the underlying eigenvalue problem. Finally the maximizing index is determined by a binary search based on the proposed iterative method for the approximation of the maximum norm. The iterative method for the maximum norm estimation inherits the linear complexity of the involved tensor arithmetic in the HT-format w.r.t. the tensor order, which is also verified by numerical tests.

The evaluation of robustness and reliability of realistic structures in the presence of uncertainty involves costly numerical simulations with a very high number of evaluations. This motivates model order reduction techniques like the proper orthogonal decomposition. When only a few quantities are of interest an approximative mapping from the high-dimensional parameter space onto each quantity of interest is sufficient. Appropriate methods for this task are for instance the polynomial chaos expansion or low-rank tensor approximations. In this work we focus on a non-linear neo-hookean deformation problem with the maximal deformation as our quantity of interest. POD and adaptive POD models of this problem are constructed and compared with respect to approximation quality and construction cost. Additionally, the adapative proper orthogonal decomposition (APOD) is introduced and compared to the regular POD. Building upon that, several hierarchical Tucker approximations (HTAs) are constructed from the reduced and unreduced models. A simple Monte Carlo method in combination with HTA and (A)POD is used to estimate the mean and variance of our quantity of interest. Furthermore, the HTA of the unreduced model is employed to find feasible snapshots for (A)POD.

We consider tensors in the Hierarchical Tucker format and suppose the tensor data to be distributed among several compute nodes. We assume the compute nodes to be in a one-to-one correspondence with the nodes of the Hierarchical Tucker format such that connected nodes can communicate with each other. An appropriate tree structure in the Hierarchical Tucker format then allows for the parallelization of basic arithmetic operations between tensors with a parallel runtime which grows like $\log(d)$, where $d$ is the tensor dimension. We introduce parallel algorithms for several tensor operations, some of which can be applied to solve linear equations $\mathcal{A}X=B$ directly in the Hierarchical Tucker format using iterative methods like conjugate gradients or multigrid. We present weak scaling studies, which provide evidence that the runtime of our algorithms indeed grows like $\log(d)$. Furthermore, we present numerical experiments in which we apply our algorithms to solve a parameter-dependent diffusion equation in the Hierarchical Tucker format by means of a multigrid algorithm.

Low rank tensor completion is a highly ill-posed inverse problem, particularly when the data model is not accurate, and some sort of regularization is required in order to solve it. In this article we focus on the calibration of the data model. For alternating optimization, we observe that existing rank adaption methods do not enable a continuous transition between manifolds of different ranks. We denote this characteristic as $\textit{instability (under truncation)}$. As a consequence of this property, arbitrarily small changes in the iterate can have arbitrarily large influence on the further reconstruction. We therefore introduce a singular value based regularization to the standard alternating least squares (ALS), which is motivated by averaging in microsteps. We prove its $\textit{stability}$ and derive a natural semi-implicit rank adaption strategy. We further prove that the standard ALS microsteps for completion problems are only stable on manifolds of fixed ranks, and only around points that have what we define as $\textit{internal tensor restricted isometry property, iTRIP}$. In conclusion, numerical experiments are provided that show improvements of the reconstruction quality up to orders of magnitude in the new Stable ALS Approximation (SALSA) compared to standard ALS and the well known Riemannian optimization RTTC., Comment: 43 pages

We consider the problem of fitting a low rank tensor $A\in\mathbb{R}^{{\mathcal I}}$, ${\mathcal I} = \{1,\ldots,n\}^{d}$, to a given set of data points $\{M_i\in\mathbb{R}\mid i\in P\}$, $P\subset{\mathcal I}$. The low rank format under consideration is the hierarchical or TT or MPS format. It is characterized by rank bounds $r$ on certain matricizations of the tensor. The number of degrees of freedom is in ${\cal O}(r^2dn)$. For a fixed rank and mode size $n$ we observe that it is possible to reconstruct random (but rank structured) tensors as well as certain discretized multivariate (but rank structured) functions from a number of samples that is in ${\cal O}(\log N)$ for a tensor having $N=n^d$ entries. We compare an alternating least squares fit (ALS) to an overrelaxation scheme inspired by the LMaFit method for matrix completion. Both approaches aim at finding a tensor $A$ that fulfils the first order optimality conditions by a nonlinear Gauss-Seidel type solver that consists of an alternating fit cycling through the directions $\mu=1,\ldots,d$. The least squares fit is of complexity ${\cal O}(r^4d\#P)$ per step, whereas each step of ADF is in ${\cal O}(r^2d\#P)$, albeit with a slightly higher number of necessary steps. In the numerical experiments we observe robustness of the completion algorithm with respect to noise and good reconstruction capability. Our tests provide evidence that the algorithm is suitable in higher dimension ($>$10) as well as for moderate ranks. Keywords: MPS, Tensor Completion, Tensor Train, TT, Hierarchical Tucker, HT, ALS., Comment: (Accepted for publication in SIAM Journal on Scientific Computing (SISC))

In this paper, we develop a multigrid method on unstructured shape-regular grids. For a general shape-regular unstructured grid of ${\cal O}(N)$ elements, we present a construction of an auxiliary coarse grid hierarchy on which a geometric multigrid method can be applied together with a smoothing on the original grid by using the auxiliary space preconditioning technique. Such a construction is realized by a cluster tree which can be obtained in ${\cal O}(N\log N)$ operations for a grid of $N$ elements. This tree structure in turn is used for the definition of the grid hierarchy from coarse to fine. For the constructed grid hierarchy we prove that the convergence rate of the multigrid preconditioned CG for an elliptic PDE is $1 - {\cal O}({1}/{\log N})$. Numerical experiments confirm the theoretical bounds and show that the total complexity is in ${\cal O}(N\log N)$.

A recurring theme in attempts to break the curse of dimensionality in the numerical approximations of solutions to high-dimensional partial differential equations (PDEs) is to employ some form of sparse tensor approximation. Unfortunately, there are only a few results that quantify the possible advantages of such an approach. This paper introduces a class $\Sigma_n$ of functions, which can be written as a sum of rank-one tensors using a total of at most $n$ parameters and then uses this notion of sparsity to prove a regularity theorem for certain high-dimensional elliptic PDEs. It is shown, among other results, that whenever the right-hand side $f$ of the elliptic PDE can be approximated with a certain rate $\mathcal{O}(n^{-r})$ in the norm of ${\mathrm H}^{-1}$ by elements of $\Sigma_n$, then the solution $u$ can be approximated in ${\mathrm H}^1$ from $\Sigma_n$ to accuracy $\mathcal{O}(n^{-r'})$ for any $r'\in (0,r)$. Since these results require knowledge of the eigenbasis of the elliptic operator considered, we propose a second "basis-free" model of tensor sparsity and prove a regularity theorem for this second sparsity model as well. We then proceed to address the important question of the extent such regularity theorems translate into results on computational complexity. It is shown how this second model can be used to derive computational algorithms with performance that breaks the curse of dimensionality on certain model high-dimensional elliptic PDEs with tensor-sparse data., Comment: 41 pages, 1 figure

During the last years, low-rank tensor approximation has been established as a new tool in scientific computing to address large-scale linear and multilinear algebra problems, which would be intractable by classical techniques. This survey attempts to give a literature overview of current developments in this area, with an emphasis on function-related tensors.