MIMS EPrints
Not a member yet
2151 research outputs found
Sort by
Algorithms for the rational approximation of matrix-valued functions
Full text linkA selection of algorithms for the rational approximation of matrix-valued functions are discussed, including variants of the interpolatory AAA method, the RKFIT method based on approximate least squares fitting, vector fitting, and a method based on low-rank approximation of a block Loewner matrix. A new method, called the block-AAA algorithm, based on a generalized barycentric formula with matrix-valued weights is proposed. All algorithms are compared in terms of obtained approximation accuracy and runtime on a set of problems from model order reduction and nonlinear eigenvalue problems, including examples with noisy data. It is found that interpolation-based methods are typically cheaper to run, but they may suffer in the presence of noise for which approximation-based methods perform better
Quadratic Realizability of Palindromic Matrix Polynomials: the Real Case
No full textLet \cL = (\cL_1,\cL_2) be a list consisting of structural data for a matrix polynomial; here \cL_1 is a sublist consisting of powers of irreducible (monic) scalar polynomials over the field \RR, and \cL_2 is a sublist of nonnegative integers.
For an arbitrary such \cL, we give easy-to-check necessary and sufficient conditions for \cL to be the list of elementary divisors and minimal indices of some real -palindromic quadratic matrix polynomial.
For a list \cL satisfying these conditions, we show how to explicitly build a real -palindromic quadratic matrix polynomial having \cL as its structural data;
that is, we provide a -palindromic quadratic realization of \cL over \RR.
A significant feature of our construction differentiates it from related work in the literature; the realizations constructed here are direct sums of blocks with low bandwidth, that transparently display the spectral and singular structural data in the original list \cL
Stochastic Rounding: Implementation, Error Analysis, and Applications
No full textStochastic rounding randomly maps a real number to one of the two nearest values in a finite precision number system. First proposed for use in computer arithmetic in the 1950s, it is attracting renewed interest. If used in floating-point arithmetic in the computation of the inner product of two vectors of length n, it yields an error bounded by \sqrt(n)u with high probability, where u is the unit roundoff, which is not necessarily the case for round to nearest. A particular attraction of stochastic rounding is that, unlike round to nearest, it is immune to the phenomenon of stagnation, whereby a sequence of tiny updates to a relatively large quantity are lost. We survey stochastic rounding, covering its mathematical properties and probabilistic error analysis, its implementation, and its use in applications, including deep learning and the numerical solution of differential equations
Model order reduction of layered waveguides via rational Krylov fitting
Full text linkNetwork-based data-driven reduced order models have recently emerged as an efficient numerical tool for forward and inverse problems of wave propagation. Currently, this technique is limited to two classes of problems: bounded inhomogeneous domains (with applications in multiscale simulation and imaging) and homogeneous halfspaces (for the solution of exterior forward problems). Here we relax the constant coefficient requirement for the latter by considering reduced order models (ROMs) of unbounded waveguides with layered inclusions, thereby giving rise to efficient discrete perfectly matched layers (PMLs) for nonhomogeneous media. Our approach is based on the solution of a nonlinear rational least squares problem using the RKFIT method [M. Berljafa and S. Güttel, SIAM J. Sci. Comput., 39(5):A2049--A2071, 2017]. We show how the solution of this least squares problem can be converted into an accurate sparse network approximation within a rational Krylov framework. Several numerical experiments are included. They indicate that RKFIT computes PMLs more accurately than previous analytic approaches and even works in regimes where the transfer functions to be approximated are highly irregular due to pronounced scattering resonances. Spectral adaptation effects allow for accurate ROMs with dimensions near or even below the Nyquist limit
Numerical Stability of Algorithms at Extreme Scale and Low Precisions
Full text linkThe largest dense linear systems that are being solved today are of order . Single precision arithmetic, which has a unit roundoff , is widely used in scientific computing, and half precision arithmetic, with , is increasingly being exploited as it becomes more readily available in hardware. Standard rounding error bounds for numerical linear algebra algorithms are proportional to , with growing at least linearly with . Therefore we are at the stage where these rounding error bounds are not able to guarantee any accuracy or stability in the computed results for some extreme-scale or low-accuracy computations. We explain how rounding error bounds with much smaller constants can be obtained. Blocked algorithms, which break the data into blocks of size , lead to a reduction in the error constants by a factor or more. Two architectural features also reduce the error constants: extended precision registers and fused multiply--add operations, either at the scalar level or in mixed precision block form. We also discuss a new probabilistic approach to rounding error analysis that provides error constants that are the square roots of those of the worst-case bounds. Combining these different considerations provides new understanding of the numerical stability of extreme scale and low precision computations in numerical linear algebra
Computational graphs for matrix functions
Full text linkMany numerical methods for evaluating matrix functions can be naturally viewed as computational graphs. Rephrasing these methods as directed acyclic graphs (DAGs) is a particularly effective way to study existing techniques, improve them, and eventually derive new ones. As the accuracy of these matrix techniques is determined by the accuracy of their scalar counterparts, the design of algorithms for matrix functions can be viewed as a scalar-valued optimization problem. The derivatives needed during the optimization can be calculated automatically by exploiting the structure of the DAG, in a fashion akin to backpropagation. The Julia package GraphMatFun.jl offers the tools to generate and manipulate computational graphs, to optimize their coefficients, and to generate Julia, MATLAB, and C code to evaluate them efficiently. The software also provides the means to estimate the accuracy of an algorithm and thus obtain numerically reliable methods. For the matrix exponential, for example, using a particular form (degree-optimal) of polynomials produces algorithms that are cheaper, in terms of computational cost, than the Padé-based techniques typically used in mathematical software. The optimized graphs and the corresponding generated code are available online
A spectral-in-time Newton-Krylov method for nonlinear PDE-constrained optimization
Full text linkWe devise a method for nonlinear time-dependent PDE-constrained optimization problems that uses a spectral-in-time representation of the residual, combined with a Newton-Krylov method to drive the residual to zero. We also propose a preconditioner to accelerate this scheme. Numerical results indicate that this method can achieve fast and accurate solution of nonlinear problems for a range of mesh sizes and problem parameters
Generating extreme-scale matrices with specified singular values or condition numbers
Full text linkA widely used form of test matrix is the randsvd matrix constructed as the product A = USV*, where U and V are random orthogonal or unitary matrices from the Haar distribution and S is a diagonal matrix of singular values. Such matrices are random but have a specified singular value distribution. The cost of forming an m-by-n randsvd matrix is m³ + n³ flops, which is prohibitively expensive at extreme scale; moreover, the randsvd construction requires a significant amount of communication, making it unsuitable for distributed memory environments. By dropping the requirement that U and V be Haar distributed and that both be random, we derive new algorithms for forming A that have cost linear in the number of matrix elements and require a low amount of communication and synchronization. We specialize these algorithms to generating matrices with specified 2-norm condition number. Numerical experiments show that the algorithms have excellent efficiency and scalability
A Catalogue of Software for Matrix Functions. Version 3.0
Full text linkA catalogue of software for computing matrix functions and their Fr\'echet derivatives is presented. For a wide variety of languages and for both open source packages and commercial products we describe what matrix function codes are available and which algorithms they implement