We consider the following four cases:

The following functions in the modules scipy.linalg and scipy.sparse.linalg help us to compute eigenvalues and eigenvectors:

For any kind of eigenvalue problem where the matrices are not symmetric or banded, we use the function eig, which is a wrapper for the LAPACK routines GEEV and GGEV (the latter for generalized eigenvalue problems). The function eigvals is syntactic sugar for a case of eig that only outputs the eigenvalues, but not the eigenvectors. To report whether we require left of right eigenvectors, we use the optional Boolean parameters left and right. By default, left is set to False and right to True, hence offering right eigenvectors.

For eigenvalue problems with non-banded real symmetric or Hermitian matrices, we use the function eigh, which is a wrapper for the LAPACK routines of the form *EVR, *GVD, and *GV. We are given the choice to output as many eigenvalues as we want, with the optional parameter eigvals. This is a tuple of integers that indicate the indices of the lowest and the highest eigenvalues required. If omitted, all eigenvalues are returned. In such a case, it is possible to perform the computation with a much faster algorithm based on divide and conquer techniques. We may indicate this choice with the optional Boolean parameter turbo (by default set to False).

If we wish to report only eigenvalues, we can set the optional parameter eigvals_only to True, or use the corresponding syntactic sugar eighvals.

The last case that we contemplate in the scipy.linalg module is that of the eigenvalue problem of a banded real symmetric or Hermitian matrix. We use the function eig_banded, making sure that the input matrices are in the AB format. This function is a wrapper for the LAPACK routines *EVX.

For extremely large matrices, the computation of eigenvalues is often computationally impossible. If these large matrices are sparse, it is possible to calculate a few eigenvalues with two iterative algorithms, namely the Implicitly Restarted Arnoldi and the Implicitly Restarted Lanczos methods (the latter for symmetric or Hermitian matrices). The module scipy.sparse.linalg has two functions, eigs and eigsh, which are wrappers to the ARPACK routines *EUPD that perform them. We also have the function lobpcg that performs another iterative algorithm, the Locally Optimal Block Preconditioned Conjugate Gradient method. This function accepts a Preconditioner, and thus has the potential to converge more rapidly to the desired eigenvalues.

We will illustrate the usage of all these functions with an interesting matrix: Andrews. It was created in 2003 precisely to benchmark memory-efficient algorithms for eigenvalue problems. It is a real symmetric sparse matrix with size 60,000 × 60,000 and 760,154 non-zero entries. It can be downloaded from the Sparse Matrix Collection at www.cise.ufl.edu/research/sparse/matrices/Andrews/Andrews.html.

For this example, we downloaded the matrix in the Matrix Market format Andrews.mtx. Note that the matrix is symmetric, and the file only provides data on or below the main diagonal. After collecting all this information, we ensure that we populate the upper triangle too:

In [1]: import numpy as np, scipy.sparse as spsp, \
   ...: scipy.sparse.linalg as spspla
In [2]: np.set_printoptions(suppress=True, precision=6)
In [3]: rows, cols, data = np.loadtxt("Andrews.mtx", skiprows=14,
   ...:                                unpack=True); \
   ...: rows-=1; \
   ...: cols-=1
In [4]: A = spsp.csc_matrix((data, (rows, cols)), \
   ...:                     shape=(60000,60000)); \
   ...: A = A + spsp.tril(A, k=1).transpose()

We compute first the top largest five eigenvalues in absolute value. We call the function eigsh, with the option which='LM'.

In [5]: %time eigvals, v = spspla.eigsh(A, 5, which='LM')
CPU times: user 3.59 s, sys: 104 ms, total: 3.69 s
Wall time: 3.13 s
In [6]: print eigvals
[ 69.202683  69.645958  70.801108  70.815224  70.830983]

We may compute the smallest eigenvalues in terms of the absolute value too, by switching to the option which='SM':

In [7]: %time eigvals, v = spspla.eigsh(A, 5, which='SM')
CPU times: user 19.3 s, sys: 532 ms, total: 19.8 s
Wall time: 16.7 s
In [8]: print eigvals
[ 10.565523  10.663114  10.725135  10.752737  10.774503]

In [9]: A = spspla.aslinearoperator(A)
In [10]: %time spspla.eigsh(A, 5, sigma=10.0, mode='cayley')
CPU times: user 2min 5s, sys: 916 ms, total: 2min 6s
Wall time: 2min 6s
In [11]: print eigvals
[ 10.565523  10.663114  10.725135  10.752737  10.774503]

There are four cases:

There are four functions in the module scipy.linalg that provide us with tools to compute any of these decompositions:

The function hessenberg gives us the first step in the computation of any Schur decomposition. This is a factorization of any square matrix A in the form A = Q ● U ● Q^H, where Q is unitary and U is an upper Hessenberg matrix (all entries are zero below the sub-diagonal). The algorithm is based on the combination of the LAPACK routines GEHRD, GEBAL (to compute U), and the BLAS routines GER, GEMM (to compute Q).

The functions schur and qz are wrappers to the LAPACK routines GEES and GGES, to compute the normal and generalized Schur decompositions (respectively) of square matrices. We choose whether to report complex or real decompositions on the basis of the optional parameter output (which we set to 'real' or 'complex'). We also have the possibility of sorting the eigenvalues in the matrix representation. We do so with the optional parameter sort, with the following possibilities: