6.2: Autosolucionadores numéricos

Como se discutió anteriormente, la teoría de la imposibilidad de Abel nos dice que no existe una fórmula algebraica general para calcular los valores propios de una\(N\times N\) matriz, para\(N\ge 5\). En la práctica, sin embargo, existen métodos numéricos, llamados autosolvers, que pueden calcular valores propios (y vectores propios) incluso para matrices muy grandes, con cientos de hilas/columnas, o mayores. ¿Cómo podría ser esto?

La respuesta es que los autosolucionadores numéricos son aproximados, no exactos. Pero aunque sus resultados no son exactos, son muy precisos, pueden acercarse a los valores propios exactos dentro de los límites de precisión fundamentales de la aritmética de punto flotante. Como estamos limitados por esos límites de precisión de punto flotante de todos modos, ¡eso es lo suficientemente bueno!

6.2.1 Croquis del Método Eigensolver

No entraremos en detalles sobre cómo funcionan los autosolucionadores numéricos, ya que eso implicaría una digresión bastante larga. Para los interesados, el siguiente trabajo brinda una buena presentación pedagógica de la asignatura:

• Bruno Lang, Solucionadores Directos para Problemas Simétricos de Autovalor. Métodos y Algoritmos Modernos de Química Cuántica, Proceedings, Segunda Ed. (2000). Enlace de descarga en PDF

Aquí hay un boceto muy breve del método básico. Similar a la eliminación gaussiana, el algoritmo contiene dos fases, una fase inicial relativamente costo/lenta y una segunda fase relativamente rápida. La primera fase, que se llama Reducción de Householder, aplica un conjunto cuidadosamente elegido de transformaciones de similitud a la matriz de entrada\(\mathbf{A}_0\):

\[\mathbf{A}_0 \;\rightarrow \;\mathbf{A}_1 = \mathbf{X}_1^{-1} \mathbf{A}_0 \mathbf{X}_1 \;\rightarrow \; \mathbf{A}_2 = \mathbf{X}_2^{-1} \mathbf{A}_1 \mathbf{X}_2, \;\; \mathrm{etc.}\]

El resultado final es una matriz\(\mathbf{A}_k\) que está en forma de Hessenberg: los elementos debajo de la primera subdiagonal son todos cero (los elementos inmediatamente debajo de la diagonal principal, es decir, a lo largo de la primera subdiagonal, pueden ser distintos de cero). Toda la fase de reducción de Householder requiere operaciones\(O(N^{3})\) aritméticas, donde\(N\) está el tamaño de la matriz.

La segunda fase del algoritmo se llama iteración QR. Usando un tipo diferente de transformación de similitud, los elementos a lo largo de la subdiagonal de la matriz de Hessenberg se reducen en magnitud. Cuando estos elementos se vuelven insignificantes, la matriz se vuelve triangular superior; en ese caso, los valores propios son simplemente los elementos a lo largo de la diagonal.

El proceso QR es iterativo, ya que reduce progresivamente la magnitud de los elementos de la matriz a lo largo de la subdiagonal. Formalmente, se requeriría un número infinito de iteraciones para reducir estos elementos a cero, ¡por eso no se viola el teorema de la imposibilidad de Abel! En la práctica, sin embargo, la iteración QR converge de manera extremadamente rápida, por lo que esta fase acaba tomando solo\(O(N^{2})\) tiempo.

De ahí que el tiempo de ejecución general para encontrar los valores propios de una matriz se escala como\(O(N^{3})\). Los vectores propios también se pueden calcular como un efecto secundario del algoritmo, sin un aumento adicional en el escalado de tiempo de ejecución.

6.2.2 Implementación de Python

Hay cuatro autosolvers numéricos principales implementados en Scipy, que se encuentran todos en el paquete scipy.linalg:

• scipy.linalg.eig devuelve los valores propios y vectores propios de una matriz.
• scipy.linalg.eigvals devuelve los valores propios (only) de una matriz.
• scipy.linalg.eigh devuelve los valores propios y vectores propios de una matriz hermitiana.
• scipy.linalg.eigvalsh devuelve los valores propios (only) de una matriz hermitiana.

La razón para tener cuatro funciones separadas es la eficiencia. Los tiempos de ejecución de las cuatro funciones escalan como\(O(N^{3})\), pero para cada uno\(N\) los tiempos de ejecución reales de eigvalsh y eigvalsh serán más cortos que eig y eigh, porque solo se le pide al solucionador propio que encuentre los valores propios y no es necesario construir los vectores propios. Además, eigvalsh es más rápido que eigvals, y eigh es más rápido que eig, porque el algoritmo eigensolver puede hacer uso de ciertos atajos numéricos que son válidos solo para matrices hermitianas.

Si pasas eigvalsh o eigh una matriz que en realidad no es hermitiana, los resultados son impredecibles; la función puede devolver un valor incorrecto sin señalar ningún error. Por lo tanto, solo debes usar estas funciones si estás seguro de que la matriz de entrada es definitivamente hermitiana (lo cual suele ser porque construiste la matriz de esa manera); si el mtrix es hermitiano, eigvalsh o eigh son ciertamente preferibles de usar, porque corren más rápido que sus Contrapartes no hermitianas.

Aquí hay un programa corto que usa eigvals para encontrar los valores propios de una\(3\times 3\) matriz:

from scipy import *
import scipy.linalg as lin

A = array([[1,3,1],[1, 3, 4],[2, 4, 2]])
lambd = lin.eigvals(A)

print(lambd)

Ejecución de los resultados del programa:

[ 7.45031849+0.j         -0.72515925+0.52865751j -0.72515925-0.52865751j]

El valor de retorno de eigvals es una matriz 1D de números complejos, que almacena los valores propios de la entrada. La función eigvalsh se comporta de manera similar, excepto que se devuelve una matriz real (ya que las matrices hermitianas tienen valores propios reales). Nota: no podemos usar lambda como nombre de una variable, porque lambda está reservada como una palabra clave especial en Python.

Aquí hay un ejemplo de uso de eig:

>>> A = array([[1,3,1],[1, 3, 4],[2, 4, 2]])
>>> lambd, Q = lin.eig(A)
>>> lambd
array([ 7.45031849+0.j    , -0.72515925+0.52865751j,   -0.72515925-0.52865751j])
>>> Q
array([[ 0.40501343+0.j        ,  0.73795979+0.j        ,  0.73795979-0.j        ],
       [ 0.65985810+0.j        , -0.51208724+0.22130102j, -0.51208724-0.22130102j],
       [ 0.63289132+0.j        ,  0.26316357-0.27377508j,  0.26316357+0.27377508j]])

La función eig devuelve un par de valores; la primera es una matriz 1D de valores propios (que nombramos lambd en el ejemplo anterior), y la segunda es una matriz 2D que contiene los vectores propios correspondientes en cada columna (que nombramos Q). Por ejemplo, se puede acceder al primer vector propio con Q [:,0]. Podemos verificar que esto es de hecho un vector propio:

>>> dot(A, Q[:,0])
array([ 3.01747903+0.j,  4.91615298+0.j,  4.71524187+0.j])
>>> lambd[0] * Q[:,0]
array([ 3.01747903+0.j,  4.91615298+0.j,  4.71524187+0.j])

La función eigh se comporta de manera similar, excepto que la matriz 1D de valores propios es real.

6.2.3 Problema Generalizado de Autovalor

A veces, también puede encontrarse con problemas generalizados de valores propios, que tienen la forma

\[\mathbf{A}\, \vec{x} = \lambda\, \mathbf{B}\, \vec{x},\]

para matrices cuadradas de igual tamaño conocidas\(\mathbf{A}\) y\(\mathbf{B}\). Llamamos\(\lambda\) un “autovalor generalizado”, y\(\vec{x}\) un “autovector generalizado”, para el par de matrices\((\mathbf{A},\mathbf{B})\). El problema del valor propio generalizado se reduce al problema del valor propio ordinario cuando\(\mathbf{B}\) es la matriz de identidad.

La manera ingenua de resolver el problema del valor propio generalizado sería calcular la inversa de\(\mathbf{B}^{-1}\), y luego resolver el problema del valor propio para\(\mathbf{B}^{-1}\mathbf{A}\). Sin embargo, resulta que el problema del valor propio generalizado puede resolverse directamente con solo ligeras modificaciones al algoritmo habitual de autosolucionador numérico. De hecho, los autosolvers de Scipy descritos en la sección anterior resolverán el problema de autovalores generalizados si pasa una matriz 2D\(\mathbf{B}\) como segunda entrada (si se omite esa segunda entrada, los autosolvers resuelven el problema de autovalores ordinarios, como se describió anteriormente).

Aquí hay un programa de ejemplo para resolver un problema de autovalor generalizado:

from scipy import *
import scipy.linalg as lin

A = array([[1,3,1],[1, 3, 4],[2, 4, 2]])
B = array([[0,2,1], [0, 1, 1], [2, 0, 1]])

lambd, Q = lin.eig(A, B)

## Verify the solution for first generalized eigenvector:
lhs = dot(A,Q[:,0])             # A . x
rhs = lambd[0] * dot(B, Q[:,0]) # lambda B . x

print(lhs)
print(rhs)

Ejecutando las impresiones del programa anterior:

[-0.16078694+0.j -0.07726949+0.j  0.42268561+0.j]
[-0.16078694+0.j -0.07726949+0.j  0.42268561+0.j]

Los autosolvers hermitianos, eigh y eigvalsh, se pueden utilizar para resolver el problema del autovalor generalizado solo si tanto las matrices como\(\mathbf{A}\) las\(\mathbf{B}\) matrices son hermitianas.