Determinar los valores y vectores propios de una matriz, aporta gran información acerca de las características y propiedades de esta, como también posee gran cantidad de aplicaciones prácticas como: Análisis de convergencia de sistemas dinámicos, PCA (Principal Component Analysis), análisis espectral, Eigenfaces, etc.
Sin embargo la determinación de los valores y vectores propios no es un problema simple. Como ya debe haber estudiado en cursos anteriores, existe un método directo basado en cálculo de las raíces del polinomio característico $p(x)$. Pero este problema resulta ser mal condicionado, esto es, a pequeñas variaciones en la matriz $A$ original, existe una gran variación en los resultados de los valores y vectores propios obtenidos (Ver polinomio de Wilkinson, texto guia).
En este notebook estudiaremos un método iterativo conocido como Power Iteration (y sus extensiones), que de modo similar a una iteración de punto fijo, permite obtener numéricamente los eigen(valores/vectores).
La motivación tras PI (Power Iteration) es que la multiplicación por matrices, tiende a "dirigir" a los vectores hacia el vector propio dominante (aquel con valor propio de mayor magnitud).
El algoritmo en cuestión es como sigue:
x = 'Initial guess'
for i in range n_iter:
u = x / ||x|| #normalization step
x = dot(A,u) #power iteration step
lamb = dot(u, dot(A, u)) #Rayleigh quotient
return x / ||x||
en donde se agrega una paso de normalización, para evitar que la magnitud del vector aumente sin límite, y el valor propio asociado se obtiene por medio del cociente de Rayleigh:
$$ \lambda = \frac{\mathbf{x}^*\,A\,\mathbf{x}}{\mathbf{x}^*\,\mathbf{x}} $$Para entender porque se de esta convergencia, considere una matriz $A \in \mathbb{R}^{m \times m}$ con valores propios reales $\lambda_1, \lambda_2, \ldots, \lambda_m$ tales que $|\lambda_1| > |\lambda_2| \geq |\lambda_3| \geq \ldots \geq |\lambda_m|$, tales que los vectores propios $\{\mathbf{v}_1, \mathbf{v}_2, \ldots, \mathbf{v}_m \}$ conforman una base de $\mathbb{R}^m$. Sea entonces $x_0$ el initial guess, este puede ser expresado como una combinación lineal de los vectores propios $\mathbf{v}_k$:
\begin{align} A x_0 &= c_1 \,A\,\mathbf{v}_1 + \cdots + c_m \,A\,\mathbf{v}_m = c_1 \lambda_1 \mathbf{v}_1 + \cdots + c_m \lambda_m \mathbf{v}_m \\ A^2 x_0 & = c_1 \lambda_1 \,A\,\mathbf{v}_1 + \cdots + c_m \lambda_m \,A\,\mathbf{v}_m = c_1 \lambda_1^2 \mathbf{v}_1 + \cdots + c_m \lambda_m^2 \mathbf{v}_m \\ \vdots &= \vdots \\ A^k x_0 &= c_1 \lambda_1^k \mathbf{v}_1 + \cdots + c_m \lambda_m^k \mathbf{v}_m \end{align}Factorizando $\lambda_1^k$ del último resultado se obtiene:
$$ \frac{A^k x_0}{\lambda_1^k} = c_1 \mathbf{v}_1 + c_2 \left(\frac{\lambda_2}{\lambda_1}\right)^k \mathbf{v}_2 + \cdots + c_m \left(\frac{\lambda_m}{\lambda_1}\right)^k \mathbf{v}_m$$Dado que $|\lambda_1|>|\lambda_i| \ \ \forall i \neq 1$, a medida que $k \rightarrow \infty$ todos los términos excepto el primero tienden a cero, con razón de convergencia $S \leq |\lambda_2/\lambda_1|$. Obteniendo como resultado un vector que es múltiplo del vector propio dominante.
Nota: Para más detalles revisar: Numerical Analysis, Tymothy Sauer, Chapter 12: Eigenvalues and Singular Values
In [1]:
import numpy as np
from scipy import linalg
from matplotlib import pyplot as plt
%matplotlib inline
In [2]:
A = np.array([[1, 0.5],[0.5, 1]])
x = np.array([1.,0.])
A = np.array([[1., 0.5,-0.1],[0.5, 1.,10.0],[2.,3.,5.]])
x = np.array([1.,0.,0.])
print("A =\n",A)
print("x =",x)
In [3]:
def power_iteration(A, x, k, verbose=False):
"""
Program 12.1 Power iteration
Computes dominant eigenvector of square matrix
Input: matrix A, initial (nonzero) vector x, number of steps k
Output: dominant eigenvalue lam, eigenvector u
"""
if verbose: print("Power Iteration Method\n%s"%('='*80))
for j in range(k):
u = x/np.linalg.norm(x)
x = np.dot(A, u)
lam = np.dot(u, x) #not really necessary to compute it at each iteration
if verbose: print("k=%d, lambda=%+.3f, u=%s"%(j,lam,str(u.T)))
u = x/np.linalg.norm(x)
if verbose: print("k=%d, lambda=%+.3f, u=%s\n"%(j+1,lam,str(u.T)))
return (lam, u)
In [4]:
# Testing algorithm
lam, u = power_iteration(A, x, 20, verbose=True)
print("lambda = {0}".format(lam))
print("u (dominant eigenvector) = {0}".format(u))
Una de las complicaciones que tiene el algoritmo anterior, es que sólo permite encontrar el valor y vectores propios dominantes. Luego ¿Cómo encontramos el resto?. Para responder esta pregunta, es necesario examinar dos propiedades importantes:
Tarea: Pruebe estas propiedades!
La idea es entonces realizar un shift $\widetilde{s}$ cercano a algún valor propio $s_k$, y computar PI sobre $(A - \widetilde{s}I)^{-1}$. Luego:
$$ |\lambda_k - \widetilde{s}| < |\lambda_i - \widetilde{s}| \leftrightarrow \bigg| \frac{1}{\lambda_k - \widetilde{s}} \bigg| > \bigg| \frac{1}{\lambda_i - \widetilde{s}} \bigg| \ \ \forall i \neq k \ $$entonces $\frac{1}{\lambda_k - \widetilde{s}}$ corresponderá con el vector propio dominante de $(A - \widetilde{s}\,I)^{-1}$. Notar que por lo enunciado en las propiedades, los vectores propios se mantienen sin alteraciones.
La idea anterior se ve reflejada en el algoritmo implementado a continuación:
In [5]:
def inverse_power_iteration(A, x, s, k, verbose=False):
"""
Program 12.2 Inverse Power iteration
Computes eigenvector of square matrix nearest to input s
Input: matrix A, initial (nonzero) vector x, shift s, number of steps k
Output: dominant eigenvalue lam, eigenvector of inv(A-sI)
"""
if verbose: print("Inverse Power Iteration Method\n%s"%('='*80))
As = A - s*np.eye(*A.shape)
for j in range(k):
u = x/np.linalg.norm(x)
x = np.linalg.solve(As, u) # Critical line!
lam = np.dot(u.T, x)
if verbose: print("k=%d, lambda=%+.3f, u=%s"%(j,1./lam+s,str(u.T)))
u = x/np.linalg.norm(x)
if verbose: print("k=%d, lambda=%+.3f, u=%s\n"%(j+1,1./lam+s,str(u.T)))
return (1./lam+s, u)
In [6]:
# Testing algoritm
lam, u = inverse_power_iteration(A, x, s=1./4, k=10, verbose=True)
print("lambda = {0}".format(lam))
print("v = {0}".format(u))
Como se analizó anteriormente, PI e Inverse PI tienen convergencia lineal con razón de convergencia $S \approx \frac{\lambda_2}{\lambda_1}$. Además sabemos que Inverse PI converge hacia el valor propio más cercano al shift, y que mientras más cercano sea el shift a tal valor, más rápido se logra la convergencia.
Entonces la idea de RQI es la siguiente: Si en cada iteración se tiene una aproximación del valor propio que andamos buscando, podemos ocupar esta aproximación como shift $s$, y dado que el shift será más cercano al valor propio, se acelerará la convergencia.
Tal valor aproximado es obtenido con el cociente de Rayleigh, y entonces el shift es actualizado con este cociente en cada iteración. Como resultado se produce el siguiente trade-off:
A continuación se presenta una implementación del RQI:
In [7]:
def rqi(A, x, k, verbose=False):
"""
Program 12.3 Rayleigh Quotient Iteration
Input: matrix A, initial (nonzero) vector x, number of steps k
Output: eigenvalue lam, eigenvector of inv(A-sI)
"""
if verbose: print("Rayleigh Quotient Iteration\n%s"%('='*80))
for j in range(k):
u = x/np.linalg.norm(x)
lam = np.dot(u.T, np.dot(A, u))
try:
x = np.linalg.solve(A -lam*np.eye(*A.shape), u)
except numpy.linalg.LinAlgError:
break
if verbose: print("k=%d, lambda=%+.3f, u=%s"%(j,lam,str(u.T)))
u = x/np.linalg.norm(x)
lam = float(np.dot(u.T, np.dot(A, u)))
if verbose: print("k=%d, lambda=%+.3f, u=%s\n"%(j+1,lam,str(u.T)))
return (lam, u)
Preguntas:
try y except en las líneas 11 y 13? ¿Que significa que el sistema no pueda ser resuelto?
In [8]:
# Testing algorithm
lam, v = rqi(A, x, k=2)
print("lambda = {0}".format(lam))
print("v = {0}".format(v))
La librería scipy tiene implementados algoritmos que permite calcular los valores y vectores propios. Las opciones posibles son:
En la librería scipy.linalg: eigvals/eigvalsh/eigvals_banded, eig/eigh/eig_banded,
En la librería scipy.sparse.linalg: eigen, eigs, eigsh.
En general siempre conviene utilizar las funciones desde scipy y no de numpy. La librería numpy hace un excelente trabajo al permitir el uso de vectores de tipo numérico, pero contiene solo algunos algoritmos numéricos y no necesariamente los más rápidos.
A continuación se muestra como utilizar algunas de estas funciones.
In [9]:
# Full matrices
from scipy import linalg as LA
N = 3
Aux = np.random.rand(N,N)
A = Aux + Aux.T # symmetric, so we'll deal with real eigs.
print(LA.eigvals(A)) # Only the eigenvalues, A not necessarily symmetric
print("*"*80)
print(LA.eigvalsh(A)) # Only the eigenvalues, A symmetric
print("*"*80)
print(LA.eig(A)) # All the eigenvalues and eigenvectors, A not necessarily symmetric
print("*"*80)
print(LA.eigh(A)) # All the eigenvalues and eigenvectors, A symmetric (faster)
print("*"*80)
lambdas, V = LA.eigh(A) # All the eigenvalues and eigenvectors, A symmetric (faster)
l1 = lambdas[0]
v1 = V[:,0]
print(l1)
print(v1)
print(np.dot(A, v1))
print(l1*v1)
Las matrices simétricas tiene una propiedad muy interesante:
En base a lo anterior se propone el siguiente algoritmo para encontrar los primeros $k$ valores/vectores propios:
def kEigenFinder(A, k, p):
m = A.shape[0]
lamb = 0.
v = np.zeros((m,1))
Lamb = []
V = []
for i in range(k):
A -= lamb*np.dot(v,v.T)
lamb,v = power_iteration(A, p)
Lamb.append(lamb)
V.append(v)
return Lamb,V
Construya una matriz simétrica de $100 \times 100$ y ejecute el kEigenFinder sobre tal matriz. Una forma fácil de construir una matriz simétrica es la matriz de covarianza:
$$\Sigma_X = \frac{1}{n-1}X^T X$$
donde $X \in \mathbb{R}^{m \times n}$, con $m$ samples y $n$ features.
Concluya acerca de la utilidad del procedimiento propuesto.
ctorres@inf.utfsm.cl) y ayudantes: Alvaro Salinas y Martín Villanueva. DI UTFSM. Abril 2016.El presente notebook ha sido creado para el curso ILI286 - Computación Científica 2, del Departamento de Informática, Universidad Técnica Federico Santa María.
El material ha sido creado por Claudio Torres ctorres@inf.utfsm.cl y Sebastian Flores sebastian.flores@usm.cl, y es distribuido sin restricciones. En caso de encontrar un error, por favor no dude en contactarnos.
[Update 2015] Se ha actualizado los notebooks a Python3 e includio el "magic" "%matplotlib inline" antes de cargar matplotlib para que los gráficos se generen en el notebook.
[Update 2016] (Martín) Modificaciones mayores al formato original. Agregado contexto: Introducción, marco teórico, explicaciones y tareas. Modificaciones menores en los algoritmos. Agregada la sección de Problema de Aplicación.
[Update 2018] (C.Torres) Using np.linalg.
[Update 2019] (C.Torres) Using * and mathbf for vectors. Fixing issue with titles of sections.
In [ ]: