Tensor Product Domain Research Articles

Approximation in the extended functional tensor train format

This work proposes the extended functional tensor train (EFTT) format for compressing and working with multivariate functions on tensor product domains. Our compression algorithm combines tensorized Chebyshev interpolation with a low-rank approximation algorithm that is entirely based on function evaluations. Compared to existing methods based on the functional tensor train format, the adaptivity of our approach often results in reducing the required storage, sometimes considerably, while achieving the same accuracy. In particular, we reduce the number of function evaluations required to achieve a prescribed accuracy by up to over 96%\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$96\\%$$\\end{document} compared to the algorithm from Gorodetsky et al. (Comput. Methods Appl. Mech. Eng. 347, 59–84 2019).

A μ-mode BLAS approach for multidimensional tensor-structured problems

In this manuscript, we present a common tensor framework which can be used to generalize one-dimensional numerical tasks to arbitrary dimension d by means of tensor product formulas. This is useful, for example, in the context of multivariate interpolation, multidimensional function approximation using pseudospectral expansions and solution of stiff differential equations on tensor product domains. The key point to obtain an efficient-to-implement BLAS formulation consists in the suitable usage of the μ-mode product (also known as tensor-matrix product or mode-n product) and related operations, such as the Tucker operator. Their MathWorks MATLAB®/GNU Octave implementations are discussed in the paper, and collected in the package KronPACK. We present numerical results on experiments up to dimension six from different fields of numerical analysis, which show the effectiveness of the approach.

Removal of spurious outlier frequencies and modes from isogeometric discretizations of second- and fourth-order problems in one, two, and three dimensions

A key advantage of isogeometric discretizations is their accurate and well-behaved eigenfrequencies and eigenmodes. For degree two and higher, however, a few spurious modes appear that possess inaccurate frequencies, denoted as “outliers”. The outlier frequencies and corresponding modes are at the root of several efficiency and robustness issues in isogeometric analysis. One example is explicit dynamics where outlier frequencies unnecessarily reduce the critical time step. Another example is wave propagation where the inaccurate outlier modes may participate in the solution. In this paper, we first investigate the spurious outlier frequencies and corresponding modes of isogeometric discretizations of second- and fourth-order model problems and provide a complete characterization. We then devise a new approach that removes all outliers modes without negatively affecting the accuracy of the discretizations. Our approach is variationally consistent and works for a range of common boundary conditions on tensor product domains. We finally demonstrate that our approach allows a much larger critical time step, irrespective of polynomial degree, providing a pathway towards efficient higher-order explicit dynamics.

Functional Tucker Approximation Using Chebyshev Interpolation

This work is concerned with approximating a trivariate function defined on a tensor-product domain via function evaluations. Combining tensorized Chebyshev interpolation with a Tucker decomposition of low multilinear rank yields function approximations that can be computed and stored very efficiently. The existing Chebfun3 algorithm [Hashemi and Trefethen, SIAM J. Sci. Comput., 39 (2017)]uses a similar format but the construction of the approximation proceeds indirectly, via a so called slice-Tucker decomposition. As a consequence, Chebfun3 sometimes uses unnecessarily many function evaluations and does not fully benefit from the potential of the Tucker decomposition to reduce, sometimes dramatically, the computational cost. We propose a novel algorithm Chebfun3F that utilizes univariate fibers instead of bivariate slices to construct the Tucker decomposition. Chebfun3F reduces the cost for the approximation in terms of the number of function evaluations for nearly all functions considered, typically by 75%, and sometimes by over 98%.

Two-scale Finite Element Discretizations for Semilinear Parabolic Equations

In this paper, to reduce the computational cost of solving semilinear parabolic equations on a tensor product domain Ω⊂ℝ^<i>d</i> with <i>d</i> = 2 or 3, some two-scale finite element discretizations are proposed and analyzed. The time derivative in semilinear parabolic equations is approximated by the backward Euler finite difference scheme. The two-scale finite element method is designed for the space discretization. The idea of the two-scale finite element method is based on an understanding of a finite element solution to an elliptic problem on a tensor product domain. The high frequency parts of the finite element solution can be well captured on some univariate fine grids and the low frequency parts can be approximated on a coarse grid. Thus the two-scale finite element approximation is defined as a linear combination of some standard finite element approximations on some univariate fine grids and a coarse grid satisfying <i>H</i> = <i>O</i> (<i>h</i>^1/2), where <i>h</i> and <i>H</i> are the fine and coarse mesh widths, respectively. It is shown theoretically and numerically that the backward Euler two-scale finite element solution not only achieves the same order of accuracy in the <i>H</i>¹ (Ω) norm as the backward Euler standard finite element solution, but also reduces the number of degrees of freedom from <i>O</i>(<i>h</i>^-<i>d</i>×<i>τ</i>^-1) to <i>O</i>(<i>h</i>^{-(<i>(d)</i>+1)/2}×<i>τ</i>^-1) where <i>τ</i> is the time step. Consequently the backward Euler two-scale finite element method for semilinear parabolic equations is more efficient than the backward Euler standard finite element method.

APPROXIMATING SMOOTH, MULTIVARIATE FUNCTIONS ON IRREGULAR DOMAINS

In this paper, we introduce a method known aspolynomial frame approximationfor approximating smooth, multivariate functions defined on irregular domains in$d$dimensions, where$d$can be arbitrary. This method is simple, and relies only on orthogonal polynomials on a bounding tensor-product domain. In particular, the domain of the function need not be known in advance. When restricted to a subdomain, an orthonormal basis is no longer a basis, but a frame. Numerical computations with frames present potential difficulties, due to the near-linear dependence of the truncated approximation system. Nevertheless, well-conditioned approximations can be obtained via regularization, for instance, truncated singular value decompositions. We comprehensively analyze such approximations in this paper, providing error estimates for functions with both classical and mixed Sobolev regularity, with the latter being particularly suitable for higher-dimensional problems. We also analyze the sample complexity of the approximation for sample points chosen randomly according to a probability measure, providing estimates in terms of the correspondingNikolskii inequalityfor the domain. In particular, we show that the sample complexity for points drawn from the uniform measure is quadratic (up to a log factor) in the dimension of the polynomial space, independently of $d$, for a large class of nontrivial domains. This extends a well-known result for polynomial approximation in hypercubes.

Compressive Hermite Interpolation: Sparse, High-Dimensional Approximation from Gradient-Augmented Measurements

We consider the sparse polynomial approximation of a multivariate function on a tensor product domain from samples of both the function and its gradient. When only function samples are prescribed, weighted $\ell^1$ minimization has recently been shown to be an effective procedure for computing such approximations. We extend this work to the gradient-augmented case. Our main results show that for the same asymptotic sample complexity, gradient-augmented measurements achieve an approximation error bound in a stronger Sobolev norm, as opposed to the $L^2$-norm in the unaugmented case. For Chebyshev and Legendre polynomial approximations, this sample complexity estimate is algebraic in the sparsity $s$ and at most logarithmic in the dimension $d$, thus mitigating the curse of dimensionality to a substantial extent. We also present several experiments numerically illustrating the benefits of gradient information over an equivalent number of function samples only.

Hierarchical matrix approximation for the uncertainty quantification of potentials on random domains

Computing statistical quantities of interest of the solution of partial differential equations on random domains is an important and challenging task in engineering. We consider the computation of these quantities by the perturbation approach. Especially, we discuss how third order accurate expansions of the mean and the correlation can numerically be computed. These expansions become even fourth order accurate for certain types of boundary variations. The correction terms are given by the solution of correlation equations in the tensor product domain, which can efficiently be computed by means of H-matrices. They have recently been shown to be an efficient tool to solve correlation equations with rough data correlations, that is, with low Sobolev smoothness or small correlation length, in almost linear time. Numerical experiments in three dimensions for higher order ansatz spaces show the feasibility of the proposed algorithm. The application to a non-smooth domain is also included.

Novel results for the anisotropic sparse grid quadrature

This article is dedicated to the anisotropic sparse grid quadrature for functions which are analytically extendable into an anisotropic tensor product domain. Taking into account this anisotropy, we end up with a dimension independent error versus cost estimate of the proposed quadrature. In addition, we provide a novel and improved estimate for the cardinality of the underlying anisotropic index set. To validate the theoretical findings, we present several examples ranging from simple quadrature problems to diffusion problems on random domains. These examples demonstrate the remarkable convergence behavior of the anisotropic sparse grid quadrature in applications.

Function Approximation on Arbitrary Domains Using Fourier Extension Frames

Fourier extension is an approximation scheme in which a function on an arbitrary bounded domain is approximated using classical Fourier series on a bounding box. On the domain of interest the Fourier series exhibits redundancy, and it has the mathematical structure of a frame rather than a basis. It is not trivial to construct approximations in this frame using function evaluations in points that belong to the domain only, but one way to do so is through a discrete least squares approximation. The corresponding system is extremely ill-conditioned, due to the redundancy in the frame, yet its solution via a regularized singular value decomposition is known to be accurate to very high (and nearly spectral) precision. Still, this computation requires $\mathcal{O}(N^3)$ operations. In this paper we describe an algorithm to compute such Fourier extension frame approximations in only $\mathcal{O}(N^2\log^2N)$ operations for general 2D domains. The cost improves to $\mathcal{O}(N\log^2N)$ operations for simpler tensor-product domains. The algorithm exploits a phenomenon called the plunge region in the analysis of time-frequency localization operators, which manifests itself here as a sudden drop in the singular values of the least squares matrix. It is known that the size of the plunge region scales like $\mathcal{O}(\log N)$ in 1D problems. In this paper we show that for most 2D domains in the fully discrete case the plunge region scales like $\mathcal{O}(N\log N)$, proving a discrete equivalent of a result that was conjectured by Widom for a related continuous problem. The complexity estimate depends on the Minkowski or box-counting dimension of the domain boundary, and as such it is larger than $\mathcal{O}(N\log N)$ for domains with fractal shape.

Statistical Solutions of Hyperbolic Conservation Laws: Foundations

We seek to define statistical solutions of hyperbolic systems of conservation laws as time-parametrized probability measures on $p$-integrable functions. To do so, we prove the equivalence between probability measures on $L^p$ spaces and infinite families of \textit{correlation measures}. Each member of this family, termed a \textit{correlation marginal}, is a Young measure on a finite-dimensional tensor product domain and provides information about multi-point correlations of the underlying integrable functions. We also prove that any probability measure on a $L^p$ space is uniquely determined by certain moments (correlation functions) of the equivalent correlation measure. We utilize this equivalence to define statistical solutions of multi-dimensional conservation laws in terms of an infinite set of equations, each evolving a moment of the correlation marginal. These evolution equations can be interpreted as augmenting entropy measure-valued solutions, with additional information about the evolution of all possible multi-point correlation functions. Our concept of statistical solutions can accommodate uncertain initial data as well as possibly non-atomic solutions even for atomic initial data. For multi-dimensional scalar conservation laws we impose additional entropy conditions and prove that the resulting \textit{entropy statistical solutions} exist, are unique and are stable with respect to the $1$-Wasserstein metric on probability measures on $L^1$.

Spectral semi-implicit and space–time discontinuous Galerkin methods for the incompressible Navier–Stokes equations on staggered Cartesian grids

In this paper two new families of arbitrary high order accurate spectral discontinuous Galerkin (DG) finite element methods are derived on staggered Cartesian grids for the solution of the incompressible Navier–Stokes (NS) equations in two and three space dimensions. The discrete solutions of pressure and velocity are expressed in the form of piecewise polynomials along different meshes. While the pressure is defined on the control volumes of the main grid, the velocity components are defined on edge-based dual control volumes, leading to a spatially staggered mesh. Thanks to the use of a nodal basis on a tensor-product domain, all discrete operators can be written efficiently as a combination of simple one-dimensional operators in a dimension-by-dimension fashion.In the first family, high order of accuracy is achieved only in space, while a simple semi-implicit time discretization is derived by introducing an implicitness factor θ∈[0.5,1] for the pressure gradient in the momentum equation. The real advantages of the staggering arise after substituting the discrete momentum equation into the weak form of the continuity equation. In fact, the resulting linear system for the pressure is symmetric and positive definite and either block penta-diagonal (in 2D) or block hepta-diagonal (in 3D). As a consequence, the pressure system can be solved very efficiently by means of a classical matrix-free conjugate gradient method. From our numerical experiments we find that the pressure system appears to be reasonably well-conditioned, since in all test cases shown in this paper the use of a preconditioner was not necessary. This is a rather unique feature among existing implicit DG schemes for the Navier–Stokes equations. In order to avoid a stability restriction due to the viscous terms, the latter are discretized implicitly using again a staggered mesh approach, where the viscous stress tensor is also defined on the dual mesh.The second family of staggered DG schemes proposed in this paper achieves high order of accuracy also in time by expressing the numerical solution in terms of piecewise space–time polynomials. In order to circumvent the low order of accuracy of the adopted fractional stepping, a simple iterative Picard procedure is introduced, which leads to a space–time pressure-correction algorithm. In this manner, the symmetry and positive definiteness of the pressure system are not compromised. The resulting algorithm is stable, computationally very efficient, and at the same time arbitrary high order accurate in both space and time. These features are typically not easy to obtain all at the same time for a numerical method applied to the incompressible Navier–Stokes equations. The new numerical method has been thoroughly validated for approximation polynomials of degree up to N=11, using a large set of non-trivial test problems in two and three space dimensions, for which either analytical, numerical or experimental reference solutions exist.

Novel and efficient computation of Hilbert–Huang transform on surfaces

Hilbert–Huang Transform (HHT) has proven to be extremely powerful for signal processing and analysis in 1D time series, and its generalization to regular tensor-product domains (e.g., 2D and 3D Euclidean space) has also demonstrated its widespread utilities in image processing and analysis. Compared with popular Fourier transform and wavelet transform, the most prominent advantage of Hilbert–Huang Transform (HHT) is that, it is a fully data-driven, adaptive method, especially valuable for handling non-stationary and nonlinear signals. Two key technical elements of Hilbert–Huang transform are: (1) Empirical Mode Decomposition (EMD) and (2) Hilbert spectra computation. HHT's uniqueness results from its capability to reveal both global information (i.e., Intrinsic Mode Functions (IMFs) enabled by EMD) and local information (i.e., the computation of local frequency, amplitude (energy) and phase information enabled by Hilbert spectra computation) from input signals. Despite HHT's rapid advancement in the past decade, its theory and applications on surfaces remain severely under-explored due to the current technical challenge in conducting Hilbert spectra computation on surfaces. To ameliorate, this paper takes a new initiative to compute the Riesz transform on 3D surfaces, a natural generalization of Hilbert transform in higher-dimensional cases, with a goal to make the theoretic breakthrough. The core of our theoretic and computational framework is to fully exploit the relationship between Riesz transform and fractional Laplacian operator, which can enable the computation of Riesz transform on surfaces via eigenvalue decomposition of Laplacian matrix. Moreover, we integrate the techniques of EMD and our newly-proposed Riesz transform on 3D surfaces by monogenic signals to compute Hilbert spectra, which include the space-frequency-energy distribution of signals defined over 3D surfaces and characterize key local feature information (e.g., instantaneous frequency, local amplitude, and local phase). Experiments and applications in spectral geometry processing and prominent feature detection illustrate the effectiveness of the current computational framework of HHT on 3D surfaces, which could serve as a solid foundation for upcoming, more serious applications in graphics and geometry computing fields.

Multidimensional Summation-by-Parts Operators: General Theory and Application to Simplex Elements

Summation-by-parts (SBP) finite-difference discretizations share many attractive properties with Galerkin finite-element methods (FEMs), including time stability and superconvergent functionals; however, unlike FEMs, SBP operators are not completely determined by a basis, so the potential exists to tailor SBP operators to meet different objectives. To date, application of high-order SBP discretizations to multiple dimensions has been limited to tensor product domains. This paper presents a definition for multidimensional SBP finite-difference operators that is a natural extension of one-dimensional SBP operators. Theoretical implications of the definition are investigated for the special case of a diagonal-norm (mass) matrix. In particular, a diagonal-norm SBP operator exists on a given domain if and only if there is a cubature rule with positive weights on that domain and the polynomial-basis matrix has full rank when evaluated at the cubature nodes. Appropriate simultaneous-approximation terms are devel...

Two-scale sparse finite element approximations

To reduce computational cost, we study some two-scale finite element approximations on sparse grids for elliptic partial differential equations of second order in a general setting. Over any tensor product domain $\Omega\subset\mathbb{R}^d$ with $d= 2,3$, we construct the two-scale finite element approximations for both boundary value and eigenvalue problems by using a Boolean sum of some existing finite element approximations on a coarse grid and some univariate fine grids and hence they are cheaper approximations. As applications, we obtain some new efficient finite element discretizations for the two classes of problem: The new two-scale finite element approximation on a sparse grid not only has the less degrees of freedom but also achieves a good accuracy of approximation.

Direct discretizations of bi-variate population balance systems with finite difference schemes of different order

The accurate and efficient simulation of bi-variate population balance systems is nowadays a great challenge since the domain spanned by the external and internal coordinates is five-dimensional. This report considers direct discretizations of this equation in tensor-product domains. In this situation, finite difference methods can be applied. The studied model includes the transport of dissolved potassium dihydrogen phosphate (KDP) and of energy (temperature) in a laminar flow field as well as the nucleation and growth of KDP particles. Two discretizations of the coupled model will be considered which differ only in the discretization of the population balance equation: a first order monotone upwind scheme and a third order essentially non-oscillatory (ENO) scheme. The Dirac term on the right-hand side of this equation is discretized with a finite volume method. The numerical results show that much different results are obtained even in the class of direct discretizations.

On the construction of sparse tensor product spaces

Let Ω1 ⊂ Rn1 and Ω2 ⊂ Rn2 be two given domains and consider on each domain a multiscale sequence of ansatz spaces of polynomial exactness r1 and r2, respectively. In this paper, we study the optimal construction of sparse tensor products made from these spaces. In particular, we derive the resulting cost complexities to approximate functions with anisotropic and isotropic smoothness on the tensor product domain Ω1×Ω2. Numerical results validate our theoretical findings.

Multiple point evaluation on combined tensor product supports

We consider the multiple point evaluation problem for an n-dimensional space of functions [???1,1[ d ?? spanned by d-variate basis functions that are the restrictions of simple (say linear) functions to tensor product domains. For arbitrary evaluation points this task is faced in the context of (semi-)Lagrangian schemes using adaptive sparse tensor approximation spaces for boundary value problems in moderately high dimensions. We devise a fast algorithm for performing m???n point evaluations of a function in this space with computational cost O(mlog d n). We resort to nested segment tree data structures built in a preprocessing stage with an asymptotic effort of O(nlog d???1 n).

Postprocessed Two-Scale Finite Element Discretizations, Part I

This paper studies a simple postprocessing of two-scale finite element discretizations of elliptic partial differential operators, where boundary value problems and eigenvalue problems on tensor-product domains $\Omega\subset\mathbb{R}^d$ are both examined and conforming piecewise d-linear elements are used. For $d\geq2$, suppose that the two-scale method is applied using a coarse mesh and some univariate fine meshes satisfying $H=\mathcal{O}(h^{2/3})$, where h and H are the fine and coarse mesh widths. It is shown for both classes of problem that the postprocessed two-scale finite element solution achieves the same order of accuracy in the $H^1(\Omega)$ norm as the postprocessed standard finite element solution, but the former solution uses only $\mathcal{O}(h^{-(2d+1)/3})$ degrees of freedom compared with the $\mathcal{O}(h^{-d})$ degrees of freedom required by the latter.

Three-scale finite element eigenvalue discretizations

Some three-scale finite element discretization schemes are proposed and analyzed in this paper for a class of elliptic eigenvalue problems on tensor product domains. With these schemes, the solution of an eigenvalue problem on a fine grid may be reduced to the solutions of eigenvalue problems on a relatively coarse grid and some partially mesoscopic grids, together with the solutions of linear algebraic systems on a globally mesoscopic grid and several partially fine grids. It is shown theoretically and numerically that this type of discretization schemes not only significantly reduce the number of degrees of freedom but also produce very accurate approximations.