# The Discretized Laplace Operator on Hyperrectangles with Zero Dirichlet Boundary Conditions

In this blog post, I present stiffness and mass matrix as well as eigenvalues and eigenvectors of the Laplace operator (Laplacian) on domains $(0, \ell)$, $(0, \ell_1) \times (0, \ell_2)$, and so on (hyperrectangles) with zero Dirichlet boundary conditions discretized with the finite difference method (FDM) and the finite element method (FEM) on equidistant grids. For the FDM discretization, we use the central differences scheme with the standard five-point stencil in 2D. For the FEM, the ansatz functions are the hat functions. The matrices, standard eigenvalue problems $A v = \sigma v$, and generalized eigenvalue problems $K w = \tau M w$ arising from the discretization lend themselves for test problems in numerical linear algebra because they are well-conditioned, not diagonal, and the matrix dimension can be increased arbitrarily.

Python code generating the matrices and their eigenpairs can be found in my git repository discrete-laplacian.

# SuperLU vs Direct Substructuring

The eigenproblem solver in my master's thesis used SuperLU, a direct solver for the solution of systems of linear equations (SLE) $Ax = b$. For the largest test problems, the eigensolver ran out of memory when decomposing the matrix $A$ which is why I replaced SuperLU with direct substructuring in an attempt to reduce memory consumption. For this blog post, I measured set-up time, solve time, and memory consumption of SuperLU and direct substructuring with real symmetric positive definite real-world matrices for SLEs with a variable number of right-hand sides, I will highlight that SuperLU was deployed with a suboptimal parameter choice, and why the memory consumption of the decomposition of $A$ is the wrong objective function when you want to avoid running out of memory.

# Master's Thesis: Projection Methods for Generalized Eigenvalue Problems

My master's thesis deals with dense and sparse solvers for generalized eigenvalue problems (GEPs) with Hermitian positive semidefinite matrices. Key results are

• structure-preserving backward error bounds computable in linear time,
• the runtime of GSVD-based dense GEP solvers is within factor 5 of the fastest GEP solver with Netlib LAPACK in my tests,
• computing the GSVD directly is up to 20 times slower than the computation by means of QR factorizations and the CS decomposition with Netlib LAPACK in my tests,
• given a pair of matrices with 2x2 block structure, I show how to minimize eigenvalue perturbation by off-diagonal blocks with the aid of graph algorithms, and
• I propose a new multilevel eigensolver for sparse GEPs that is able to compute up to 1000 eigenpairs on a cluster node with two dual-core CPUs and 16 GB virtual memory limit for problems with up to 150,000 degrees of freedom in less than eleven hours.

The revised edition of the thesis with fixed typos is here (PDF), the source code is available here, and the abstract is below. In February, I already gave a talk on the preliminary thesis results; more details can be found in the corresponding blog post.

###### Abstract

This thesis treats the numerical solution of generalized eigenvalue problems (GEPs) $Kx = \lambda Mx$, where $K$, $M$ are Hermitian positive semidefinite (HPSD). We discuss problem and solution properties, accuracy assessment of solutions, aspect of computations in finite precision, the connection to the finite element method (FEM), dense solvers, and projection methods for these GEPs. All results are directly applicable to real-world problems.

We present properties and origins of GEPs with HPSD matrices and briefly mention the FEM as a source of such problems.

With respect to accuracy assessment of solutions, we address quickly computable and structure-preserving backward error bounds and their corresponding condition numbers for GEPs with HPSD matrices. There is an abundance of literature on backward error measures possessing one of these features; the backward error in this thesis provides both.

In Chapter 3, we elaborate on dense solvers for GEPs with HPSD matrices. The standard solver reduces the GEP to a standard eigenvalue problem; it is fast but requires positive definite mass matrices and is only conditionally backward stable. The QZ algorithm for general GEPs is backward stable but it is also much slower and does not preserve any problem properties. We present two new backward stable and structure preserving solvers, one using deflation of infinite eigenvalues, the other one using the generalized singular value decomposition (GSVD). We analyze backward stability and computational complexity. In comparison to the QZ algorithm, both solvers are competitive with the standard solver in our tests. Finally, we propose a new solver combining the speed of deflation with the ability of GSVD-based solvers to handle singular matrix pencils.

Finally, we consider black-box solvers based on projection methods to compute the eigenpairs with the smallest eigenvalues of large, sparse GEPs with Hermitian positive definite matrices (HPD). After reviewing common methods for spectral approximation, we briefly mention ways to improve numerical stability. We discuss the automated multilevel substructuring method (AMLS) before analyzing the impact of off-diagonal blocks in block matrices on eigenvalues. We use the results of this thesis and insights in recent papers to propose a new divide-and-conquer eigensolver and to suggest a change that makes AMLS more robust. We test the divide-and-conquer eigensolver on sparse structural engineering matrices with 10,000 to 150,000 degrees of freedom.

2010 Mathematics Subject Classification. 65F15, 65F50, 65Y04, 65Y20.

Edit: Revised master's thesis from April 2016 (PDF)

# Advice: Implementing a Solver using LAPACK

For my master's thesis I implemented multiple solvers for structured generalized eigenvalue problems using LAPACK. In this post, I will briefly discuss a method to simplify the memory management and ways to catch programming errors as early as possible when implementing a solver for linear algebra problems that uses LAPACK. The advice in this post only supplements good programming practices like using version control systems and automated tests.

# Talk: Master's Thesis

Yesterday I gave a talk on the preliminary results of my master's thesis on Projection Methods for Generalized Eigenvalue Problems. Below is the announcement for the talk; it mentions only complex (Hermitian) matrices but all results are immediately applicable to real matrices. The slides are here (PDF).

###### Projection Methods for Generalized Eigenvalue Problems

In this talk, we will discuss backward errors, dense solvers, and a new projection method for generalized eigenvalue problems (GEPs) with Hermitian positive semidefinite (HPSD) matrices.

We will present properties and origins of such GEPs. We will also address quickly computable and structure preserving backward error bounds for these kinds of GEPs. There is an abundance of literature on backward error measures possessing one of these features but only recently, the author came across a backward error providing both.

We will elaborate on dense solvers for GEPs with HPSD matrices. The standard solver for GEPs with Hermitian matrices is fast but requires positive definite mass matrices and is only conditionally backward stable; the QZ algorithm for general GEPs is backward stable but it is also magnitudes slower and does not preserve any problem properties. In the talk, we will present two new backward stable and structure preserving solvers, one using deflation, the other one using the generalized singular value decomposition (GSVD). In comparison to the QZ algorithm, both solvers are competitive with the standard solver in our tests.

Finally, we will touch on a new solver for large, sparse GEPs.

# Evaluating Approximate Solutions using the Backward Error

Given a linear algebra problem and an approximate solution, can we say something about the quality of the given approximate solution without knowing an exact solution?

# Improved Error Bounds for "Accuracy and Stability of Numerical Algorithms"

Higham's Accuracy and Stability of Numerical Algorithms is regularly needed for my work. The second edition of the book was released in 2002 so some of the bounds do not reflect the state of the art in 2015. I collected some papers with improved bounds which are relevant for me.

# Spectral Norm Bounds for Hermitian Matrices

In the mathematical literature, numerous bounds can be found for the spectral norm of a matrix. Can we improve these bounds if the matrices are known to be Hermitian or even real symmetric semidefinite?

# Performance and Accuracy of xPTEQR

xPTEQR is a LAPACK function for the computation of the  eigendecomposition of a symmetric positive definite tridiagonal matrix. I compare the performance and accuracy of xPTEQR to the other symmetric eigensolvers in LAPACK 3.5.0.

In layman’s terms, there are several algorithms available on a computer that calculate the eigenvectors and eigenvalues of matrices with special properties (symmetric, positive definite, tridiagonal). xPTEQR is the name of the computer implementation of one of those algorithms and I am not aware of any measurements with it. In this blog post I compare xPTEQR in terms of speed and accuracy to other good algorithms. In the figures below, xPTEQR can be found under the label "SVD" (xPTEQR is the name of the implementation, SVD is the algorithm). Moreover, in every figure it holds that lower is better.