FEAST: A New Eigensolver for the nAG Library

The nAG Library contains a considerable array of state-of-the-art routines for solving eigenvalue problems. From the dense QR and MR³ algorithms of LAPACK (see Chapter F08) to the ARPACK solvers for sparse eigenproblems (Chapter F12), we have things well covered. It might be a surprise, then, to hear about a new suite of eigensolvers, introduced at Mark 27.3.

The FEAST Eigenvalue Solver was developed by Professor Eric Polizzi of the University of Massachusetts, along with numerous collaborators. It works in a completely different way from other eigensolvers (more on this shortly) and, as a result, has some very useful properties.

• It searches for only those eigenvalues lying within a user-specified contour in the complex plane.
• It can be used to solve standard, generalized and polynomial eigenvalue problems.
• Solvers are available for complex, real, symmetric, Hermitian, and non-symmetric problems.
• It can be used for both large, sparse and small, dense problems.

It is the first property above that really differentiates FEAST from other eigensolvers. This is particularly useful when applied to large, sparse eigenproblems when a subset of the spectrum is required. Solvers such as ARPACK are typically only able to find the largest/smallest magnitude eigenvalues or a specified number of eigenvalues close to a given complex value. Only the FEAST eigensolver can search within a specific contour within which the number of eigenvalues might be unknown. We are aware of several areas in which such functionality is useful, including modelling the oscillations of railway tracks, simulation of fluid flows and electronic structure computations. (If you discover another interesting application, please let us know!)

FEAST works by exploiting the rich mathematics of calculus in the complex plane. In particular, it relies on Cauchy’s residue theorem, which states that the integral of a function, \(f(z)\), around a closed contour can be found by considering poles of \(f\) lying within the contour. The trick, then, is to choose \(f\) such that its poles correspond to the required eigenvalues. The FEAST algorithm achieves this by evaluating:

where \(\gamma\) is a closed contour, \(A\) is the \(n\times n\) matrix whose eigenvalues are sought, and \(Y\) is an initial guess at the eigenvector subspace. The number of columns of \(Y\) is typically taken to be far less than \(n\). This means that the matrix \(Q\) can then be used by FEAST to form a reduced-size eigenvalue problem for the eigenvalues within the contour.

We will demonstrate FEAST using the Python interfaces to the nAG Library. The FEAST routines themselves are found in the sparseig submodule. Let’s use the matrix qc324 from the Matrix Market repository. First, we will import a few modules and read in the matrix.

import numpy as np
import scipy.io
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import math
from naginterfaces.library import machine, sparse, sparseig

A = scipy.io.mmread('qc324.mtx.gz')

Figure 1: Gershgorin disks for qc324

In many applications, something is already known about the eigenspace and where to seek eigenvalues. But an alternative approach, which we will illustrate here, is to use Gershgorin discs to guide us. Each Gershgorin disc is centred on one of the diagonal entries of \(A\) and has a radius given by the sum of the absolute values of the off-diagonals in the corresponding row. It can be shown that each eigenvalue of \(A\) lies within at least one of the Gershgorin discs. Figure 1 shows the Gershgorin disks for qc324. It was generated using the following code snippet.

D = A.diagonal(0)
n = A.shape[0]

fig, ax = plt.subplots()
for i in range(n):
    x, y = D[i].real, D[i].imag
    ax.add_artist( plt.Circle((x, y), sum(abs(A.getrow(i).data)) - abs(D[i]), alpha=0.5))

Figure 2: Contour choice for qc324

Let’s use FEAST to search for eigenvalues in the contour shown in Figure 2, which is an ellipse centred at \(1.4-0.1i\), with semi-major and semi-minor axes 0.4 and 0.24 respectively (note that FEAST also allows user-defined contours made up of combinations of lines and arcs, and we have examples of these bundled with the Library). Since the Gershgorin discs are intersecting, it is possible that there are no eigenvalues within this contour — FEAST will be able to tell us if this is the case.

There is some setup to do before calling the solver. We need to initialize a data handle, which stores data for the suite. Then we set an option determining the eccentricity of the ellipse. Finally, we call a contour setting routine.

# Initialize the FEAST data handle
handle = sparseig.feast_init()

# Use option setting routine to specify axis ratio for ellipse
sparseig.feast_option(handle, 'Ellipse Contour Ratio = 60')

# Generate contour
sparseig.feast_gen_contour(handle, complex(1.4, -0.1), 0.4)

We are almost ready to call a solver to find the eigenvalues. There are various solvers available, but we will use feast_complex_symm_solve, which is designed for complex symmetric matrices. The solvers use reverse communication, returning to the calling program to solve linear systems and perform matrix multiplications. So first we will declare some functions to perform these operations.

def form_matrix(ze, A):
    """ Form the sparse matrix ze - A """
    n = A.shape[0]
    az_tmp = np.concatenate((-A.data,ze*np.ones(n, dtype=np.complex128)))
    irowz_tmp = np.concatenate((A.row+1,np.arange(1,n+1)))
    icolz_tmp = np.concatenate((A.col+1,np.arange(1,n+1)))
    # Sort elements into correct coordinate storage format
    az, irowz, icolz, _ = sparse.complex_gen_sort(n, az_tmp, irowz_tmp, icolz_tmp, 'S', 'R')
    return az, irowz, icolz

def incomplete_lu(az, irowz, icolz, n):
    """ Form the incomplete LU factorization of a sparse matrix """
    nnaz = np.size(az)
    az = np.concatenate((az, np.zeros(2*nnaz)))
    irowz = np.concatenate((irowz, np.zeros(2*nnaz, dtype=int)))
    icolz = np.concatenate((icolz, np.zeros(2*nnaz, dtype=int)))
    az, irowz, icolz, ipivp, ipivq, istr, idiag, _, _
		= sparse.complex_gen_precon_ilu(nnaz, az, irowz, icolz, n, 0, 0.0, 'P', 'N')
    return az, irowz, icolz, ipivp, ipivq, istr, idiag

def solve_lin_sys(
        nnz_az, az, irowz, icolz, ipivp, ipivq, istr, idiag, y):
    """ Solve a sparse linear system """
    n = np.size(y, 0)
    m0 = np.size(y, 1)
    w = np.zeros(n, dtype=np.complex128)
    for i in range(m0):
        w[0:n] = y[0:n, i]
        y[0:n, i] = 1.0
        y[0:n, i], _, _ = sparse.complex_gen_solve_ilu(
            'RGMRES', nnz_az, az, irowz, icolz, ipivp,
            ipivq, istr, idiag, w, 30,
            math.sqrt(machine.precision()), 500, y[0:n, i])

def sparse_matmul(A, x, z):
    """ Perform sparse matrix multiplication """
    n = np.size(x, 0)
    m0 = np.size(x, 1)
    for i in range(m0):
        x[0:n, i] = sparse.complex_gen_matvec(
            'N', A.data, A.row + 1,
            A.col + 1, z[0:n,i], 'N')

The code calling the reverse communication solver is shown below. Depending on the value of irevcm, we need to either: form \(ze – A\) (where \(ze\) is a complex scalar returned by the solver) and factorize it; use the factorization to solve a linear system; or perform a matrix multiplication.

irevcm = 0
ze = 0.0 + 0.0j
m0 = 30
eps = 0.0
itera = 0
x = np.zeros((n, m0), dtype=np.complex128)
y = np.zeros((n, m0), dtype=np.complex128)
z = np.zeros((n, m0), dtype=np.complex128)
d = np.zeros(m0, dtype=np.complex128)
resid = np.zeros(m0, dtype=np.float64)

# Call the reverse communication solver
while True:
    irevcm, ze, m0, nconv, eps, itera = \
        sparseig.feast_complex_symm_solve(
            handle, irevcm, ze, x,
            y, m0, d, z, eps, itera, resid)

    if irevcm == 0:
        break

    if irevcm == 1:
        # Form the sparse matrix ze - A
        az, irowz, icolz = form_matrix(ze, A)
        nnz_az = np.size(az)
        # Form incomplete LU factorization of ze - A
        az, irowz, icolz, ipivp, ipivq, istr, idiag = \
            incomplete_lu(az, irowz, icolz, n)

    if irevcm == 2:
        # Solve (ze - A) y = w, with m0 righthand sides
        solve_lin_sys(nnz_az, az, irowz, icolz, ipivp, ipivq,
                      istr, idiag, y)

    if irevcm == 3:
        sparse_matmul(A, x, z) # Compute x <- A z

    if irevcm == 4:
	# If we were solving a generalized eigenvalue problem
	# then here we would compute x <- Bz
        x[:,:] = z[:,:] # Compute x <- z

The eigenvalues that FEAST found within the contour are shown in Figure 3.

Figure 3: Eigenvalues of qc324