Weighted high dimensional data reduction of finite element features: an application on high pressure of an abdominal aortic aneurysm

In this work we propose a low rank approximation of areal, particularly three dimensional, data utilizing additional weights. This way we enable effective compression if additional information indicates that parts of the data are of higher interest than others. The guiding example are high fidelity finite element simulations of an abdominal aortic aneurysm, i.e. a deformed blood vessel. The additional weights encapsulate the areas of high stress, which we assume indicates the rupture risk of the aorta. The stress values on the grid are modeled as a Gaussian Markov random field and we define our approximation as a basis of vectors that solve a series of optimization problems. Each of these problems describes the minimization of an expected weighted quadratic loss. We provide an effective numerical heuristic to compute the basis under general conditions, which relies on the sparsity of the precision matrix to ensure acceptable computing time even for large grids. We explicitly explore two such bases on the surface of a high fidelity finite element grid and show their efficiency for compression. Finally, we utilize the approach as part of a larger model to predict the van Mises stress in areas of interest using low and high fidelity simulations.

No funds, grants, or other support was received.

Authors and Affiliations

1. Department of Statistics, Faculty of Mathematics, Informatics and Statistics, Ludwig-Maximilians-Universität München, Ludwigstr. 33, 80539, Munich, Germany

   Christoph Striegel & Göran Kauermann

2. Institute for Computational Mechanics, Department of Mechanical Engineering, Technical University of Munich, Munich, Germany

   Jonas Biehler

Corresponding author

Correspondence to Christoph Striegel.

Code and Data

full code and data to reproduce the results is available at: https://drive.google.com/drive/folders/1c68u1bATuOZEIAXuOYzIFb3urSNQwUav?usp=sharing.

Reformulation as Tensor Decomposition

reformulation_tensor_decomposition.pdf contains a reformulation of the first basis vector problem as a tensor decomposition.

Simulation Study

simulation_study.pdf contains as small simulation study with sampled weight vectors.

Histogram of the number of neighbors for the high fidelity grid points

neighbors_count.pdf.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix

1.1 Miscellaneous numerics

There are a number of crucial numerical concerns that need to be addressed in order to make the basis computation problem laid out in the previous parts solvable in practice. Firstly, we only know the precision matrix, \(\varvec{Q}\), but not its inverse the GMRF covariance matrix \(\varvec{\varSigma }\), as seen in (8). \(\varvec{Q}\), the standard GMRF precision matrix (Rue and Held 2005) is sparse. Moreover, this matrix is not invertible in a strict sense, as the constant vector is an eigenvector with eigenvalue 0. We can however fix this by adding an identity matrix, i.e. in this work we use

where \(\epsilon \) is a small constant, i.e. \(10^{-4}\) and \(\varvec{I}_{m_{y}}\) the \(m_{y}\) dimensional identity matrix.

It is not useful nor common to actually invert \({\tilde{\varvec{Q}}}\) as this is computationally expensive and results in a dense matrix which consumes a large amount of memory. Instead we compute a sparse Cholesky decomposition which is enough to rapidly compute the matrix / vector operations involving \(\varvec{\varSigma }\) in the formulae for the function value and gradient, like (8) and (9). The decomposition is not sufficient to actually compute the full Hessian matrix as seen in (10). However, this is not necessary as the ability to compute the product of the Hessian with a given vector is sufficient in our case.

There are specialized routines for the case of GMRF type precision matrices which do not only build on the highly sparse nature but also the (up to permutation) band like shape of the matrix in order to achieve a computational cost of up to \(O(n^{3/2})\). This results in a large advantage with regard to computation time in comparison to the ordinary Cholesky decomposition routines which also do not result in sparse matrices. We refer to Davis (2006) as a comprehensive reference. For the computation we use the R package Matrix, which is based on the CHOLMOD (Davis 2021) library. On our machine with 64 GB of RAM and 12 core AMD Ryzen \(9 \quad 3900\)X processor the computation of the decomposition for our 80,000 dimensional precision matrix takes just 40 s.

1.2 Low fidelity compression error

See Fig. 11.

Squared compression error for the low fidelity data using a simple (unweighted) basis of GMRF eigenvectors. Error defined analogously to (18)

Full size image

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions