## Numerical Linear Algebra Group Activities 2018

The Numerical Linear Algebra Group, many of whom are in the October 2018 photo below, had a busy year. This post summarizes what we got up to. Publications are not included here, but many of them can be found on MIMS EPrints under the category Numerical Analysis.

As well as the new group members mentioned below, Stephanie Lai joined the group in September as Executive Assistant to Nick Higham. Stephanie, who holds an MMath degree from the University of Manchester, helps to maintain this website, which launched this autumn. Martin Lotz left the group in August 2018 and is now an Associate Professor at The University of Warwick.

## Software

Together with Jack Dongarra’s team at the University of Tennessee the group is one of the two partners involved in the development of PLASMA: Parallel Linear Algebra Software for Multicore Architectures.

PhD students Weijian Zhang, Steven Elsworth and Jonathan Deakin continued to develop Etymo—a technology for developing knowledge graphs, including Etymo Scholar, a discovery engine for AI research.

We make our research codes available as open source, principally on GitHub; see the repositories of Connolly, Fasi, Higham, Liu, PraneshRelton, Sego, Tisseur, ZhangZounon.

We also put MATLAB software on MATLAB Central File Exchange and on our own web sites, e.g., the Rational Krylov Toolbox (RKToolbox).

## PhD Students

We welcomed new PhD students Michael Connolly, Xiaobo Liu, and Conor Rogers.

Weijian Zhang successfully defended his PhD on Evolving Graphs and similarity-based Graphs with Applications in October 2018.

## Postdoctoral Research Associates (PDRAs)

Theo Mary joined us in January 2018 to work jointly on the ICONIC project.

Mawussi Zounon, who was working on the Parallel Numerical Linear Algebra for Extreme Scale Systems (NLAFET) project, moved over in November 2018 to the position of KTP Associate in a three-year Knowledge Transfer Partnership with NAG.

Jakub Sistek returned to the Czech Academy of Sciences after the completion of the project Programming Model INTERoperability ToWards Exascale (INTERTWinE).

Prashanth Nadukandi completed his Marie Skłodowska-Curie Individual Fellowship in September 2018 and is now a Technical Advisor in Advanced Mathematics at the Repsol Technology Lab, Móstoles, Spain.

Maksims Abalenkovs completed his work on the EPSRC project SERT: Scale-free, Energy-aware, Resilient and Transparent Adaptation of CSE Applications to Mega-core Systems.

## Presentations at Conferences and Workshops

UK and Republic of Ireland Section of SIAM Annual Meeting, National Oceanography Centre, Southampton, January 11, 2018: Blanchard, Nadukandi.

Waseda Research Institute for Science and Engineering Institute for Mathematical Science Kickoff Meeting, Waseda University, Tokyo, March 6, 2018: Higham.

SIAM Conference on Parallel Processing for Scientific Computing, Waseda University, Tokyo, Japan, March 7-10, 2018: Hammarling, Higham, Mary, Zounon.

Conference on Scientific Computing and Approximation, Purdue University, March 30-31, 2018: Higham.

SIAM Conference on Applied Linear Algebra, May 4-8, 2018, Hong Kong Baptist University: Elsworth, Güttel, Zemaityte.

Structured Matrix Days 2018, May 14-15, 2018, ENS Lyon, France: Mary, Tisseur.

Mathematical Modelling, Numerical Analysis and Scientific Computing, Kácov, Czech Republic, May 27-June 1, 2018: Higham, Tisseur.

International Conference on Approximation and Matrix Functions, May 31 and June 1, 2018, Université de Lille, France: Güttel.

Bohemian Matrices and Applications, University of Manchester, June 20-22, 2018: Higham.

International Supercomputing conference (ISC18), Frankfurt, June 24-28, 2018: Hammarling, Zounon.

10th International Workshop on Parallel Matrix Algorithms and Applications (PMAA18), June 27-29, 2018, ETH Zurich: Zounon.

6th IMA Conference on Numerical Linear Algebra and Optimization, June 27-29 2018, University of Birmingham: Blanchard, Mary, Tisseur.

SIAM Annual Meeting, Portland, USA, July 9-13, 2018: Blanchard, Higham, Tisseur.

JuliaCon 2018, London, August 7-11, 2018: Higham.

Sparse Day Meeting 2018, September 26-27 2018, Toulouse, France: Mary.

## Recognition and Service

• Nick Higham gave the Samuel D. Conte Distinguished Lecture at the Department of Computer Science, Purdue University, USA, in March 2018.
• Massimilano Fasi was presented with the SIAM Student Chapter Certificate of Recognition 2018.

## The Paterson–Stockmeyer method for evaluating polynomials and rational functions of matrices

According to Moler and Van Loan, truncating the Taylor series expansion to the exponential at 0 is the least effective of their nineteen dubious ways to compute the exponential of a matrix. Such an undoubtedly questionable approach can, however, become a powerful tool for evaluating matrix functions when used in conjunction with other strategies, such as, for example, the scaling and squaring technique. In fact, truncated Taylor series are just a special case of a much larger class of rational approximants, known as the Padé family, on which many state-of-the-art algorithms rely in order to compute matrix functions.

A customary choice for evaluating these approximants at a matrix argument is the Paterson–Stockmeyer method, an evaluation scheme that was originally proposed as an asymptotically optimal algorithm for evaluating polynomials of matrices, but generalizes quite naturally to rational functions, which are nothing but the solutions of a linear system whose coefficients and right-hand side are polynomials of the same matrix. This technique exploits a clever rewriting of a polynomial in $A$ as a polynomials in $A^s$, for some positive integer $s$, and overall requires about $2\sqrt{k}$ matrix products to evaluate a polynomial of degree $k$.

As shown in the figure, when the Paterson–Stockmeyer scheme is used the number of matrix multiplications required to evaluate a polynomial of degree $k$ grows slower than $k$ itself, with the result that evaluating polynomials of different degree will asymptotically have the same cost. For example, evaluating a polynomial of any degree between 43 and 49 requires 12 matrix multiplications, thus there is little point in considering an approximant of degree 43 when evaluating that of degree 49 has roughly the same cost but will in all likelihood deliver a more accurate result.

When designing algorithms to evaluate functions of matrices, one is interested in finding the optimal degrees, those marked with a red circle in the figure above, since they guarantee maximal accuracy for a given computational cost. When fixed precision is considered, finding all such degrees is not a problem: a backward error analysis can be used to determine the maximum degree that will ever be needed, $m_{max}$ say, and then looking at the plot is enough to find all the optimal degrees smaller than $m_{max}$. In order to deal with arbitrary precision arithmetic, however, a different strategy is needed, as depending on the working precision and the required tolerance, approximants of arbitrarily high degree may be needed. The new Eprint Optimality of the Paterson–Stockmeyer Method for Evaluating Matrix Polynomials and Rational Matrix Functions studies the cost of the Paterson–Stockmeyer method for polynomial evaluation and shows that a degree is optimal if and only if it is a quarter-square, that is, a number of the form $\lfloor n^2/4 \rfloor$ for some nonnegative integer $n$, where $\lfloor \cdot \rfloor$ is the floor function.

Similar results can be obtained for Paterson–Stockmeyer-like algorithms for evaluating diagonal Padé approximants, rational functions whose numerator and denominator have same degree. In that case, one can show that an approximant is optimal if and only if the degree of numerator and denominator is an eight-square, an integer of the form $\lfloor n^2/8 \rfloor$ for some $n \in \mathbb{N}$. In particular, for diagonal Padé approximants to the exponential, due to a symmetry of the coefficients of numerator and denominator, faster algorithms can be developed, and an explicit formula—not as nice as that in the two previous cases—can be derived for the optimal orders of these approximants.

## Theo Mary awarded the Gilles Kahn prize

The Gilles Kahn prize was awarded to Theo Mary for his PhD thesis Block Low-Rank multifrontal solvers: complexity, performance, and scalability (pdf here). The thesis was prepared at the University of Toulouse, France and co-supervised by Patrick Amestoy and Alfredo Buttari. Theo is currently a Research Associate in the Numerical Linear Algebra group.

The Gilles Kahn prize is awarded each year by the SIF, the French Society of Computer Science, for an excellent PhD thesis in the field of computer science.

In the thesis, Theo investigated the use of low-rank approximations inside direct, factorization-based methods to reduce the computational cost of the solution of large, sparse systems of linear equations. His work included theoretically computing the asymptotic complexity reduction achieved by these methods and implementing them on parallel computers to translate this theoretical reduction into significant performance gains on both shared- and distributed-memory architectures. The potential and efficiency of these methods was assessed on several industrial applications, notably ones arising in geosciences and structural mechanics.

The official award ceremony will take place at the SIF congress (February  6-7, 2019, in Bordeaux). Theo has also been invited to present his work to the entire scientific community at the SIF congress.

## Weijian Zhang Awarded PhD for Work on Evolving Graphs

Congratulations to our PhD student Weijian Zhang for being awarded his PhD, which was supervised by Nick Higham. We asked him a few questions about his thesis.

Please can you introduce yourself and tell us a bit about your experience before attending University of Manchester?
I was born in Liaoning, China and studied fine arts at the Fine Arts High School Affiliated to China Central Academy of Fine Arts. I graduated with a First Class MMath degree at the University of Manchester in 2014.

My thesis title is “Evolving Graphs and similarity-based Graphs with Applications“. We look at ways to traverse an evolving graph. In particular, we examine the influence of temporal information on node centrality. In the process, we developed EvolvingGraphs.jl, a Julia software package for analysing time-dependent networks. We also developed Etymo Scholar, a search system for discovering interesting research papers.

Why did you choose to study your PhD in Manchester?
I was interested in doing something that relates to both maths and computer science. Numerical Linear Algebra looked like a perfect fit for me. I was very fortunate to be able to join Professor Nick Higham’s research group, which is a world-renowned research group in Numerical Linear Algebra.

How did you find Manchester?
Manchester is a vibrant city and the University of Manchester offers fantastic support for research students. All in all, it is a great place for doing research.

I have just joined Conversocial, a leading social customer service solution provider, as a Machine Learning Engineer / Data Scientist in London.

## Computing the matrix exponential in arbitrary precision

Edmond Laguerre was the first to discuss, en passant in a paragraph of a long letter to Hermite, the exponential of a matrix. In fact, the French mathematician confined himself to defining $e^X$ via the Taylor series expansion of $e^z$ at $0$ and remarking that the scalar identity $e^{x+y} = e^x e^y$ does not generalise to matrices. Surprisingly perhaps, truncating the Taylor expansion is still today, just over 150 years later, one of the most effective way of approximating the matrix exponential, along with the use of a family of rational approximants discovered a couple of decades later by Henri Padé, a student of Hermite’s.

Despite the fact that the exponential is an entire function, Taylor and Padé approximants converge only locally, and may require an approximant of exceedingly high order for matrices of large norm. Therefore, scaling the input is often necessary, in order to ensure fast and accurate computations. This can be effectively achieved by means of a technique called scaling and squaring, first proposed in 1967 by Lawson, who noticed that by exploiting the identity $e^A = \bigl(e^{2^{-s}A}\bigr)^{2^s}$ one can compute the matrix exponential by evaluating one of the approximants above at $2^{-s}A$ and then squaring the result $s$ times.

The choice of $s$ and of the order of the approximant hinges on a delicate trade-off that has been thoroughly studied in the last 50 years. The most efficient algorithms are tailored to double precision, and rely on the computation of precision-dependent constants, a relatively expensive task that needs to be performed only once during the design stage of the algorithm, and produces algorithms that are very efficient at runtime.

These techniques, however, do not readily extend to arbitrary precision environments, where the working precision at which the computation will be performed is known only at runtime. In the EPrint An arbitrary precision scaling and squaring algorithm for the matrix exponential, Nick Higham and I show how the scaling and squaring method can be adapted to the computation of the matrix exponential in precisions higher than double.

The algorithm we propose combines a new analysis of the forward truncation error of Padé approximants to the exponential with an improved strategy for selecting the parameters required to perform the scaling and squaring step. We evaluate at runtime a bound that arises naturally from our analysis in order to check whether a certain choice of those parameters will deliver an accurate evaluation of the approximant or not. Moreover, we discuss how this technique can be implemented efficiently for truncated Taylor series and diagonal Padé approximants.

Top left: forward error of several algorithms on a test set of non-Hermitian matrices. Top right: corresponding performance profile. Bottom: legend for the other two plots. Variants of our algorithm are in the top row, existing implementations for computing the matrix exponential in arbitrary precision are in the bottom row.

We evaluate experimentally several variants of the proposed algorithm in a wide range of precisions, and find that they all behave in a forward stable fashion. In particular, the most accurate of our implementations outperforms existing algorithms for computing the matrix exponential in high precision.

## Partnership with NAG on Numerical Linear Algebra for Emerging Computer Architectures

Professors Jack Dongarra, Nick Higham, and Françoise Tisseur are delighted to be working on a new three-year project with The Numerical Algorithms Group (NAG) to incorporate new numerical linear algebra routines into NAG software in order to exploit emerging computer architectures.

This Knowledge Transfer Partnership (KTP), funded by NAG and the Technology Strategy Board of Innovate UK, will

• develop new linear algebra kernels for multicore architectures,
• develop multiprecision codes that can exploit the half precision arithmetic increasingly available in hardware, and
• incorporate new sparse matrix functionality into the NAG Library.

The project builds on a long history of collaboration between NAG and the Manchester Numerical Linear Algebra Group.

In particular, they have successfully collaborated on two previous KTPs (2010–2013 and 2012–2013) as well as on joint MSc and PhD projects.