## Solving Block Low-Rank Linear Systems by LU Factorization is Numerically Stable

In many applications requiring the solution of a linear system $Ax=b$, the matrix $A$ possesses a block low-rank (BLR) property: most of its off-diagonal blocks are of low numerical rank and can therefore be well approximated by low-rank matrices. This property arises for example in the solution of discretized partial differential equations, because the numerical rank of a given off-diagonal block is closely related to the distance between the two subdomains associated with this block. BLR matrices have been exploited to significantly reduce the cost of solving $Ax=b$, in particular in the context of sparse direct solvers such as MUMPS, PaStiX, and STRUMPACK.

However, the impact of these low-rank approximations on the numerical stability of the solution in floating-point arithmetic has not been previously analyzed. The difficulty of such an analysis lies in the fact that, unlike for classical algorithms without low-rank approximations, there are two kinds of errors to analyze: floating-point errors (which depend on the unit roundoff $u$), and low-rank truncation errors (which depend on the the low-rank threshold $\varepsilon$, the parameter controlling the accuracy of the blockwise low-rank approximations). Moreover, the two kinds of error cannot easily be isolated in the analysis, because the BLR compression and factorization stages are often interlaced.

In our recent preprint, with Nick Higham, we present rounding error analysis for the solution of a linear system by LU factorization of BLR matrices. Assuming that a stable pivoting scheme is used, we prove backward stability: the relative backward error is bounded by a modest constant times $\varepsilon$, and does not depend on the unit roundoff $u$ as long as $u$ is safely smaller than $\varepsilon$. This is a very desirable theoretical guarantee that can now be given to the users, who can therefore control the numerical behavior of BLR solvers simply by setting $\varepsilon$ to the target accuracy. We illustrate this key result in the figure below, which shows a strong correlation between the backward error and the $\varepsilon$ threshold for 26 matrices from the SuiteSparse collection coming from a wide range of real-life applications.

## Nick Higham Delivers Invited Talk at the ICIAM Congress in Valencia

Photo provided by @ICIAMnews

Professor Nick Higham delivered an invited talk at the International Congress on Industrial and Applied Mathematics (ICIAM) on July 19, 2019, in Valencia, Spain. His talk “Exploiting Low Precision Arithmetic in the Solution of Linear Systems” reported work by Higham and his colleagues over the last three years to use the fast half precision arithmetic available on accelerators such as GPUs to speed up the solution of linear systems.  The slides are available here.

The talk was summarized in the Friday 19th July edition of the ICIAM2019 newsletter:

“Nick Higham (University of Manchester) advocated in his invited lecture using arithmetics of different precision at different stages of computations in order to design algorithms that are faster, require less communications and consume less energy. The main motivation is that last-generation GPUs may be up to eight times faster when they perform arithmetic operations in half precision than when they do in single precision.

The general philosophy is doing the bulk of the computations in half precision, and then perform some kind of clean-up refinement of the solution in higher precisions. As an archetypal example to illustrate this line of thought Higham chose how to accelerate the solution of linear systems via Gaussian elimination. He showed that if the LU factorization is computed in half precision, followed by iterative refinement using a mix of half, double and quadruple precision, the solution can be sped up significantly without accuracy loss. This works in principle for systems with moderate condition number, but even ill-conditioned systems can be dealt with by performing the iterative refinement via GMRES on a suitably pre-conditioned system.”

Photo credit Françoise Tisseur

## Françoise Tisseur Delivers the Olga Taussky-Todd Lecture at the ICIAM Congress in Valencia

Professor Françoise Tisseur delivered the Olga Taussky-Todd Lecture at the International Congress on Industrial and Applied Mathematics (ICIAM) on July 15, 2019, in Valencia, Spain. The honour of giving this lecture is conferred every four years on a woman who has made outstanding contributions in applied mathematics and/or scientific computation.

Her talk “Challenges in the numerical solution of nonlinear eigenvalue problems” (slides available here in PDF form) concerned eigenvalue problems defined in terms of matrix-valued functions $F$ of a scalar parameter, which occur in many areas of science and engineering. She described recent advances in numerical methods for solving such problems, which include methods for approximating $F$ by a rational or polynomial matrix function and methods for linearizing rational and polynomial eigenproblems. Pitfalls that the methods must overcome were illustrated through simple but realistic practical examples.

Her talk followed the opening ceremony, which included traditional Valencian dancing and an address to the approximately 4000 attendees by King Felipe VI of Spain (video clip), who spoke of the good health of Spanish mathematics and the importance of mathematics for driving innovation.

## Massimiliano Fasi awarded PhD

Congratulations to Massimiliano Fasi for being awarded his PhD, which was supervised by Nick Higham. We asked him a few questions about his thesis, title Computing Matrix Functions in Arbitrary Precision Arithmetic.

Please can you introduce yourself and tell us a bit about your experience before attending University of Manchester?
I was born in Assisi, a smallish town in central Italy chiefly known for having been home to Saint Francis. I thought I would become a philologist, and during my teen years I studied mostly linguistics and ancient languages, while taking plenty of courses at the local music conservatory. After much pondering, I decided to start my university career with a scientific discipline and I chose computer science. I found that I really enjoyed the subject, and after graduating from the University of Perugia I pursued a Master’s degree at the University of Bologna and, at the same time, one at the École Normale Supérieure of Lyon.

My thesis investigates the computation of matrix functions in arbitrary precision arithmetic, and is in fact a collection of results that are somewhat connected to this topic. We started by revisiting several well-known techniques for two very general problems: the evaluation of rational functions at a matrix argument and the solution of rational matrix equations. We developed new numerical methods and provided some theoretical insights, and used these as building blocks to design multiprecision algorithms to evaluate matrix functions that are of interest in applications. We focused mainly on the matrix exponential and the matrix logarithm, but the machinery we developed is more general and can be used for a much broader range of problems.

Why did you choose to study your PhD in Manchester?
Numerical analysis was my favourite module as an undergraduate student. I decided to do an internship with the lecturer teaching that course, and the subject happened to be the matrix Lambert W function. I’ve been interested in matrix functions ever since, and when the time came to choose where to do a PhD, I thought that it would be great to join the strongest linear algebra group I was aware of. I was lucky enough to be admitted at The University of Manchester, where I could do research with Professor Nick Higham, who is a leading expert in this field.

How did you find Manchester?
Manchester is a lively city and is well connected to the rest of the world. The postgraduate community I was part of is very dynamic, and somewhat larger than I was used to. A PhD can be rather difficult at times, but the many people I had the chance to meet here made my Mancunian experience unforgettable.

My long-term plan is to remain in academia, and I will try to follow the path that leads to a permanent position. For the moment, I am a research associate in the group here in Manchester, and am trying my best to deepen my understanding of the current trends in numerical linear algebra.

Massimiliano Fasi

## Nick Higham awarded the LMS Naylor Prize and Lectureship

Professor Nick Higham has been awarded the London Mathematics Society’s Naylor Prize and Lectureship for his leadership in numerical linear algebra, numerical stability analysis, and communication of mathematics.

The Naylor Prize and Lectureship is awarded every two years in memory of Dr V. D. Naylor. The grounds for the award may include work in, and influence on, and contributions to applied mathematics and/or the applications of mathematics, and lecturing gifts.

The full prize citation is available here.

## Determinants of Bohemian matrix families

The adjective “Bohemian” was used for the first time in a linear algebra context by Robert Corless and Steven Thornton to describe the eigenvalues of matrices whose entries are taken from a finite discrete set, usually of integers. The term is a partial acronym for “BOunded Height Matrix of Integers”, but the origin of the term was soon forgotten, and the expression “Bohemian matrix” is now widely accepted.

As Olga Taussky observed already in 1960, the study of matrices with integer elements is “very vast and very old”, with early work of Sylvester and Hadamard that dates back to the second half of the nineteenth century. These names are the first two in a long list of mathematicians that worked on what is now known as the “Hadamard conjecture”: for any positive integer $n$ multiple of 4, there exists an $n$ by $n$ matrix $H$, with entries $-1$ and $+1$, such that $HH^T = nI$.

If this is the best-known open problem surrounding Bohemian matrices, it is far from being the only one. During the 3-day workshop “Bohemian Matrices and Applications” that our group hosted in June last year, Steven Thornton released the Characteristic Polynomial Database, which collects the determinants and characteristic polynomials of billions of samples from certain families of structured as well as unstructured Bohemian matrices. All the available data led Steven to formulate a number of conjectures regarding the determinants of several families of Bohemian upper Hessenberg matrices.

Gian Maria Negri Porzio and I attended the workshop, and set ourselves the task of solving at least one of these open problems. In our recent preprint, we enumerate all the possible determinants of Bohemian upper Hessenberg matrices with ones on the subdiagonal. We consider also the special case of families with main diagonal fixed to zero, whose determinants turn out to be related to some generalizations of Fibonacci numbers. Many of the conjectures stated in the Characteristic Polynomial Database follow from our results.

## Highlights of Advances in Numerical Linear Algebra Conference

by Sven Hammarling, Nick Higham and Françoise Tisseur

The conference Advances in Numerical Linear Algebra: Celebrating the Centenary of the Birth of James H. Wilkinson, took place at the University of Manchester, May 29-30, 2019.  The purpose of the conference was to discuss recent developments and future challenges in numerical linear algebra and to celebrate the centenary of the birth of James H. Wilkinson FRS (1919-1986).

The 60 attendees heard reminiscences about Wilkinson and his work from several attendees who knew him.  Sven Hammarling opened the proceedings with personal reflections on Wilkinson. Margaret Wright discussed some of the treasures in lecture notes from courses that Wilkinson taught at Stanford University (1977-1982). We have recently added these notes to our Wilkinson web page. Nick Higham announced the availability of the Argonne tapes, which are videos of Wilkinson and Cleve Moler lecturing at an Eigensystem Workshop held at Argonne National Laboratory, Illinois, USA, in June 1973.

We were pleased to welcome two of Jim Wilkinson’s nephews, John Liebman and Danny Liebman, together with John’s wife Liz.

Attendees were able to wish speaker Cleve Moler a happy 80th birthday and share two birthday cakes with him.

Photos from the conference appear below. Click on a photo to enlarge it.

Slides of the talks are available at the conference web page. The first day’s lectures were professionally videoed and the videos are available on the NLA group’s YouTube channel: here is a link to the playlist.

Thank you to all the speakers and attendees for their participation.

We gratefully acknowledge sponsorship from the Royal Society, The Alan Turing Institute, The QJMAM Fund for Applied Mathematics, The Numerical Algorithm Group and National Physical Laboratory.