8.3: Uso de matrices dispersas

Se deben usar formatos de matriz dispersa, en lugar de matrices 2D convencionales, cuando se trata de matrices dispersas de tamaño\(10^3\times 10^3\) o más grandes. (Para matrices de tamaño\(100\times 100\) o menos, las diferencias en el rendimiento y el uso de memoria suelen ser insignificantes).

Por lo general, es bueno elegir el formato CSR o CSC, dependiendo de qué operaciones matemáticas pretendas realizar. Se puede construir la matriz ya sea (i) directamente, por medio de una tupla (data, idx) como se describió anteriormente, o (ii) creando una matriz LIL, rellenando los elementos deseados y luego convirtiéndola a formato CSR o CSC.

Si se trata de una matriz dispersa que es “fuertemente diagonal” (es decir, los elementos distintos de cero ocupan un número muy pequeño de diagonales, como una matriz tridiagonal), entonces puede considerar usar el formato DIA. La principal ventaja del formato DIA es que es muy fácil de construir, al suministrar una entrada (data, offsets) a la función dia_matrix, como se describió anteriormente. Sin embargo, el formato generalmente no se desempeña significativamente mejor que CSR/CSC; y si la matriz no es fuertemente diagonal, su rendimiento es mucho peor.

Otra forma común de construir una matriz dispersa es usar las funciones scipy.sparse.diags o scipy.sparse.spdiags. Estas funciones permiten especificar el contenido de la matriz en términos de sus diagonales, así como qué formato disperso usar. Las dos funciones tienen convenciones de llamada ligeramente diferentes; consulte la documentación para más detalles.

8.3.1 El método de punto

Cada matriz dispersa tiene un método de punto, que calcula el producto de la matriz con una entrada (en forma de matriz 1D o 2D, o matriz dispersa), y devuelve el resultado. Para los productos de matriz dispersa, este método debe usarse en lugar de la función independiente, llamada de manera similar punto, que se usa para productos de matriz no dispersa. El método de matriz dispersa hace uso de la dispersión de matriz para acelerar el cálculo del producto. Por lo general, el formato CSR es más rápido en esto que los otros formatos escasos.

Por ejemplo, supongamos que queremos calcular el producto\(\mathbf{A} \vec{x}\), donde

\[\mathbf{A} = \begin{bmatrix}0 & 1.0 & 0 & 0 & 0 \\ 0 & 2.0 & -1.0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 6.6 & 0 & 0 & 0 & 1.4\end{bmatrix}, \quad \vec{x} = \begin{bmatrix}1\\1\\2\\3 \\5\end{bmatrix}.\]

Esto podría lograrse con el siguiente programa:

from scipy import *
import scipy.sparse as sp

data = [1.0, 2.0, -1.0, 6.6, 1.4]
rows = [0, 1, 1, 3, 3]
cols = [1, 1, 2, 0, 4]
A = sp.csr_matrix((data, [rows, cols]), shape=(4,5))
x = array([1.,1.,2.,3.,5.])

y = A.dot(x)

print(y)

Ejecutar este programa da el resultado esperado:

[  1.    0.    0.   13.6]

8.3.2 spsolve

La función spsolve, proporcionada en el módulo scipy.sparse.linalg, es un solucionador para un sistema lineal escaso de ecuaciones. Se necesitan dos entradas,\(\mathbf{A}\) y\(\mathbf{b}\), donde\(\mathbf{A}\) debe ser una matriz dispersa;\(\mathbf{b}\) puede ser dispersa, o una matriz 1D o 2D ordinaria. Devuelve el\(\mathbf{x}\) que resuelve las ecuaciones lineales\(\mathbf{x}\). El valor de retorno puede ser una matriz convencional o una matriz dispersa, dependiendo de si\(\mathbf{b}\) es escasa.

Para problemas escasos, siempre debes usar spsolve en lugar de scipy.linalg.solve (el solucionador habitual para problemas no dispersos). Aquí hay un programa de ejemplo que muestra spsolve en acción:

from scipy import *
import scipy.sparse as sp
import scipy.sparse.linalg as spl

## Make a sparse matrix A and a vector x
data = [1.0, 1.0, 1.0, 1.0, 1.0]
rows = [0, 1, 1, 3, 2]
cols = [1, 1, 2, 0, 3]
A = sp.csr_matrix((data, [rows, cols]), shape=(4,4))
b = array([1.0, 5.0, 3.0, 4.0])

## Solve Ax = b
x = spl.spsolve(A, b)
print(" x = ", x)

## Verify the solution:
print("Ax = ", A.dot(x))
print(" b = ", b)

Ejecutar el programa da:

 x =  [ 4.  1.  4.  3.]
Ax =  [ 1.  5.  3.  4.]
 b =  [ 1.  5.  3.  4.]

8.3.3 eigs

Para problemas de valores propios que involucran matrices dispersas, normalmente uno no intenta encontrar todos los valores propios (y vectores propios). Las matrices dispersas suelen ser tan grandes que resolver el problema del valor propio completo tomaría un tiempo imprácticamente largo, incluso si recibimos una aceleración de la aritmética de matriz dispersa. Por suerte, en la mayoría de las situaciones solo necesitamos encontrar un subconjunto de los valores propios (y vectores propios). Por ejemplo, después de discretizar la ecuación de onda de Schrödinger 1D, normalmente solo nos interesan los varios valores propios de energía más bajos.

La función eigs, proporcionada en el módulo scipy.sparse.linalg, es un solucionador propio para matrices dispersas. A diferencia de los autosolvers que hemos discutido anteriormente, como scipy.linalg.eig, la función eigs solo devuelve un subconjunto especificado de los valores propios y vectores propios.

La función eigsh es similar a eigs, excepto que está especializada para matrices hermitianas. Ambas funciones hacen uso de una biblioteca de solucionador propio numérico de bajo nivel llamada ARPACK, que también es utilizada por GNU Octave, Matlab y muchas otras herramientas numéricas. No discutiremos cómo funciona el algoritmo.

La primera entrada a eigs o eigsh es la matriz para la que queremos encontrar los valores propios. También se aceptan varias otras entradas opcionales. Aquí están los más utilizados:

• El parámetro opcional llamado k especifica el número de valores propios a encontrar (el valor predeterminado es 6).
• El parámetro opcional llamado M, si se suministra, especifica una matriz derecha para un problema de autovalor generalizado.
• El parámetro opcional llamado sigma, si se suministra, debe ser un número; significa encontrar los k valores propios que están más cerca de ese número.
• El parámetro opcional llamado que especifica qué valores propios encontrar, usando un criterio diferente de sigma: 'LM' significa encontrar los valores propios con las magnitudes más grandes, 'SM' significa encontrar aquellos con las magnitudes más pequeñas, etc. especificar simultáneamente ambos sigma y cuál. Al encontrar valores propios pequeños, suele ser mejor usar sigma en lugar de cuál (ver la discusión en la siguiente sección).
• El parámetro opcional llamado return_eigenvectors, si es True (el valor predeterminado), significa devolver tanto los valores propios como los vectores propios. Si es False, la función devuelve solo los valores propios.

Para obtener la lista completa de entradas, consulte la documentación completa de funciones para eigs y eigsh.

Por defecto, eigs y eigsh devuelven dos valores: una matriz 1D (que es compleja para eigs y real para eigsh) que contiene los valores propios encontrados, y una matriz 2D donde cada columna es uno de los vectores propios correspondientes. Si el parámetro opcional llamado return_eigenvectors se establece en False, entonces solo se devuelve la matriz 1D de valores propios.

En la siguiente sección, veremos un ejemplo de uso de eigsh.