# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# <center>
#     <h1> ILI286 - Computación Científica II  </h1>
#     <h2> Valores y Vectores Propios </h2>
#     <h2> [[S]cientific [C]omputing [T]eam](#acknowledgements)</h2>
#     <h2> Version: 1.11</h2>
# </center>

# # Tabla de Contenidos
# * [Introducción](#intro)
# * [<NAME>](#teo)
# * [Algoritmos e Implementaciones](#alg)
#     * [Power Iteration](#pi)
#     * [Inverse Power Iteration](#invpi)
#     * [Rayleigh Quotient Iteration](#rq)
#     * [SciPy Eigenvalue](#sp)
#     * [Problema de Aplicación](#problema)
# * [Acknowledgements](#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' />
# # <NAME>
#
# 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{x^T A x}{x^T 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 $\{v_1, v_2, \ldots, 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 $v_k$:
#
# \begin{align}
# A x_0 &= c_1 A v_1 + \cdots + c_m A v_m = c_1 \lambda_1 v_1 + \cdots + c_m \lambda_m v_m \\
# A^2 x_0 & = c_1 \lambda_1 A v_1 + \cdots + c_m \lambda_m A v_m = c_1 \lambda_1^2 v_1 + \cdots + c_m \lambda_m^2 v_m \\
# \vdots &= \vdots \\
# A^k x_0 &= c_1 \lambda_1^k v_1 + \cdots + c_m \lambda_m^k v_m
# \end{align}
#
# Factorizando $\lambda_1^k$ del último resultado se obtiene:
#
# $$ \frac{A^k x_0}{\lambda_1^k} = c_1 v_1 + c_2 \left(\frac{\lambda_2}{\lambda_1}\right)^k v_2 + \cdots + c_m \left(\frac{\lambda_m}{\lambda_1}\right)^k 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

import numpy as np
from scipy import linalg
from matplotlib import pyplot as plt
# %matplotlib inline

# ### Matriz y vector de prueba

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)


# <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.

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))


# <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:
#
# 1. 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.
# 2. 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:

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)
        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))


# <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_:
#
# 1. La convergencia pasa a ser cuadrática (de modo general) y cúbica para matrices simétricas.
# 2. 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:

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:** 
# 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?

# Testing algorithm
lam, v = rqi(A, x, k=2)
print("lambda = {0}".format(lam))
print("v = {0}".format(v))

# <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.

# 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)

# <div id='problema' />
# ## Problema de Aplicación
#
# Las matrices simétricas tiene una propiedad muy interesante:
#
# * Los vectores propios de las matrices simétricas son ortogonales entre sí.
#
# En base a lo anterior se propone el siguiente algoritmo para encontrar los primeros $k$ valores/vectores propios:
#
# ```python
# 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    
# ```
#
# 1. Justifique la validez de tal procedimiento.
# 2. 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](https://en.wikipedia.org/wiki/Covariance_matrix):
# $$\Sigma_X = \frac{1}{n-1}X^T X$$
# donde $X \in \mathbb{R}^{m \times n}$, con $m$ _samples_ y $n$ _features_.
#
# 3. Concluya acerca de la utilidad del procedimiento propuesto. 

# <div id='acknowledgements' />
# # Acknowledgements
# * _Material creado por profesor <NAME>_ (`<EMAIL>`) _y ayudantes: <NAME> y <NAME>. DI UTFSM. Abril 2016._
#
# ***
# ### DISCLAIMER ###
#
# El presente notebook ha sido creado para el curso **ILI286 - Computación Científica 2**, del [Departamento de Informática](http://www.inf.utfsm.cl/), [Universidad Técnica Federico Santa María](http://www.utfsm.cl/). 
#
# El material ha sido creado por <NAME> <<EMAIL>> y <NAME> <<EMAIL>>, 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.
# ***

StarcoderdataJupyterScriptsDedupFiltered