Computing the Wave-Kernel Matrix Functions

The wave-kernel functions $\cosh{\sqrt{A}}$ and $\mathrm{sinhc}{\sqrt{A}}$ arise in the solution of second order differential equations such as $u''(t) - Au(t) = b(t)$ with initial conditions at $t=0$. Here, $A$ is an arbitrary square matrix and $\mathrm{sinhc}{z} = \sinh(z)/z$. The square root in these formulas is illusory, as both functions can be expressed as power series in $A$, so there are no questions about existence of the functions.

How can these functions be computed efficiently? In Computing the Wave-Kernel Matrix Functions (SIAM J. Sci. Comput., 2018) Prashanth Nadukandi and I develop an algorithm based on Padé approximation and the use of double angle formulas. The amount of scaling and the degree of the Padé approximant are chosen to minimize the computational cost subject to achieving backward stability for $\cosh{\sqrt{A}}$ in exact arithmetic.

In the derivation we show that the backward error of any approximation to $\cosh{\sqrt{A}}$ can be explicitly expressed in terms of a hypergeometric function. To bound the backward error we derive and exploit a new bound for $\|A^k\|^{1/k}$ in terms of the norms of lower powers of $A$; this bound is sharper than one previously obtained by Al-Mohy and Higham.

Numerical experiments show that the algorithm behaves in a forward stable manner in floating-point arithmetic and is superior in this respect to the general purpose Schur–Parlett algorithm applied to these functions.

The fundamental regions of the

function cosh(sqrt(z)), needed for the backward error analysis

underlying the algorithm.

Our Alumni – Lijing Lin

In this blog post, we asked one of our alumni, Lijing Lin, a few questions about her time with the Numerical Linear Algebra Group.

Please can you introduce yourself and tell us a bit about your experience before attending University of Manchester?

I obtained my BSc from Nanjing University of Aeronautics and Astronautics and MSc from Fudan University in China, before coming to Manchester to study for my PhD in 2007.

What was your PhD thesis on?

The title of my thesis is Roots of Stochastic Matrices and Fractional Matrix Powers. Computing roots of stochastic matrices arises from Markov chain models in finance and healthcare where a transition over a certain time interval is needed but only a transition over a longer time interval may be available. Besides developing new theories, we also developed a package for computing stochastic roots. Fractional matrix powers are more general functions than matrix roots. We developed a new algorithm for computing arbitrary real powers of matrices.

Why did you choose to study your PhD in Manchester?

I had developed an interest in doing research in Numerical Linear Algebra during my MSc. The NLA group in Manchester is renowned for world-leading expertise in this area, and is one of the best places in the world to study and do research.

How did you find Manchester?

I have studied, worked and lived in Manchester for over 11 years now. It is exciting, diverse and welcoming–a city that keeps growing and never stops surprising me.

After graduating, I continued working in Manchester as a Research Associate. With a solid background in NLA, my research now has moved toward machine learning, probabilistic modelling, and statistics.

I am currently a Turing PDRA in predictive healthcare. We are building prognostic models that allow consideration of “what if” scenarios to explore the effects of interventions, e.g. how would a person’s risk of getting heart attack change if he started or quit smoking now.

Celebrating the Centenary of James H. Wilkinson’s Birth

by Sven Hammarling and Nick Higham

September 27, 2019 is the 100th anniversary of the birth of James Hardy Wilkinson—the renowned numerical analyst who died in 1986. We are marking this special anniversary year in several ways:

The tag wilkinson lists all the posts in this series.

Jack Dongarra Awarded SIAM/ACM Prize in Computational Science and Engineering

Congratulations to Jack Dongarra who recently received the SIAM/ACM Prize in Computational Science and Engineering.

Jack Dongarra will receive the SIAM/ACM Prize in Computational Science and Engineering at the SIAM Conference on Computational Science and Engineering (CSE19) held February 25 – March 1, 2019 in Spokane, Washington. He will receive the award and deliver his prize lecture, “The Singular Value Decomposition: Anatomy of an Algorithm, Optimizing for Performance,” on February 28, 2019.

SIAM and the Association for Computing Machinery (ACM) jointly award the SIAM/ACM Prize in Computational Science and Engineering every two years at the SIAM Conference on Computational Science and Engineering for outstanding contributions to the development and use of mathematical and computational tools and methods for the solution of science and engineering problems. With this award, SIAM and ACM recognize Dongarra for his key role in the development of software and software standards, software repositories, performance and benchmarking software, and in community efforts to prepare for the challenges of exascale computing, especially in adapting linear algebra infrastructure to emerging architectures.

When asked about his research for which the prize was awarded, Dongarra said “I have been involved in the design and development of high performance mathematical software for the past 35 years, especially regarding linear algebra libraries for sequential, parallel, vector, and accelerated computers. Of course, the work that led to this award could not have been achieved without the help, support, collaboration, and interactions of many people over the years. I have had the good fortune of working on a number of high profile projects: in the area of mathematical software, EISPACK, LINPACK, LAPACK, ScaLAPACK, ATLAS and today with PLASMA, MAGMA, and SLATE; community de facto standards such as the BLAS, MPI, and PVM; performance analysis and benchmarking tools such as the PAPI, LINPACK benchmark, the Top500, and HPCG benchmarks; and the software repository netlib, arguably the first open source repository for publicly available mathematical software.”

This article was extracted from SIAM News. Further information is available here.

Professor Jack Dongarra

Numerical Algorithms for High-Performance Computational Science, Royal Society, London, April 8-9, 2019

This 2-day scientific discussion meeting is being held at the Royal Society, London, UK, and is organized by Jack Dongarra, Laura Grigori and Nick Higham.

The programme consists of invited talks and contributed posters. The invited speakers are

1. Guillaume Aupy, INRIA Bordeaux
2. Erin Carson, Charles University, Prague
3. George Constantinides, Imperial College
4. Steve Furber, FRS, University of Manchester
5. Mike Heroux, Sandia National Laboratories
6. Tony Hey, Science and Technology Facilities Council
7. David Keyes, King Abdullah University of Science and Technology, Saudi Arabia
8. Doug Kothe, Director of the DOE Exascale Computing Project (ECP)
9. Satoshi Matsuoka, RIKEN Center for Computational Science, Tokyo
10. Tim Palmer, FRS, Oxford University
11. Jack Poulson, Hodge Star Scientific Computing, Toronto
12. Anna Scaife, University of Manchester
13. John Shalf, Lawrence Berkeley National Laboratory
14. Rick Stevens, Argonne National Laboratory
15. Michela Taufer, University of Delaware
16. Kathy Yelick, Lawrence Berkeley National Laboratory,

The deadline for submission of poster abstracts is Friday March 1, 2019 (extended from Monday 4 February 2019).

Attendance is free, but places are limited and advance registration is essential.

A new preconditioner exploiting low-rank factorization error

The solution of a linear system $Ax = b$ is a fundamental task in scientific computing. Two main classes of methods to solve such a system exist.

• Direct methods compute a factorization of matrix $A$, such as LU factorization, to then directly obtain the solution $x=U^{-1}L^{-1}b$ by triangular substitution; they are very reliable but also possess a high computational cost, which limits the size of problems that can be tackled.
• Iterative methods compute a sequence of iterates $x_k$ converging towards the solution $x$; they are inexpensive but their convergence and thus reliability strongly depends on the matrix properties, which limits the scope of problems that can be tackled.

A current major challenge in the field of numerical linear algebra is therefore to develop methods that are able to tackle a large scope of problems of large size.

To accelerate the convergence of iterative methods, one usually uses a preconditioner, that is, a matrix $M$ ideally satisfying three conditions: (1) $M$ is cheap to compute; (2) $M$ can be easily inverted; (3) $M^{-1}$ is a good approximation to $A^{-1}$. With such a matrix $M$, the preconditioned system $M^{-1}Ax=M^{1}b$ is then cheap to solve with an iterative method and often requires a small number of iterations only. An example of a widely used class of preconditioners is when $M$ is computed as a low-accuracy LU factorization.

Unfortunately, for many important problems it is quite difficult to find a preconditioner that is both of good quality and cheap to compute, especially when the matrix $A$ is ill conditioned, that is, when the ratio between its largest and smallest singular values is large.

In our paper A New Preconditioner that Exploits Low-rank Approximations to Factorization Error, with Nick Higham, which recently appeared in SIAM Journal of Scientific Computing, we propose a novel class of general preconditioners that builds on an existing, low-accuracy preconditioner $M=A-\Delta A$.

This class of preconditioners is based on the following key observation: ill-conditioned matrices that arise in practice often have a small number of small singular values. The inverse of such a matrix has a small number of large singular values and so is numerically low rank. This observation suggests that the error matrix $E = M^{-1}A - I = M^{-1}\Delta A \approx A^{-1}\Delta A$ is of interest, because we may expect $E$ to retain the numerically low-rank property of $A^{-1}$.

In the paper, we first investigate theoretically and experimentally whether $E$ is indeed numerically low rank; we then describe how to exploit this property to accelerate the convergence of iterative methods by building an improved preconditioner $M(I+\widetilde{E})$, where $\widetilde{E}$ is a low-rank approximation to $E$. This new preconditioner is equal to $A-M(E-\widetilde{E})$ and is thus almost a perfect preconditioner if $\widetilde{E}\approx E$. Moreover, since $\widetilde{E}$ is a low-rank matrix, $(I+\widetilde{E})^{-1}$ can be cheaply computed via the Sherman–Morrison–Woodbury formula, and so the new preconditioner can be easily inverted.

We apply this new preconditioner to three different types of approximate LU factorizations: half-precision LU factorization, incomplete LU factorization (ILU), and block low-rank (BLR) LU factorization. In our experiments with GMRES-based iterative refinement, we show that the new preconditioner can achieve a significant reduction in the number of iterations required to solve a variety of real-life $Ax=b$ problems.

Research Associate in Numerical Linear Algebra Group – 2 posts

The Numerical Linear Algebra Group at the University of Manchester is seeking two Research Associates to work with Professor Nick Higham on developing and analyzing numerical linear algebra algorithms for current and future high-performance computers.

Topics for investigation include linear equations, linear least squares problems, eigenvalue problems, the singular value decomposition, correlation matrix problems, and matrix function evaluation. This work will exploit multiprecision arithmetic (particularly the fast half precisions available on some recent and forthcoming processors) and techniques such as acceleration and randomization. It will involve using rounding error analysis, statistical analysis, and numerical experiments to obtain new understanding of algorithm accuracy and efficiency.

One of the posts is associated with Professor Higham’s Royal Society Research Professorship. The other post is associated with the EPSRC project Inference, Computation and Numerics for Insights into Cities (ICONIC), which involves Imperial College (Professor Mark Girolami), the University of Manchester (Professor Nick Higham), the University of Oxford (Professor Mike Giles), and the University of Strathclyde (Professor Des Higham).

The closing date is February 11, 2019. For the advert and more details see here.