MIMS EPrints
Not a member yet
2151 research outputs found
Sort by
A Proposed API for Batched Basic Linear Algebra Subprograms
Get PDFThis paper proposes an API for Batched Basic Linear Algebra Subprograms (Batched BLAS).
We focus on many independent BLAS operations on small matrices that are grouped together as a single routine, called Batched BLAS routine, with the aim of providing more efficient, but portable, implementations of algorithms on high-performance manycore architectures
(like multi/manycore CPU processors, GPUs, and coprocessors)
Theory and Algorithms for Periodic Functions of Matrices, with Applications
Get PDFTheoretical aspects of periodic functions of matrices and issues arising from the multivalued nature of their inverse functions are studied. Several algorithms for computing periodic and multivalued functions of matrices are developed.
We illustrate the use of matrix functions in the analysis of complex networks---an application that has recently been of very high interest. The relative importance of nodes in the whole network can be expressed via functions of the adjacency matrix. There are two functions, which have proven popular in practice. The first one is the exponential, which has the advantage of being parameter-free. The second one is the resolvent function, which can be the more computationally efficient, but it depends on a parameter. We give a prescription for selecting this parameter aiming to match the rankings of the exponential counterpart.
We define a new matrix function, the matrix unwinding function, corresponding to the scalar unwinding number of Corless, Hare, and Jeffrey introduced in 1996. The matrix unwinding function is shown to be an important tool for deriving identities involving the matrix logarithm and fractional matrix powers. We propose an algorithm for computing the matrix unwinding function based on the Schur--Parlett method with a special reordering. The matrix unwinding function is shown to be useful for computing the matrix exponential using an idea of argument reduction.
We study theoretical and computational aspects of matrix inverse trigonometric and inverse hyperbolic functions. Conditions for existence are given and principal values are defined and shown to be unique primary matrix functions. We derive various functional identities, with care taken to specify choices of signs and branches. An important tool for the derivations is the matrix unwinding function. We derive a new algorithm employing a Schur decomposition and variable-degree rational approximation for computing the principal inverse cosine (acos). It is shown how it can also be used to compute the matrix asin, acosh, and asinh. In numerical experiments the algorithm is found to behave in a forward stable fashion.
Finally, we consider argument reduction in computing the sine and cosine, and their hyperbolic counterparts. New algorithms for these functions are given, which use the matrix unwinding function with multiple angle algorithms for the sine and cosine. An argument reduction algorithm for computing general periodic functions of matrices is presented. Numerical experiments illustrate the computational saving that can accrue from applying argument reduction
BETTER IMAGING FOR LANDMINE DETECTION: AN EXPLORATION OF 3D FULL-WAVE INVERSION FOR GROUND-PENETRATING RADAR
Get PDFHumanitarian clearance of minefields is most often carried out by hand, conventionally using a a metal detector and a probe. Detection is a very slow process, as every piece of detected metal must treated as if it were a landmine and carefully probed and excavated, while many of them are not. The process can be safely sped up by use of Ground-Penetrating Radar (GPR) to image the subsurface, to verify metal detection results and safely ignore any objects which could not possibly be a landmine. In this thesis, we explore the possibility of using Full Wave Inversion (FWI) to improve GPR imaging for landmine detection. Posing the imaging task as FWI means solving the large-scale, non-linear and ill-posed optimisation problem of determining the physical parameters of the subsurface (such as electrical permittivity) which would best reproduce the data. This thesis begins by giving an overview of all the mathematical and implementational aspects of FWI, so as to provide an informative text for both mathematicians (perhaps already familiar with other inverse problems) wanting to contribute to the mine detection problem, as well as a wider engineering audience (perhaps already working on GPR or mine detection) interested in the mathematical study of inverse problems and FWI.We present the first numerical 3D FWI results for GPR, and consider only surface measurements from small-scale arrays as these are suitable for our application. The FWI problem requires an accurate forward model to simulate GPR data, for which we use a hybrid finite-element boundary-integral solver utilising first order curl-conforming N\'{e}d\'{e}lec (edge) elements. We present a novel `line search' type algorithm which prioritises inversion of some target parameters in a region of interest (ROI), with the update outside of the area defined implicitly as a function of the target parameters. This is particularly applicable to the mine detection problem, in which we wish to know more about some detected metallic objects, but are not interested in the surrounding medium. We may need to resolve the surrounding area though, in order to account for the target being obscured and multiple scattering in a highly cluttered subsurface.We focus particularly on spatial sensitivity of the inverse problem, using both a singular value decomposition to analyse the Jacobian matrix, as well as an asymptotic expansion involving polarization tensors describing the perturbation of electric field due to small objects. The latter allows us to extend the current theory of sensitivity in for acoustic FWI, based on the Born approximation, to better understand how polarization plays a role in the 3D electromagnetic inverse problem. Based on this asymptotic approximation, we derive a novel approximation to the diagonals of the Hessian matrix which can be used to pre-condition the GPR FWI problem
A rational deferred correction approach to PDE-constrained optimization
Get PDFThe accurate and efficient solution of time-dependent PDE-constrained optimization problems is a challenging task, in large part due to the very high dimension of the matrix systems that need to be solved. We devise a new deferred correction method for coupled systems of time-dependent PDEs, allowing one to iteratively improve the accuracy of low-order time stepping schemes. We consider two variants of our method, a splitting and a coupling version, and analyze their convergence properties. We then test our approach on a number of PDE-constrained optimization problems. We obtain solution accuracies far superior to that achieved when solving a single discretized problem, in particular in cases where the accuracy is limited by the time discretization. Our approach allows for the direct reuse of existing solvers for the resulting matrix systems, as well as state-of-the-art preconditioning strategies
Refined saddle-point preconditioners for discretized Stokes problems
Get PDFThis paper is concerned with the implementation of
efficient solution algorithms for elliptic problems with constraints. We establish theory that shows that including a simple scaling within well-established block diagonal preconditioners for Stokes problems can result in significantly faster convergence when applying the preconditioned MINRES method. The codes used in the numerical studies are available online
Less is More I: a pessimistic view of piecewise smooth bifurcation theory
Get PDFThe analysis of piecewise smooth bifurcations reveals an alarming proliferation of cases as the dimension of phase space increases. This suggests that a different approach needs to be taken when trying to describe bifurcations. In particular, it may not be helpful to analyze particular bifurcations at the level of detail that is standard for smooth systems
Theory and algorithms for matrix problems with positive semidefinite constraints
Get PDFThis thesis presents new theoretical results and algorithms for two matrix problems with positive semidefinite constraints: it adds to the well-established nearest correlation matrix problem, and introduces a class of semidefinite Lagrangian subspaces.
First, we propose shrinking, a method for restoring positive semidefiniteness of an indefinite matrix that computes the optimal parameter \a_* in a convex combination of and a chosen positive semidefinite target matrix. We describe three algorithms for computing \a_*, and then focus on the case of keeping fixed a positive semidefinite leading principal submatrix of an indefinite approximation of a correlation matrix, showing how the structure can be exploited to reduce the cost of two algorithms. We describe how weights can be used to construct a natural choice of the target matrix and that they can be incorporated without any change to computational methods, which is in contrast to the nearest correlation matrix problem. Numerical experiments show that shrinking can be at least an order of magnitude faster than computing the nearest correlation matrix and so is preferable in time-critical applications.
Second, we focus on estimating the distance in the Frobenius norm of a symmetric matrix to its nearest correlation matrix \Ncm(A) without first computing the latter. The goal is to enable a user to identify an invalid correlation matrix relatively cheaply and to decide whether to revisit its construction or to compute a replacement. We present a few currently available lower and upper bounds for \dcorr(A) = \normF{A - \Ncm(A)} and derive several new upper bounds, discuss the computational cost of all the bounds, and test their accuracy on a collection of invalid correlation matrices. The experiments show that several of our bounds are well suited to gauging the correct order of magnitude of \dcorr(A), which is perfectly satisfactory for practical applications.
Third, we show how Anderson acceleration can be used to speed up the convergence of the alternating projections method for computing the nearest correlation matrix, and that the acceleration remains effective when it is applied to the variants of the nearest correlation matrix problem in which specified elements are fixed or a lower bound is imposed on the smallest eigenvalue. This is particularly significant for the nearest correlation matrix problem with fixed elements because no Newton method with guaranteed convergence is available for it. Moreover, alternating projections is a general method for finding a point in the intersection of several sets and this appears to be the first demonstration that these methods can benefit from Anderson acceleration.
Finally, we introduce semidefinite Lagrangian subspaces, describe their connection to the unique positive semidefinite solution of an algebraic Riccati equation, and show that these subspaces can be represented by a subset and a Hermitian matrix that is a generalization of a quasidefinite matrix. We further obtain a semidefiniteness-preserving version of an optimization algorithm introduced by Mehrmann and Poloni [\textit{SIAM J.\ Matrix Anal.\ Appl.}, 33(2012), pp.\ 780--805] to compute a pair (\mathcal{I}_{\opt},X_{\opt}) with M = \max_{i,j} \abs{(X_{\opt})_{ij}} as small as possible, which improves numerical stability in several contexts
Near-optimal perfectly matched layers for indefinite Helmholtz problems
Get PDFA new construction of an absorbing boundary condition for indefinite Helmholtz problems on unbounded domains is presented. This construction is based on a near-best uniform rational interpolant of the inverse square root function on the union of a negative and positive real interval, designed with the help of a classical result by Zolotarev. Using Krein's interpretation of a Stieltjes continued fraction, this interpolant can be converted into a three-term finite difference discretization of a perfectly matched layer (PML) which converges exponentially fast in the number of grid points. The convergence rate is asymptotically optimal for both propagative and evanescent wave modes. Several numerical experiments and illustrations are included
Restoring Definiteness via Shrinking, with an Application to Correlation Matrices with a Fixed Block
Get PDFIndefinite estimates of positive semidefinite matrices arise in many data analysis applications involving covariance matrices and correlation matrices. We develop a method for restoring positive semidefiniteness of an indefinite estimate based on the process of shrinking, which finds a convex linear combination of the original matrix and a target positive semidefinite matrix . We describe three \alg s for computing the optimal shrinking parameter \alpha_* = \min \{\alpha \in [0,1] : \mbox{S(\alpha) is positive semidefinite}\}. One algorithm is based on the bisection method, with the use of Cholesky factorization to test definiteness, a second employs Newton's method, and a third finds the smallest eigenvalue of a symmetric definite generalized eigenvalue problem. We show that weights that reflect confidence in the individual entries of can be used to construct a natural choice of the target matrix . We treat in detail a problem variant in which a positive semidefinite leading principal submatrix of remains fixed, showing how the fixed block can be exploited to reduce the cost of the bisection and generalized eigenvalue methods. Numerical experiments show that when applied to estimates of correlation matrices shrinking can be at least an order of magnitude faster than computing the nearest correlation matrix
Block Kronecker Linearizations of Matrix Polynomials and their Backward Errors
Get PDFWe introduce a new family of strong linearizations of matrix polynomials---which we call "block Kronecker pencils"---and perform a backward stability analysis of complete polynomial eigenproblems. These problems are solved by applying any backward stable algorithm to a block
Kronecker pencil, such as the staircase algorithm for singular pencils or the QZ algorithm for regular pencils. This analysis allows us to identify those block Kronecker pencils that yield a computed complete eigenstructure which is exactly that of a slightly perturbed matrix polynomial. These favorable pencils include the famous Fiedler linearizations, which are just a very particular case of block Kronecker pencils. Thus, our analysis offers the first proof available in the literature of global backward stability for Fiedler pencils. In addition, the theory developed for block Kronecker pencils is much simpler than the theory available for Fiedler pencils, especially in the case of rectangular matrix polynomials. The global backward error analysis in this work presents for the first time the following key properties: it is a rigorous analysis valid for finite perturbations (i.e., it is not a first order analysis), it provides precise bounds, it is valid simultaneously for a large class of linearizations, and it establishes a framework that may be generalized to other classes of linearizations. These features
are related to the fact that block Kronecker pencils are a particular case of the new family of "strong block minimal bases pencils," which include certain perturbations of block Kronecker pencils; this allows us to extend the results in this paper to other contexts