MIMS EPrints
Not a member yet
2151 research outputs found
Sort by
Improving the Complexity of Block Low-Rank Factorizations with Fast Matrix Arithmetic
No full textWe consider the LU factorization of an matrix represented as a block low-rank (BLR) matrix: most of its off-diagonal blocks are approximated by matrices of small rank , which reduces the asymptoticcomplexity of computing the LU factorization of down to \O(n^2r). In this article, our aim is to further reduce this complexity by exploiting fast matrix arithmetic, that is, the ability to multiply two full-rank matrices together for \O(n^\w) flops, where \w<3. This is not straightforward: simply accelerating the intermediate operations performed in the standard BLR factorization algorithm does not suffice to reduce the quadratic complexity in , because these operations are performed on matrices whose size is too small. To overcome this obstacle, we devise a new BLR factorization algorithm that, by recasting the operations so as to work on intermediate matrices of larger size, can exploit more efficiently fast matrix arithmetic. This new algorithm achieves an asymptotic complexity of \O(n^{(\w+1)/2}r^{(\w-1)/2}), which represents an asymptotic improvement compared to the standard BLR factorization as soon as \w<3. In particular, for Strassen's algorithm, \w\approx2.81 yields a complexity of \O(n^{1.904}r^{0.904}). Our numerical experiments are in good agreement with this theoretical result
Backward error and condition number of a generalized Sylvester equation, with application to the stochastic Galerkin method
Full text linkThe governing equation of the stochastic Galerkin method can be formulated as a generalized
Sylvester equation. Therefore, developing solvers for it is attracting a lot of attention from
the uncertainty quantification community. In this regard Krylov subspace based iterative
solvers, that are used for standard linear systems are being used for the generalized Sylvester
equations as well. This is achieved by converting the generalized Sylvester equation to a
standard linear system using the Kronecker product. Accordingly the residual norm is used
as a stopping criterion for the iterations, and the condition number of linear systems is used
for the generalized Sylvester equations as well. For a linear system a small residual norm
implies a small backward error, and hence using residual norm as a stopping criterion is
justified. In this work we prove that this need not be the case for the generalized Sylvester
equation. We introduce two definitions for the backward error, and then derive an upper
bound on each of them. We also verify the predictions of the analysis using numerical
experiments. For the special case of the stochastic Galerkin method we show that the upper
bound on the backward error can be computed with minimal computational overhead, and
hence it can be used as a stopping criterion in the iterative solvers. For the matrices
stemming from the stochastic Galerkin method we numerically demonstrate that the actual
backward error can be up to 3 orders of magnitude higher than the relative residual. Finally
by taking into account the structure of the equation we derive an expression for the condition
number, and discuss an algorithm for its computation in the special case of the stochastic
Galerkin method
Singular Reduction of the 2-Body Problem on the 3-Sphere and the 4-Dimensional Spinning Top
Full text linkWe consider the dynamics and symplectic reduction of the 2-body problem on a sphere of arbitrary dimension. It suffices to consider the case for when the sphere is 3-dimensional. As the 3-sphere is a group it acts on itself by left and right multiplication, which together generate the action of the SO(4) symmetry. This gives rise to a notion of left and right momenta for the problem, and allows for a reduction in stages, first by the left and then the right, or vice versa. The intermediate reduced spaces obtained by left or right reduction are shown to be coadjoint orbits of the special Euclidean group SE(4). The full reduced spaces are generically 4-dimensional and we describe these spaces and their singular strata.
The dynamics of the 2-body problem descend through a double cover to give a dynamical system on SO(4) which, after reduction and for a particular choice of hamiltonian, coincides with that of a 4-dimensional spinning top with symmetry. This connection allows us to ‘hit two birds with one stone’ and derive results about both the spinning top and the 2-body problem simultaneously. We provide the equations of motion on the reduced spaces and fully classify the relative equilibria and discuss their stability
An Updated Set of Nonlinear Eigenvalue Problems
No full textThis is the description of the new changes made in version 4.0 of the NLEVP MATLAB toolbox. The collection and its organization are described in a separate paper.
A user's guide describes how to install and use the NLEVP MATLAB toolbox
The block rational Arnoldi method
Full text linkThe block version of the rational Arnoldi method is a widely used procedure for generating an orthonormal basis of a block rational Krylov space. We study block rational Arnoldi decompositions associated with this method and prove an implicit Q theorem. We relate these decompositions to nonlinear eigenvalue problems. We show how to choose parameters to prevent a premature breakdown of the method and improve its numerical stability. We explain how rational matrix-valued functions are encoded in rational Arnoldi decompositions and how they can be evaluated numerically. Two different types of deflation strategies are discussed. Numerical illustrations using the MATLAB Rational Krylov Toolbox are included
ABBA: Adaptive Brownian bridge-based symbolic aggregation of time series
Full text linkA new symbolic representation of time series, called ABBA, is introduced. It is based on an adaptive polygonal chain approximation of the time series into a sequence of tuples, followed by a mean-based clustering to obtain the symbolic representation. We show that the reconstruction error of this representation can be modelled as a random walk with pinned start and end points, a so-called Brownian bridge. This insight allows us to make ABBA essentially parameter-free, except for the approximation tolerance which must be chosen. Extensive comparisons with the SAX and 1d-SAX representations are included in the form of performance profiles, showing that ABBA is able to better preserve the essential shape information of time series at compression rates comparable to other algorithms. A Python implementation is provided
Accurately Computing the Log-Sum-Exp and Softmax Functions
Full text linkEvaluating 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 \new{and shifting can introduce subtractive cancellation}. 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, so can safely be used, but that a division-free variant of softmax can suffer from loss of accuracy
Mixed Precision Block Fused Multiply-Add: Error Analysis and Application to GPU Tensor Cores
Full text linkTensor Tomography
Full text linkRich tomography is becoming increasingly popular since we have seen a substantial increase in computational power and storage. Instead of measuring one scalar for each ray, multiple measurements are needed per ray for various imaging modalities. This advancement has allowed the design of experiments and equipment which facilitate a
broad spectrum of applications. We present new reconstruction results and methods for several imaging modalities including x-ray diffraction strain tomography, Photoelastic tomography and Polarimetric Neutron Magnetic Field Tomography (PNMFT). We begin with a survey of the
Radon and x-ray transforms discussing several procedures for inversion. Furthermore we highlight the Singular Value Decomposition (SVD) of the Radon transform and consider some stability results for reconstruction in Sobolev spaces.
We then move onto define the Non-Abelian Ray Transform (NART), Longitudinal Ray Transform (LRT), Transverse Ray Transform (TRT) and the Truncated Transverse Ray Transform (TTRT) where we highlight some results on the complete
inversion procedure, SVD and mention stability results in Sobolev spaces. Thereafter we derive some relations between these transforms. Next we discuss the imaging modalities in mind and relate the transforms to their specific inverse problems, primarily being linear. Specifically, NART arises in the formulation of PNMFT where we want to image magnetic structures within magnetic materials with the use of polarized neutrons.
After some initial numerical studies we extend the known Radon inversion presented by experimentalists, reconstructing fairly weak magnetic fields, to reconstruct PNMFT data up to phase wrapping. We can recover the strain field tomographically for a polycrystalline material using
diffraction data and deduce that a certain moment of that data corresponds to the TRT. Quite naturally the whole strain tensor can be reconstructed from diffraction
data measured using rotations about six axes. We develop an innovative explicit plane-by-plane filtered back-projection reconstruction algorithm for the TRT, using data from rotations about three orthogonal axes and state the reasoning why two-axis data is insufficient. For the first time we give the first published results of TRT
reconstruction. To complete our discussion we present Photoelastic tomography which relates to the TTRT and implement the algorithm discussing the difficulties that arise in reconstructing data.
Ultimately we return to PNMFT highlighting the nonlinear inverse problem due to phase wrapping. We propose an iterative reconstruction algorithm, namely the
Modified Newton Kantarovich method (MNK) where we keep the Jacobian (Fréchet derivative) fixed at the first step. However, this is shown to fail for large angles
suggesting to develop the Newton Kantarovich (NK) method where we update the Jacobian at each step of the iteration process