Numerical Python by Robert Johansson

Author:Robert Johansson
Language: eng
Format: epub
Publisher: Apress, Berkeley, CA

Eigenvalue Problems

Sparse eigenvalue and singular-value problems can be solved using the sp.linalg.eigs and sp.linalg.svds functions, respectively. For real symmetric or complex hermitian matrices, the eigenvalues (which in this case are real) and eigenvectors can also be computed using sp.linalg.eigsh. These functions do not compute all eigenvalues or singular values, but rather compute a given number of eigenvalues and vectors (the default is six). Using the keyword argument k with these functions, we can define how many eigenvalues and vectors should be computed. Using the which keyword argument, we can specify which k values are to be computed. The options for eigs are largest amplitude LM, smallest amplitude SM, largest real part LR, smallest real part SR, largest imaginary part LI, and smallest imaginary part SI. For svds only LM and SM are available.

For example, to compute the lowest four eigenvalues for the sparse matrix of the one-dimensional Poisson problem (of system size 10´10), we can use sp.linalg.eigs(A, k=4, which=’LM’):

In [54]: N = 10

In [55]: A = sp.diags([1, -2, 1], [1, 0, -1], shape=[N, N], format=’csc’)

In [56]: evals, evecs = sp.linalg.eigs(A, k=4, which=’LM’)

In [57]: evals

Out[57]: array([-3.91898595+0.j, -3.68250707+0.j, -3.30972147+0.j, -2.83083003+0.j])

The return value of sp.linalg.eigs (and sp.linalg.eigsh) is a tuple (evals, evecs) whose first element is an array of eigenvalues (evals), and the second element is an array (evecs) of shape , whose columns are the eigenvectors corresponding to the eigenvalues in evals. Thus, we expect that the dot product between A and a column in evecs is equal to the same column in evecs scaled by the corresponding eigenvalue in evals. We can directly confirm that this is indeed the case:

In [58]: np.allclose(A.dot(evecs[:,0]), evals[0] * evecs[:,0])

Out[58]: True

For this particular example, the sparse matrix A is symmetric, so instead of sp.linalg.eigs we could use sp.linalg.eigsh instead, and in doing so we obtain an eigenvalue array with real-valued elements:

In [59]: evals, evecs = sp.linalg.eigsh(A, k=4, which=’LM’)

In [60]: evals

Out[60]: array([-3.91898595, -3.68250707, -3.30972147, -2.83083003])

By changing the argument which=’LM’ (for largest magnitude) to which=’SM’ (smallest magnitude), we obtain a different set of eigenvalues and vector (those with smallest magnitude).

In [61]: evals, evecs = sp.linalg.eigs(A, k=4, which=’SM’)

In [62]: evals

Out[62]: array([-0.08101405+0.j, -0.31749293+0.j, -0.69027853+0.j, -1.16916997+0.j])

In [63]: np.real(evals).argsort()

Out[63]: array([3, 2, 1, 0])

Note that although we requested and obtained the four eigenvalues with smallest magnitude in the previous example, those eigenvalues and vectors are not necessarily sorted within each other (although they are in this particular case). Obtaining sorted eigenvalues is often desirable, and this is easily achieved with a small but convenient wrapper function that sorts the eigenvalues using NumPy’s argsort method. Here we give such a function, sp_eigs_sorted, which returns the eigenvalues and eigenvectors sorted by the real part of the eigenvalue.

In [64]: def sp_eigs_sorted(A, k=6, which=’SR’):

...: """ compute and return eigenvalues sorted by the real part """

...: evals, evecs = sp.linalg.eigs(A, k=k, which=which)

...: idx = np.real(evals).argsort()

...: return evals[idx], evecs[idx]

In [65]: evals, evecs = sp_eigs_sorted(A, k=4, which=’SM’)

In [66]: evals

Out[66]: array([-1.16916997+0.j, -0.69027853+0.j, -0.31749293+0.j, -0.08101405+0.j])

As a less trivial example using sp.linalg.eigs and the wrapper function sp_eigs_sorted, consider the spectrum of lowest eigenvalues of the linear combination of random sparse matrices M1 and M2.

Download

Copyright Disclaimer:
This site does not store any files on its server. We only index and link to content provided by other sites. Please contact the content providers to delete copyright contents if any and email us, we'll remove relevant links or contents immediately.