No Arabic abstract
In this paper, we develop efficient and accurate algorithms for evaluating $varphi(A)$ and $varphi(A)b$, where $A$ is an $Ntimes N$ matrix, $b$ is an $N$ dimensional vector and $varphi$ is the function defined by $varphi(x)equivsumlimits^{infty}_{k=0}frac{z^k}{(1+k)!}$. Such matrix function (the so-called $varphi$-function) plays a key role in a class of numerical methods well-known as exponential integrators. The algorithms use the scaling and modified squaring procedure combined with truncated Taylor series. The backward error analysis is presented to find the optimal value of the scaling and the degree of the Taylor approximation. Some useful techniques are employed for reducing the computational cost. Numerical comparisons with state-of-the-art algorithms show that the algorithms perform well in both accuracy and efficiency.
Evaluating the log-sum-exp function or the softmax function is a key step in many modern data science algorithms, notably in inference and classification. Because of the exponentials that these functions contain, the evaluation is prone to overflow and underflow, especially in low precision arithmetic. Software implementations commonly use alternative formulas that avoid overflow and reduce the chance of harmful underflow, employing a shift or another rewriting. Although mathematically equivalent, these variants behave differently in floating-point arithmetic. We give rounding error analyses of different evaluation algorithms and interpret the error bounds using condition numbers for the functions. We conclude, based on the analysis and numerical experiments, that the shifted formulas are of similar accuracy to the unshifted ones and that the shifted softmax formula is typically more accurate than a division-free variant.
Recent studies have demonstrated improved skill in numerical weather prediction via the use of spatially correlated observation error covariance information in data assimilation systems. In this case, the observation weighting matrices (inverse error covariance matrices) used in the assimilation may be full matrices rather than diagonal. Thus, the computation of matrix-vector products in the variational minimization problem may be very time-consuming, particularly if the parallel computation of the matrix-vector product requires a high degree of communication between processing elements. Hence, we introduce a well-known numerical approximation method, called the fast multipole method (FMM), to speed up the matrix-vector multiplications in data assimilation. We explore a particular type of FMM that uses a singular value decomposition (SVD-FMM) and adjust it to suit our new application in data assimilation. By approximating a large part of the computation of the matrix-vector product, the SVD-FMM technique greatly reduces the computational complexity compared with the standard approach. We develop a novel possible parallelization scheme of the SVD-FMM for our application, which can reduce the communication costs. We investigate the accuracy of the SVD-FMM technique in several numerical experiments: we first assess the accuracy using covariance matrices that are created using different correlation functions and lengthscales; then investigate the impact of reconditioning the covariance matrices on the accuracy; and finally examine the feasibility of the technique in the presence of missing observations. We also provide theoretical explanations for some numerical results. Our results show that the SVD-FMM technique has potential as an efficient technique for assimilation of a large volume of observational data within a short time interval.
In simulations of fluid motion time accuracy has proven to be elusive. We seek highly accurate methods with strong enough stability properties to deal with the richness of scales of many flows. These methods must also be easy to implement within current complex, possibly legacy codes. Herein we develop, analyze and test new time stepping methods addressing these two issues with the goal of accelerating the development of time accurate methods addressing the needs of applications. The new methods are created by introducing inexpensive pre-filtering and post-filtering steps to popular methods which have been implemented and tested within existing codes. We show that pre-filtering and post-filtering a multistep or multi-stage method results in new methods which have both multiple steps and stages: these are general linear methods (GLMs). We utilize the well studied properties of GLMs to understand the accuracy and stability of filtered method, and to design optimal new filters for popular time-stepping methods. We present several new embedded families of high accuracy methods with low cognitive complexity and excellent stability properties. Numerical tests of the methods are presented, including ones finding failure points of some methods. Among the new methods presented is a novel pair of alternating filters for the Implicit Euler method which induces a third order, A-stable, error inhibiting scheme which is shown to be particularly effective.
We propose an efficient, accurate and robust implicit solver for the incompressible Navier-Stokes equations, based on a DG spatial discretization and on the TR-BDF2 method for time discretization. The effectiveness of the method is demonstrated in a number of classical benchmarks, which highlight its superior efficiency with respect to other widely used implicit approaches. The parallel implementation of the proposed method in the framework of the deal.II software package allows for accurate and efficient adaptive simulations in complex geometries, which makes the proposed solver attractive for large scale industrial applications.
Fourier extension is an approximation method that alleviates the periodicity requirements of Fourier series and avoids the Gibbs phenomenon when approximating functions. We describe a similar extension approach using regular wavelet bases on a hypercube to approximate functions on subsets of that cube. These subsets may have a general shape. This construction is inherently associated with redundancy which leads to severe ill-conditioning, but recent theory shows that nevertheless high accuracy and numerical stability can be achieved using regularization and oversampling. Regularized least squares solvers, such as the truncated singular value decomposition, that are suited to solve the resulting ill-conditioned and skinny linear system generally have cubic computational cost. We compare several algorithms that improve on this complexity. The improvements benefit from the sparsity in and the structure of the discrete wavelet transform. We present a method that requires $mathcal O(N)$ operations in 1-D and $mathcal O(N^{3(d-1)/d})$ in $d$-D, $d>1$. We experimentally show that direct sparse QR solvers appear to be more time-efficient, but yield larger expansion coefficients.