<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <center>
    <h1> ILI286 - Computación Científica II  </h1>
    <h2> Valores y Vectores Propios </h2>
    <h2> <a href="#acknowledgements"> [S]cientific [C]omputing [T]eam </a> </h2>
    <h2> Version
Step1:  Matriz y vector de prueba
Step3:  <div id='pi' />
Power Iteration
A continuación se entrega el código del algoritmo clásico de Power Iteration. Pruebe cambiando las matrices y los parámetros del algoritmo.
Step5:  <div id='invpi' />
Inverse Power Iteration
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
Step7:  <div id='rq' />
Rayleigh Quotient Iteration
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
Step8:  Preguntas
Step9:  <div id='sp' />
$\texttt{SciPy}$ Eigenvalue
La librería scipy tiene implementados algoritmos que permite calcular los valores y vectores propios. Las opciones posibles son
<ASSISTANT_TASK:>
Python Code:
import numpy as np
from scipy import linalg
from matplotlib import pyplot as plt
%matplotlib inline
Explanation: <center>
    <h1> ILI286 - Computación Científica II  </h1>
    <h2> Valores y Vectores Propios </h2>
    <h2> <a href="#acknowledgements"> [S]cientific [C]omputing [T]eam </a> </h2>
    <h2> Version: 1.14</h2>
</center>
Tabla de Contenidos
Introducción
Marco Teórico
Algoritmos e Implementaciones
Power Iteration
Inverse Power Iteration
Rayleigh Quotient Iteration
SciPy Eigenvalue
Problema de Aplicación
Acknowledgements
<div id='intro' />
Introducción
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).
<div id='teo' />
Marco Teórico
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:
python
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
<div id='alg' />
Algoritmos e Implementaciones
Librerías utilizadas durante la clase
End of explanation
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)
Explanation: Matriz y vector de prueba
End of explanation
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)    
# Testing algorithm
lam, u = power_iteration(A, x, 20, verbose=True)
print("lambda = {0}".format(lam))
print("u (dominant eigenvector) = {0}".format(u))
Explanation: <div id='pi' />
Power Iteration
A continuación se entrega el código del algoritmo clásico de Power Iteration. Pruebe cambiando las matrices y los parámetros del algoritmo.
End of explanation
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)
# 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))
Explanation: <div id='invpi' />
Inverse Power Iteration
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:
Los valores propios de la matriz inversa $A^{-1}$ son los recíprocos de los valores propios de $A$, es decir: $\lambda_1^{-1}, \lambda_2^{-1}, \ldots , \lambda_m^{-1}$. Los vectores propios de se mantienen inalterados.
Los valores propios de la matriz con shift $A - sI$ son: $\lambda_1-s, \lambda_2-s, \ldots, \lambda_m-s$. Del mismo modo, los vectores propios se mantienen inalterados.
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:
End of explanation
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)
Explanation: <div id='rq' />
Rayleigh Quotient Iteration
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:
La convergencia pasa a ser cuadrática (de modo general) y cúbica para matrices simétricas.
Sin embargo, se paga el costo de tener que resolver un sistema de ecuaciones diferentes en cada iteración.
A continuación se presenta una implementación del RQI:
End of explanation
# Testing algorithm
lam, v = rqi(A, x, k=2)
print("lambda = {0}".format(lam))
print("v = {0}".format(v))
Explanation: Preguntas: 
1. ¿Porque es necesario el try y except en las líneas 11 y 13? ¿Que significa que el sistema no pueda ser resuelto?
2. Como puede observar RQI no recibe shift como parámetro. ¿A cuál valor/vector propio convergerá? ¿Como forzar/guiar a que tienda hacia un valor/vector propio distinto?
End of explanation
# 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)
Explanation: <div id='sp' />
$\texttt{SciPy}$ Eigenvalue
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.
End of explanation
<END_TASK>