<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step2:  (AAEVALEVEC)=
2.3 Algoritmos y aplicaciones de eigenvalores y eigenvectores de una matriz
```{admonition} Notas para contenedor de docker
Step3:  ```{margin}
Estos son los índices que corresponden a la entrada $a_{12}$ de $A$.
```
Step4:  ```{margin}
La matriz de rotación de Jacobi es en este caso de tamaño $2 \times 2$.
```
Step5:  Por construcción $A$ y $J^TAJ$ son similares y por tanto tienen el mismo espectro
Step7:  Ver eigvals y eigvalsh.
````{admonition} Comentarios
En matrices mayores a dos dimensiones el método de rotaciones de Jacobi para matrices simétricas requiere eliminación repetida (volver a hacer cercanas a cero) del par de entradas de $A$ seleccionadas de iteraciones previas pues en cada iteración se destruyen (vuelven a ser diferentes de cero) en general. 
El método de rotaciones de Jacobi para matrices simétricas utiliza como criterios de paro
Step8:  El primer sweep considera las entradas
Step9:  Entrada $a_{31}$ rotación en el plano $(1,3)$
Step10:  Entrada $a_{21}$ rotación en el plano $(1,2)$
Step11:  Entrada $a_{42}$ rotación en el plano $(2,4)$
Step12:  Entrada $a_{32}$ rotación en el plano $(2,3)$
Step13:  Entrada $a_{43}$ rotación en el plano $(3,4)$
Step14:  Segundo sweep
Entrada $a_{41}$ rotación en el plano $(1,4)$
Step15:  Entrada $a_{31}$ rotación en el plano $(1,3)$
Step16:  Entrada $a_{21}$ rotación en el plano $(1,2)$
Step17:  Entrada $a_{42}$ rotación en el plano $(2,4)$
Step18:  Entrada $a_{32}$ rotación en el plano $(2,3)$
Step19:  Entrada $a_{43}$ rotación en el plano $(3,4)$
Step20:  ````{margin}
Obsérvese que las entradas diferentes en la diagonal son más pequeñas que la matriz original $A$.
````
Step21:  ```{admonition} Ejercicio
Step22:  Se toma $T_0$ como $A$
Primera iteración
Step23:  ```{margin} 
$T_{k-1} = Q_kR_k$.
Ver numpy.linalg.qr para factorización QR con NumPy.
```
Step24:  ```{margin}
Matriz $Q$ ortogonal.
```
Step25:  ```{margin}
Matriz $R$ triangular superior.
```
Step26:  ```{margin}
$T_k = R_k Q_k$.
```
Step27:  Segunda iteración
Step28:  Tercera iteración
Step29:  ````{margin}
Obsérvese que las entradas diferentes en la diagonal son más pequeñas que la matriz original $A$.
````
Step30:  (COMENTARIOSALGQR)=
````{admonition} Comentarios
Si $A = Q\Lambda Q^T$, ver {ref}descomposición espectral de una matriz simétrica &lt;DESCESP&gt;, el método anterior converge a una matriz diagonal $\Lambda$ en la que se encuentran los eigenvalores de $A$ si se satisfacen las siguientes condiciones
Step31:  Primera iteración
Step32:  ```{margin}
$z^{(k)} = Az^{(k-1)}$.
```
Step33:  Segunda iteración
Step34:  ```{margin}
$z^{(k)} = Az^{(k-1)}$.
```
Step35:  Tercera iteración
Step36:  ...
Novena iteración
Step37:  Y es más cercano al eigenvector $\left [\begin{array}{c} 1\ 1 \end{array} \right ]$ (escalado por un factor).
Lo anterior funciona pues si $q_1, \dots, q_n$ son eigenvectores de $A$ y $z^{(0)} = \displaystyle \sum_{i=1}^n \alpha_i q_i$ (combinación lineal) entonces
Step38:  ```{margin}
Compárese con
Step39:  ```{admonition} Comentarios
El método de la potencia aproxima al eigenvector asociado al eigenvalor con módulo máximo. Esto en contraste con el algoritmo QR y el método de rotaciones de Jacobi que aproximan todos los eigenvalores y eigenvectores. En la práctica es posible que sólo se requieran unos cuantos eigenvalores-eigenvectores.
Típicamente funciona el método de la potencia pero puede fallar por razones como
Step40:  ```{margin}
Matriz que se le añadirá a $Q$, tiene en sus columnas el valor igual a $\frac{1}{6}$ para el ejemplo dado.
```
Step41:  ```{margin}
Se obtiene $P$ a partir de $Q + \frac{1}{n}ed^T$.
```
Step42:  ```{admonition} Comentario
Una propiedad que tienen las matrices estocásticas por columnas es $P^k$ es estocástica por columnas, $k \in \mathbb{N}$.
```
Ejemplo
Step43:  ```{margin}
Ver sección stochastic matrices and Markov chains del libro Matrix Analysis and Applied Linear Algebra para caso general de matrices que convergen a un vector de distribución estacionario y que tienen mismos valores en sus entradas por renglones o columnas.
```
Step44:  ```{admonition} Comentarios
Otras propiedades que tienen las matrices estocásticas por columnas son
Step45:  Una pregunta natural es si dada una matriz estocástica por columnas y un estado inicial $x_0$, los vectores de la cadena de Markov se estabilizan, esto es, la sucesión ${x_0, x_1, x_2, \dots }$ converge a algún vector $x$.
```{admonition} Definición
Dada una matriz estocástica por columna $P$ se le nombra $x$ un vector estacionario de P si y sólo si  $x$ es un vector de probabilidad con la propiedad
Step46:  Matriz irreducible y estocástica por columnas
Para asegurar que la matriz $P$ anterior tiene un único eigenvalor $\lambda_1$ igual a $1$ se utiliza la idea de teleportation y además se busca que la matriz sea irreducible.
```{admonition} Definición
Nombramos a una matriz irreducible como aquella para la que su grafo asociado siempre existe una ruta dirigida que une dos nodos distintos. Esto se nombra grafo fuertemente conectado.
```
```{margin}
En general la modificación es de la forma
Step47:  Tomando un valor de $\alpha = 0.85$ se tiene
Step48:  Y sus eigenvalores son
Step49:  sus eigenvectores son
Step50:  Verifiquemos que si definimos un vector de probabilidad $x_0$ con norma Euclidiana unitaria entonces al aplicar el método de la potencia
Step51:  Obsérvese que la aproximación es un múltiplo de la primer columna de la matriz siguiente
Step52:  Aproximación
Step53:  z_k es un vector de probabilidad
Step54:  Ejemplo
Imaginemos tenemos $4$ páginas $web$ conectadas como se muestra a continuación
Step55:  Definamos el vector de probabilidad inicial para aplicar el método de la potencia. Para ello, supongamos que la surfer visita aleatoriamente alguna de las páginas con la misma probabilidad, esto es
Step56:  $Q$ en una matriz irreducible y estocástica con la que podríamos trabajar, lo siguiente se realiza sólo para mostrar la aplicación de la transformación descrita.
Step57:  Aplicamos el método de la potencia
Step58:  z_k es un vector de probabilidad
Step59:  Verificamos que se aproxima al eigenvector asociado al eigenvalor igual a $1$
Step60:  Obsérvese que la aproximación es un múltiplo de la primer columna de la matriz anterior
Step61:  Interpretación
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(precision=3, suppress=True)
def sign(x):
    
    Helper function for computing sign of real number x.
    
    if x >=0:
        return 1
    else:
        return -1
def compute_cos_sin_Jacobi_rotation(Ak, idx1, idx2):
    
    Helper function for computing entries of Jacobi rotation.
    Args:
        Ak (numpy ndarray): Matrix of iteration k in Jacobi rotation method.
        idx1 (int): index for rows in Jacobi rotation matrix.
        idx2 (int): index for columns in Jacobi rotation matrix.
    Returns:
        c (float): value of cos of theta for Jacobi rotation matrix.
        s (float): value of sin of theta for Jacobi rotation matrix.
    
    if np.abs(Ak[idx1,idx2]) > np.finfo(float).eps:
        tau = (Ak[idx2, idx2] - Ak[idx1, idx1])/(2*Ak[idx1, idx2])
        t_star = sign(tau)/(np.abs(tau) + np.sqrt(1+tau**2))
        c = 1/np.sqrt(1+t_star**2)
        s = c*t_star
    else: #no rotation is performed
        c = 1
        s = 0
    return (c,s)
A = np.array([[1,2],
              [2,1]])
A_k = A.copy()
print(A_k)
Explanation: (AAEVALEVEC)=
2.3 Algoritmos y aplicaciones de eigenvalores y eigenvectores de una matriz
```{admonition} Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar &lt;ruta a mi directorio&gt; por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker y &lt;versión imagen de docker&gt; por la versión más actualizada que se presenta en la documentación.
docker run --rm -v &lt;ruta a mi directorio&gt;:/datos --name jupyterlab_optimizacion -p 8888:8888 -d palmoreck/jupyterlab_optimizacion:&lt;versión imagen de docker&gt;
password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_optimizacion
Documentación de la imagen de docker palmoreck/jupyterlab_optimizacion:&lt;versión imagen de docker&gt; en liga.
```
Nota generada a partir de liga.
```{admonition} Al final de esta nota la comunidad lectora:
:class: tip
Pondrá en uso los conocimientos revisados en las notas de {ref}operaciones y transformaciones básicas del Álgebra Lineal Numérica &lt;OTBALN&gt; y {ref}eigenvalores y eigenvectores &lt;EIGVALEIGVEC&gt; para comprender los algoritmos básicos y típicos que se han desarrollado en el análisis numérico para el cálculo de eigenvalores y eigenvectores de una matriz simétrica.
Comprenderá la relación entre el método de pagerank para el ordenamiento o ranking de una página web utilizada por Google y el problema del cálculo de eigenvalores y eigenvectores de una matriz.
```
Para los algoritmos revisados en esta nota asumimos que la matriz $A \in \mathbb{R}^{n \times n}$ es simétrica. Esto tiene como consecuencia que sus eigenvalores son reales y existe una base de eigenvectores ortonormal, ver {ref}descomposición espectral de una matriz simétrica &lt;DESCESP&gt;. Para la sección del método de pagerank se asume una matriz sin una estructura particular.
(ROTJACMATSIM)=
Rotaciones de Jacobi para matrices simétricas $A \in \mathbb{R}^{n \times n}$
Este método produce una secuencia de transformaciones ortogonales de la forma $J_k^TAJ_k$ con el objetivo de hacer "más diagonal" a la matriz $A \in \mathbb{R}^{n \times n}$.
Si la matriz $A$ es simétrica y $J_0$ es una transformación de rotación de Jacobi, ver {ref}transformaciones de rotación &lt;TROT&gt;, entonces el esquema iterativo:
$$A_{k+1} = (J_0 J_1 \cdots J_k)^T A(J_0 J_1 \cdots J_k)$$
converge a una matriz diagonal en la que se encuentran los eigenvalores de $A$.
```{margin}
Los pasos de un algoritmo representan una guía para la implementación, no implica que se tengan que implementar uno a continuación del otro como se describe. Si una implementación respeta la lógica y al mismo algoritmo, entonces pueden seguirse los pasos de una forma distinta.
```
```{margin}
Los criterios de paro en un método iterativo ayudan a contestar preguntas del tipo ¿hemos resuelto el problema de forma aproximada?
```
Algoritmo: método de rotaciones de Jacobi para matrices simétricas
Dados $A$ simétrica y $tol >0$ definir $A_0 = A$, $Q_0 = I_n$.
Repetir el siguiente bloque para $k=0,1,2,\dots$
Elegir un par de índices $(idx1,idx2)$ con alguna de las metodologías descritas en el bloque siguiente de comentarios.
Calcular las entradas $\cos(\theta),\sin(\theta)$ de la matriz de rotación $J_k$.
$A_{k+1} = J_k^T A_k J_k$
$Q_{k+1} = Q_{k}J_k$.
hasta convergencia: satisfacer criterio de paro en el que se utiliza $tol$ y $maxsweeps$.
La matriz $J_k$ se utiliza para eliminar un par de entradas (simétricas) en la matriz $A_k$, esto  preserva la simetría de la matriz original. En las columnas de la matriz $Q_{k}$ se encuentran aproximaciones a los eigenvectores de $A$ y en la diagonal de $A_{k}$ se tienen aproximaciones a los eigenvalores de $A$.
```{admonition} Comentarios
En el método se hace mención de metodologías que ayudan a elegir los índices del renglón y columna del par de entradas de $A$ que serán eliminadas (hacer cercanas a cero). Algunas de éstas son:
1.Elegir $(idx1,idx2)$ tales que $|a_{idx1,idx2}| = \displaystyle \max_{i \neq j}|a_{ij}|$.
2.Ordenamiento cíclico por renglones: elegir $(idx1, idx2)$ en el conjunto $(1,2),(1,3),\dots,(1,n),(2,3),(2,4)\dots,(n-1,n)$.
```
```{admonition} Observación
:class: tip
Obsérvese que $A_{k+1}$ y $A_{k}$ son matrices ortogonalmente similares, ver {ref}similitud &lt;SIMILITUD&gt;.
```
Para encontrar la forma que debe tener $J_k$ es suficiente considerar el caso $2 \times 2$ y se asume que $a_{12} \neq 0$ pues si $a_{12} =0$ entonces no hay que realizar rotación:
$$
\begin{eqnarray}
J_k^TAJ_k &=&
\left [
\begin{array}{cc}
c & -s\
s & c
\end{array}
\right ]
\left [
\begin{array}{cc}
a_{11} & a_{12}\
a_{12} & a_{22}
\end{array}
\right ]
\left [
\begin{array}{cc}
c & s\
-s & c
\end{array}
\right ]
\nonumber \
&=&
\left [
\begin{array}{cc}
c^2a_{11} -2 csa_{12} + s^2a_{22} & c^2a_{12} - cs(a_{22}-a_{11})-s^2a_{12}\nonumber \
c^2a_{12} - cs(a_{22}-a_{11})-s^2a_{12} & c^2a_{22}+2csa_{12}+s^2a_{11} \nonumber
\end{array}
\right ]
\nonumber 
\end{eqnarray}
$$
donde: $c$ y $s$ representan a $\cos(\theta), \sin(\theta)$ respectivamente y $\theta$ ángulo para rotar.
Si se desea que la entrada $(1,2)$ (equivalentemente por simetría la $(2,1)$) sea cero se debe cumplir:
$$c^2a_{12} - cs(a_{22}-a_{11})-s^2a_{12}=0.$$
Asignando la variable $t = \frac{s}{c}$ (tangente de $\theta$) se obtiene la ecuación cuadrática:
$$1 - t\frac{(a_{22}-a_{11})}{a_{12}} - t^2 = 0.$$
Equivalentemente:
$$t^2 + t\frac{(a_{22}-a_{11})}{a_{12}} - 1 = t^2 + 2\tau t -1 = 0$$
```{margin}
Las funciones $\frac{1}{\tau+\sqrt{\tau^2+1}}$, $\frac{1}{\tau-\sqrt{\tau^2+1}}$ son estrictamente decrecientes para $\tau \geq 0$ y $\tau < 0$ respectivamente. Un valor de $\tau$ corresponde a un único ángulo $\theta \in \left [-\frac{\pi}{4}, \frac{\pi}{4} \right ]$.
```
donde: $\tau = \frac{a_{22}-a_{11}}{2a_{12}}$. Las raíces de la ecuación anterior están dadas por:
$$
\begin{eqnarray}
t_1^ &=& -\tau + \sqrt{\tau^2+1}&=&\frac{1}{\tau+\sqrt{\tau^2+1}} \nonumber \
t_2^ &=& -\tau - \sqrt{\tau^2+1}&=&\frac{1}{\tau-\sqrt{\tau^2+1}} \nonumber
\end{eqnarray}
$$
Se sugiere utilizar la raíz de menor magnitud para disminuir errores por redondeo por lo que:
$$t^* = \frac{\text{signo}(\tau)}{|\tau| + \sqrt{1+\tau^2}}$$
donde:
$$
\text{signo}(x) = \begin{cases}
1 \text{ si } x>=0\
-1 \text{ en otro caso}
\end{cases}
$$
Las relaciones entre coseno, seno y tangente permiten obtener sus valores correspondientes:
$$c = \frac{1}{\sqrt{1+t^{*2}}},$$
$$s = ct^*$$
y así tener completamente definida a la matriz $J_k$.
Ejemplo
Considera: 
$$
A = 
\left [
\begin{array}{cc}
1&2\
2&1\
\end{array}
\right ].
$$
Eliminar las entradas $(1,2)$ y $(2,1)$ con una matriz $J$ de rotación de Jacobi:
End of explanation
idx_1 = 0
idx_2 = 1
c, s = compute_cos_sin_Jacobi_rotation(A_k, idx_1, idx_2)
Explanation: ```{margin}
Estos son los índices que corresponden a la entrada $a_{12}$ de $A$.
```
End of explanation
J = np.array([[c, s],
              [-s, c]])
print(J.T@A@J)
Explanation: ```{margin}
La matriz de rotación de Jacobi es en este caso de tamaño $2 \times 2$.
```
End of explanation
print(np.linalg.eigvalsh(A))
print(np.linalg.eigvalsh(J.T@A@J))
Explanation: Por construcción $A$ y $J^TAJ$ son similares y por tanto tienen el mismo espectro:
End of explanation
def compute_Jacobi_rotation(Ak, idx1, idx2):
    
    Compute Jacobi rotation matrix.
    Args:
        Ak (numpy ndarray): Matrix of iteration k in Jacobi rotation method.
        idx1 (int): index for rows in Jacobi rotation matrix.
        idx2 (int): index for columns in Jacobi rotation matrix.
    Returns:
        J (numpy ndarray): Jacobi rotation matrix.
    
    c,s = compute_cos_sin_Jacobi_rotation(Ak, idx1, idx2)
    m,n = Ak.shape
    J = np.eye(m)
    J[idx1, idx1] = J[idx2, idx2] = c
    J[idx1, idx2] = s
    J[idx2, idx1] = -s
    return J
A = np.array([[1,2,3,4],
              [2,-2,4,5],
              [3, 4, 6, 7],
              [4, 5, 7,-8.0]])
print(A)
A_k = A.copy()
Explanation: Ver eigvals y eigvalsh.
````{admonition} Comentarios
En matrices mayores a dos dimensiones el método de rotaciones de Jacobi para matrices simétricas requiere eliminación repetida (volver a hacer cercanas a cero) del par de entradas de $A$ seleccionadas de iteraciones previas pues en cada iteración se destruyen (vuelven a ser diferentes de cero) en general. 
El método de rotaciones de Jacobi para matrices simétricas utiliza como criterios de paro:
La cantidad $\text{off}(A) = \sqrt{\displaystyle \sum_{i=1}^n \sum_{j=1, j\neq i}^n a_{ij}^2}$ (norma de Frobenius de $A$ sin la diagonal, ver {ref}Normas matriciales &lt;NMAT&gt;).
Número máximo de sweeps. Un sweep es igual a $\frac{n(n-1)}{2}$ y corresponde al número máximo de entradas de la matriz que forman la parte triangular superior de $A$ sin contar a la diagonal que se asumen diferentes de cero. No existe teoría rigurosa para el número de sweeps pero una heurística encontrada por Brent y Luk, p. 13, 1985 menciona que el número máximo es proporcional a $\mathcal{O}(\log(n))$ y en la práctica se utilizan entre $6$ y $10$, ver H. Rutishauser, The Jacobi method for real symmetric matrices, 1966.
El criterio de paro es de la forma:
while off(A_k) &gt; tol ||A_k||_F &amp;&amp; sweeps &lt; max_sweeps
con $tol$ una cantidad pequeña y positiva (comúnmente menor o igual a $10^{-8}$), ||A_k||_F norma de Frobenius de $A_k$ y sweeps contador de los sweeps.
````
Ejemplo
Considera: 
$$
A = 
\left [
\begin{array}{cccc}
1 & 2 & 3 & 4\
2 & -2 & 4 & 5\
3 & 4 & 6 & 7\
4 & 5 & 7 & -8
\end{array}
\right ].
$$
Utilizando la metodología de ordenamiento cíclico por renglones realizar dos sweeps del método de rotaciones de Jacobi para matrices simétricas.
End of explanation
idx_1 = 0
idx_2 = 3
J = compute_Jacobi_rotation(A_k, idx_1, idx_2)
A_k = J.T@A_k@J
print(A_k)
Explanation: El primer sweep considera las entradas: $(1,2), (1,3), (1,4), (2,3), (2,4), (3,4)$. 
A continuación se toman en el orden $(1,4), (1,3), (1,2), (2,4), (2,3)$ y finalmente $(3,4)$ para mostrar que el orden de selección de las entradas en una forma decreciente no importa.
```{margin}
Se sugiere haber revisado los ejemplos de {ref}transformaciones de rotación &lt;TROT&gt;.
```
Entrada $a_{41}$ rotación en el plano $(1,4)$:
End of explanation
idx_1 = 0
idx_2 = 2
J = compute_Jacobi_rotation(A_k, idx_1, idx_2)
A_k = J.T@A_k@J
print(A_k)
Explanation: Entrada $a_{31}$ rotación en el plano $(1,3)$:
End of explanation
idx_1 = 0
idx_2 = 1
J = compute_Jacobi_rotation(A_k, idx_1, idx_2)
A_k = J.T@A_k@J
print(A_k)
Explanation: Entrada $a_{21}$ rotación en el plano $(1,2)$:
End of explanation
idx_1 = 1
idx_2 = 3
J = compute_Jacobi_rotation(A_k, idx_1, idx_2)
A_k = J.T@A_k@J
print(A_k)
Explanation: Entrada $a_{42}$ rotación en el plano $(2,4)$:
End of explanation
idx_1 = 1
idx_2 = 2
J = compute_Jacobi_rotation(A_k, idx_1, idx_2)
A_k = J.T@A_k@J
print(A_k)
Explanation: Entrada $a_{32}$ rotación en el plano $(2,3)$:
End of explanation
idx_1 = 2
idx_2 = 3
J = compute_Jacobi_rotation(A_k, idx_1, idx_2)
A_k = J.T@A_k@J
print(A_k)
Explanation: Entrada $a_{43}$ rotación en el plano $(3,4)$:
End of explanation
idx_1 = 0
idx_2 = 3
J = compute_Jacobi_rotation(A_k, idx_1, idx_2)
A_k = J.T@A_k@J
print(A_k)
Explanation: Segundo sweep
Entrada $a_{41}$ rotación en el plano $(1,4)$:
End of explanation
idx_1 = 0
idx_2 = 2
J = compute_Jacobi_rotation(A_k, idx_1, idx_2)
A_k = J.T@A_k@J
print(A_k)
Explanation: Entrada $a_{31}$ rotación en el plano $(1,3)$:
End of explanation
idx_1 = 0
idx_2 = 1
J = compute_Jacobi_rotation(A_k, idx_1, idx_2)
A_k = J.T@A_k@J
print(A_k)
Explanation: Entrada $a_{21}$ rotación en el plano $(1,2)$:
End of explanation
idx_1 = 1
idx_2 = 3
J = compute_Jacobi_rotation(A_k, idx_1, idx_2)
A_k = J.T@A_k@J
print(A_k)
Explanation: Entrada $a_{42}$ rotación en el plano $(2,4)$:
End of explanation
idx_1 = 1
idx_2 = 2
J = compute_Jacobi_rotation(A_k, idx_1, idx_2)
A_k = J.T@A_k@J
print(A_k)
Explanation: Entrada $a_{32}$ rotación en el plano $(2,3)$:
End of explanation
idx_1 = 2
idx_2 = 3
J = compute_Jacobi_rotation(A_k, idx_1, idx_2)
A_k = J.T@A_k@J
Explanation: Entrada $a_{43}$ rotación en el plano $(3,4)$:
End of explanation
print(A)
print(A_k)
Explanation: ````{margin}
Obsérvese que las entradas diferentes en la diagonal son más pequeñas que la matriz original $A$.
````
End of explanation
A = np.array([[7, 2],
              [2, 4]])
Explanation: ```{admonition} Ejercicio
:class: tip
Realizar dos sweeps más para la matriz simétrica del ejemplo anterior y sólo dos sweeps para la matriz:
$$
A = 
\left [
\begin{array}{ccc}
1 & 0 & 2\
0 & 2 & 1\
2 & 1 & 1
\end{array}
\right ]
$$
Para ambos ejercicios calcular en cada iteración realizada $\text{off}(A)$.
```
```{margin}
Los algoritmos que pueden ejecutarse en una forma por bloques se les nombra blocking algorithms. Tales algoritmos aprovechan niveles altos de BLAS en su implementación.
```
```{admonition} Comentarios
Algunas características del método de rotaciones de Jacobi para matrices simétricas:
Tiene una convergencia cuadrática, esto es, la matriz $A_k$ al ir aproximándose hacia una diagonal, el número de dígitos de precisión en cada iteración se duplica. Sin embargo el algoritmo QR (versión no simple) es más eficiente en cuanto al número de iteraciones.
La precisión que tiene en el cálculo de eigenvalores es alta (error relativo pequeño).
Es un algoritmo que puede realizarse con cómputo en paralelo si se utilizan ordenamientos específicos como el cíclico por renglones. También este algoritmo puede realizarse por bloques.
```
(ALGQR)=
Algoritmo QR o QR iteration (versión simple) para matrices simétricas $A \in \mathbb{R}^{n \times n}$
Este método produce una secuencia de transformaciones ortogonales de la forma $Q_k^TAQ_k$ con el objetivo de hacer "más diagonal" a la matriz $A \in \mathbb{R}^{n \times n}$.
Si la matriz $A$ es simétrica y $Q_0$ es una matriz con columnas linealmente independientes entonces el esquema iterativo:
$$T_k = (Q_0 Q_1 \cdots Q_k)^T A(Q_0 Q_1 \cdots Q_k)$$
casi siempre converge a una matriz diagonal en la que se encuentran los eigenvalores de $A$, ver {ref}comentarios algoritmo QR &lt;COMENTARIOSALGQR&gt; para suposiciones que garantizan la convergencia.
Factorización QR
Sólo en esta subsección asumimos $A \in \mathbb{R}^{m \times n}$ con $m >= n$. 
Existen $Q \in \mathbb{R}^{m \times m}$ ortogonal y $R \in \mathbb{R}^{m \times n}$ triangular superior tales que $A=QR$:
<img src="https://dl.dropboxusercontent.com/s/92uiu8rwi64he0t/A_QR.png?dl=0" heigth="350" width="350">
```{admonition} Comentarios
La factorización $QR$ puede escribirse en una forma "delgada" obteniéndose la factorización thin QR o *reduced* :
<img src="https://dl.dropboxusercontent.com/s/u0u365q5372xihl/A_thin_QR.png?dl=0" heigth="300" width="300">
En este caso $Q_1$ no se nombra ortogonal sino con columnas ortonormales, ver {ref}matriz ortogonal y matriz con columnas ortonormales &lt;MATORTMATCOLORTONO&gt;, $R_1$ es triangular superior. Si el rank(A) es igual a $n$ ($A$ se nombra de rank completo o full rank) y $R_1$ tiene entradas en la diagonal positivas, la factorización thin QR es única.
La factorización thin o reduced $QR$ de $A$ realiza un cambio de bases entre las columnas de $A$ y las columnas de $Q_1$ teniendo las entradas de $R$ las coordenadas del cambio entre bases.
La factorización QR puede calcularse con reflexiones de Householder o rotaciones Givens, ver {ref}transformaciones de reflexión &lt;TREF&gt; y {ref}transformaciones de rotación &lt;TROT&gt;.
```
Algoritmo: algoritmo QR para matrices simétricas (versión simple)
Dados $A \in \mathbb{R}^{n \times n}$ simétrica, $tol > 0$, $Q_0 \in \mathbb{R}^{n \times n}$ ortogonal definir $T_0 = Q_0^TAQ_0$.
Repetir el siguiente bloque para $k=1,2,\dots$
$T_{k-1} = Q_kR_k$ (factorización $QR$ de $T_{k-1}$).
$T_k = R_k Q_k$
hasta convergencia: satisfacer criterio de paro en el que se utiliza $tol$.
````{admonition} Observaciones
:class: tip
Obsérvese que $T_k = R_kQ_k = Q_k^T(Q_k R_k)Q_k = Q_k ^TT_{k-1}Q^k$ por lo que $T_k$ y $T_{k-1}$ son ortogonalmente similares, ver {ref}similitud &lt;SIMILITUD&gt;.
Es fácil verificar que $T_k = (Q_0 Q_1 \cdots Q_k)^T A(Q_0 Q_1 \cdots Q_k)$ por lo que $T_k$ es ortogonalmente similar a $A$.
````
Ejemplo
End of explanation
T_k_minus_1 = A
Explanation: Se toma $T_0$ como $A$
Primera iteración
End of explanation
q,r = np.linalg.qr(T_k_minus_1)
Explanation: ```{margin} 
$T_{k-1} = Q_kR_k$.
Ver numpy.linalg.qr para factorización QR con NumPy.
```
End of explanation
print(q)
Explanation: ```{margin}
Matriz $Q$ ortogonal.
```
End of explanation
print(r)
Explanation: ```{margin}
Matriz $R$ triangular superior.
```
End of explanation
T_k = r@q
print(T_k)
Explanation: ```{margin}
$T_k = R_k Q_k$.
```
End of explanation
T_k_minus_1 = T_k
q,r = np.linalg.qr(T_k_minus_1)
print(q)
print(r)
T_k = r@q
print(T_k)
Explanation: Segunda iteración
End of explanation
T_k_minus_1 = T_k
q,r = np.linalg.qr(T_k_minus_1)
print(q)
print(r)
T_k = r@q
Explanation: Tercera iteración
End of explanation
print(A)
print(T_k)
print(np.linalg.eigvals(A))
Explanation: ````{margin}
Obsérvese que las entradas diferentes en la diagonal son más pequeñas que la matriz original $A$.
````
End of explanation
A = np.array([[1.5, 0.5],
              [0.5, 1.5]])
Explanation: (COMENTARIOSALGQR)=
````{admonition} Comentarios
Si $A = Q\Lambda Q^T$, ver {ref}descomposición espectral de una matriz simétrica &lt;DESCESP&gt;, el método anterior converge a una matriz diagonal $\Lambda$ en la que se encuentran los eigenvalores de $A$ si se satisfacen las siguientes condiciones:
$A \in \mathbb{R}^{n \times n}$ es simétrica.
$|\lambda_1| > |\lambda_2| > \cdots > |\lambda_n|.$
El coseno del ángulo entre las columnas de $Q_0$ y las columnas de $Q$ que contienen los eigenvectores de $A$ es diferente de cero (no ortogonales), ver {ref}normas vectoriales y matriciales &lt;NVM&gt;. 
El criterio de paro es de la forma:
while off(T_k) &gt; tol ||T_k||_F &amp;&amp; iterations &lt; max_iters
con tol una cantidad pequeña y positiva (comúnmente menor o igual a $10^{-8}$), ||T_k||_F norma de Frobenius de $T_k$, iterations un contador de las iteraciones y $\text{off}(T) = \sqrt{\displaystyle \sum_{i=1}^n \sum_{j=1, j\neq i}^n t_{ij}^2}$ (norma de Frobenius de $T$ sin la diagonal, ver {ref}Normas matriciales &lt;NMAT&gt;).
````
```{admonition} Ejercicio
:class: tip
Realiza tres iteraciones del algoritmo QR para aproximar los eigenvalores de la matriz:
$$
A = 
\left [
\begin{array}{cccc}
8 & 4 & 2 & 1 \
4 & 8 & 2 & 1 \
2 & 2 & 8 & 1 \
1 & 1 & 1 & 8
\end{array}
\right ]
$$
```
(MPOTMATSIM)=
Método de la potencia para matrices simétricas $A \in \mathbb{R}^{n \times n}$
Si la matriz $A$ tiene un único eigenvalor $\lambda_1$ que tiene máximo módulo y $z^{(0)}$ es un vector distinto de cero entonces el esquema iterativo:
$$z^{(k)} = A z^{(k-1)}$$
casi siempre converge a un múltiplo del eigenvector de $A$ asociado a $\lambda_1$, ver {ref}comentarios del método de la potencia &lt;COMENTARIOSMETPOTENCIA&gt; para suposiciones que garantizan la convergencia.
Ejemplo
End of explanation
z_k_minus_1 = np.array([0,1])
Explanation: Primera iteración
End of explanation
z_k = A@z_k_minus_1
print(z_k)
Explanation: ```{margin}
$z^{(k)} = Az^{(k-1)}$.
```
End of explanation
z_k_minus_1 = z_k
Explanation: Segunda iteración
End of explanation
z_k = A@z_k_minus_1
print(z_k)
Explanation: ```{margin}
$z^{(k)} = Az^{(k-1)}$.
```
End of explanation
z_k_minus_1 = z_k
z_k = A@z_k_minus_1
print(z_k)
Explanation: Tercera iteración
End of explanation
z_k = np.array([127.5, 128.5])
print(z_k)
Explanation: ...
Novena iteración
End of explanation
A = np.array([[1.5, 0.5],
              [0.5, 1.5]])
n = 2
q_k = np.array([0,1])
max_iter = 7
q_k_iter = np.zeros((n, max_iter))
lambda_k_iter = np.zeros(max_iter)
for k in range(max_iter):
    z_k = A@q_k
    q_k = z_k/np.linalg.norm(z_k)
    string = "$q^{(%s)}$"%k
    plt.annotate(string,(q_k[0],q_k[1]),fontsize=10)
    q_k_plot = np.row_stack((np.zeros(n), q_k))
    q_k_iter[:,k] = q_k
    lambda_k = q_k.T@A@q_k
    lambda_k_iter[k] = lambda_k
    plt.plot(q_k_plot[:,0], q_k_plot[:,1])
plt.axhline(color='k')
plt.axvline(color='k')
plt.title("Método de la potencia normalizando: cada vector aproximación tiene norma Euclidiana 1")
plt.show()
Explanation: Y es más cercano al eigenvector $\left [\begin{array}{c} 1\ 1 \end{array} \right ]$ (escalado por un factor).
Lo anterior funciona pues si $q_1, \dots, q_n$ son eigenvectores de $A$ y $z^{(0)} = \displaystyle \sum_{i=1}^n \alpha_i q_i$ (combinación lineal) entonces:
$$
\begin{eqnarray}
z^{(k)} &=& A z^{(k-1)} = A^2 z^{(k-2)} = \cdots = A^k z^{(0)} \nonumber \
&=& A^k \displaystyle \sum_{i=1}^n \alpha_i q_i = \displaystyle \sum_{i=1}^n \alpha_i A^k q_i = \displaystyle \sum_{i=1}^n \alpha_i \lambda_i ^k q_i \nonumber \
&=& \lambda_1^k \left (\alpha_1 q_1 + \displaystyle \sum_{i=2}^n \left ( \frac{\lambda_i}{\lambda_1} \right ) ^k \alpha_i q_i \right ) \nonumber
\end{eqnarray}
$$
Y como $\left | \frac{\lambda_i}{\lambda_1} \right | < 1$ para $i > 1$ entonces al ir avanzando en las iteraciones lo anterior converge a un múltiplo de $q_1$.
```{admonition} Observación
:class: tip
Se añade normalización del eigenvector en cada iteración para evitar problemas de redondeo como overflow y underflow.
```
```{margin}
Este mismo algoritmo puede utilizarse para matrices no simétricas. La normalización puede ser con la norma $2$ o norma infinito (incluso funciona la componente que tiene máxima magnitud tomando todo y su signo).
```
Algoritmo: método de la potencia para matrices simétricas
Dados $A \in \mathbb{R}^{n \times n}$ simétrica, $q^{(0)} \in \mathbb{R}^n$ de norma Euclidiana unitaria y $tol > 0$.
Repetir el siguiente bloque para $k=1,2,\dots$
$z^{(k)} = Aq^{(k-1)}$.
$q^{(k)} = \frac{z^{(k)}}{||z^{(k)}||_2}$.
$\lambda ^{(k)} = q^{(k)T} Aq^{(k)}$.
hasta convergencia: satisfacer criterio de paro en el que se utiliza $tol$.
```{margin}
En el caso no simétrico se requiere que la matriz sea diagonalizable para la convergencia del método, ver {ref}matriz diagonalizable &lt;DIAGONALIZABLE&gt;.
```
(COMENTARIOSMETPOTENCIA)=
````{admonition} Comentarios
El método anterior converge de forma lineal al eigenvalor de máximo módulo $\lambda_1$ y a un eigenvector asociado a tal eigenvalor, $q_1$, si se satisfacen las siguientes condiciones:
$A \in \mathbb{R}^{n \times n}$ es simétrica.
$\lambda_1$ cumple: $|\lambda_1| > |\lambda_2| \geq \cdots \geq |\lambda_n|$ esto es, $\lambda_1$ es el único eigenvalor de máximo módulo.
El coseno del ángulo entre $q^{(0)}$ y $q_1$ es diferente de cero (no ortogonales), ver {ref}normas vectoriales y matriciales &lt;NVM&gt;. Esto no es un requerimiento fuerte pues al usar un sistema de punto flotante y tener errores por redondeo garantiza que en la mayoría de los casos prácticos tal coseno es diferente de cero.
El criterio de paro es de la forma: 
while ||Aq_k - lambda_kq_k||_2 &gt; tol &amp;&amp; iterations &lt; max_iters
con tol una cantidad pequeña y positiva (comúnmente menor o igual a $10^{-8}$), iterations un contador de iteraciones.
````
Ejemplo
End of explanation
eigval, eigvec = np.linalg.eig(A)
print("eigenvalor:")
print(eigval)
print("eigenvectores:")
print(eigvec)
print(q_k)
print(lambda_k_iter)
Explanation: ```{margin}
Compárese con:
```
End of explanation
l = [[0, 1/3, 0, 0, 0, 0],
     [1/3, 0, 0, 0, 0, 0],
     [0, 1/3, 0, 0, 1/3, 1/2],
     [1/3, 0, 0, 0, 1/3, 0],
     [1/3, 1/3, 0, 0, 0, 1/2],
     [0, 0, 1, 0, 1/3, 0]]
Q = np.array(l)
print(Q)
_,n = Q.shape
e = np.ones(n)
d = np.array([0, 0, 0, 1, 0, 0]) #definition of d
Explanation: ```{admonition} Comentarios
El método de la potencia aproxima al eigenvector asociado al eigenvalor con módulo máximo. Esto en contraste con el algoritmo QR y el método de rotaciones de Jacobi que aproximan todos los eigenvalores y eigenvectores. En la práctica es posible que sólo se requieran unos cuantos eigenvalores-eigenvectores.
Típicamente funciona el método de la potencia pero puede fallar por razones como:
No hay convergencia a un eigenvalor-eigenvector con entradas en $\mathbb{C}$ a partir de una matriz con entradas en $\mathbb{R}$ (pero si se utiliza una matriz $A$ y vectores con entradas en $\mathbb{C}$ se recupera la convergencia).
Puede haber más de un eigenvalor que tenga el mismo módulo y la convergencia podría darse para una combinación lineal de los eigenvectores asociados a tales eigenvalores.
La ventaja que tiene el método de la potencia es que sólo requiere una operación del tipo matriz-vector lo cual para matrices grandes y ralas (también nombradas sparse: que tienen muchas entradas iguales a cero) es ventajoso.
Sólo obtiene al eigenvalor con máximo módulo $\lambda_1$ y un eigenvector asociado, por lo que si se desean otras parejas de eigenvalores-eigenvectores se debe recurrir a métodos adicionales por ejemplo deflation.
La convergencia de este método depende mucho de la cercanía entre $|\lambda_1|$ y $|\lambda_2|$: si son muy cercanos la convergencia es lenta. Esto es posible mejorar al utilizar shifts.
```
```{admonition} Ejercicio
:class: tip
Realiza tres iteraciones del método de la potencia para calcular $\lambda_{\text{max}}(A)$ de la matriz:
$$
\left [
\begin{array}{ccc}
4 & -1 & 1 \
-1 & 3 & -2 \
1 & -2 & 3
\end{array}
\right ]
$$
con vector inicial 
$$ q^{(0)}= \left [
\begin{array}{c}
1 \
0 \
0
\end{array}
\right ]
$$
```
Iteración inversa o método de la potencia inversa para matrices simétricas $A \in \mathbb{R}^{n \times n}$
Este método originalmente obtiene al eigenvalor con módulo mínimo con su eigenvector asociado, sin embargo en una forma más general si se utiliza en el método de la potencia la matriz $(A-\mu I)^{-1}$ en lugar de la matriz $A$, con $\mu \in \mathbb{R}$ aproximación a un eigenvalor $\lambda$ de $A$, entonces se obtiene el método de la iteración inversa o de la potencia inversa:
```{margin}
Al igual que para el método de la potencia, este método funciona para matrices no simétricas si la matriz $A$ es diagonalizable.
```
```{margin}
La factorización $LU$ de $A$ puede calcularse con scipy.linalg.lu. Tal factorización utiliza {ref}transformaciones de Gauss &lt;TGAUSS&gt;.
```
Algoritmo: iteración inversa o método de la potencia inversa para matrices simétricas
Dados $A \in \mathbb{R}^{n \times n}$ simétrica, $q^{(0)} \in \mathbb{R}^n$ de norma Euclidiana unitaria, $\mu \in \mathbb{R}$ aproximación a un eigenvalor $\lambda$ de $A$ y $tol > 0$.
Calcular la factorización $LU$ de $A-\mu I$.
Repetir el siguiente bloque para $k=1,2,\dots$
Resolver el sistema de ecuaciones lineales $(A-\mu I)z^{(k)} = q^{(k-1)}$ para $z^{(k)}$ utilizando los factores $LU$ del paso 1.
$q^{(k)} = \frac{z^{(k)}}{||z^{(k)}||_2}$.
hasta convergencia: satisfacer criterio de paro en el que se utiliza $tol$.
````{admonition} Comentarios
Una vez que se conoce una aproximación a un eigenvalor, este método es muy útil para calcular su eigenvector asociado.
La convergencia del método anterior está asegurada si tal eigenvalor está bien separado de los demás.
Al igual que en el método de la potencia si se desea aproximar un eigenvalor de $A$ cercano a otro entonces la convergencia puede ser lenta.
El criterio de paro es de la forma: 
while ||Aq_k - lambda_kq_k||_2 &gt; tol &amp;&amp; iterations &lt; max_iters
con tol una cantidad pequeña y positiva (comúnmente menor o igual a $10^{-8}$), iterations un contador de iteraciones, lambda_k = q_k ^T A q_k.
````
(ITERCRAYMATSIM)=
Iteración por el cociente de Rayleigh para matrices simétricas $A \in \mathbb{R}^{n \times n}$
El cociente de Rayleigh de $x$ se define como:
$$
r(x) = \frac{x^TAx}{x^Tx}
$$
```{admonition} Ejercicio
:class: tip
Verificar que el cociente de Rayleigh de $x$ minimiza $f(\lambda) = ||(A - \lambda I)x||_2$.
```
Este cociente se utiliza para acelerar la convergencia en el método de la potencia o el de iteración inversa. En el caso de la iteración inversa se tiene el siguiente algoritmo.
Algoritmo: método de la potencia inversa con cociente de Rayleigh para matrices simétricas
Dados $A \in \mathbb{R}^{n \times n}$ simétrica, $x^{(0)} \in \mathbb{R}^n - {0}$ y $tol > 0$.
Repetir el siguiente bloque para $k=0,1,2,\dots$
Calcular $\mu^{(k)} = r(x^{(k)})$.
Resolver el sistema de ecuaciones lineales $(A-\mu^{(k)}I)z^{(k+1)} = x^{(k)}$ para $z^{(k+1)}$.
$x^{(k+1)} = \frac{z^{(k+1)}}{||z^{(k+1)}||_2}$
hasta convergencia: satisfacer criterio de paro en el que se utiliza $tol$.
````{admonition} Comentarios
Si $x$ es un vector que aproxima bien a un eigenvector de $A$ el cociente $r(x)$ es buena aproximación al eigenvalor asociado.
El algoritmo anterior tiene una tasa de convergencia cúbica si $x^{(0)}$ aproxima a un eigenvector de $A$.
El criterio de paro es de la forma: 
while ||Ax_k - mu_kx_k||_2 &gt; tol &amp;&amp; iterations &lt; max_iters
con tol una cantidad pequeña y positiva (comúnmente menor o igual a $10^{-8}$), iterations un contador de iteraciones.
````
Aplicaciones y usos: Page Ranking en un buscador de páginas web
En esta sección no se asume una estructura particular de la matriz $A \in \mathbb{R}^{n \times n}$.
```{sidebar} Un poco de historia...
En $1998$ S. Brin, L. Page, R. Motwani y T. Winograd publicaron el algoritmo Page Rank en el artículo The Page Rank Citation Ranking: Bringing Order to the Web. Tales autores crearon Google, a web search engine, que asigna a cada página una medida de importancia, quality measure, y a partir de ésta las ordena.
El $29$ de septiembre  de $2005$ una búsqueda en Google con la palabra university resultó en más de $2 \times 10^9$ de resultados.
```
En los inicios del internet muchos buscadores de páginas web eran deficientes para responder con resultados que fueran de interés a partir de una consulta que devolvía miles o millones de resultados. Hoy en día tenemos el buscador de páginas web Google que de forma exitosa filtra los resultados ordenando a las páginas web de acuerdo a una medida de importancia relativa overall para los millones de resultados que devuelve la consulta. 
Google utiliza el método de nombre PageRank para rankear cada página basándose en el grafo de la web. Tal ranking provee una quality measure que le nombra pagerank. A grandes rasgos, el método asigna a una página web un alto rank si es referenciada o tiene ligas que apuntan a ella, inlinks, por páginas que tienen un alto rank. La formulación matemática de lo anterior conduce a un problema de cálculo del eigenvalor con módulo máximo y su eigenvector asociado como se describe a continuación.
Descripción del método PageRank
El método resulta de la información del número de links que van hacia y desde una página Web. Para esto, se definen $O_i$ como el conjunto de outlinks o las páginas que son referenciadas por $i$ y el conjunto $I_i$ como el conjunto de inlinks o las páginas que hacen referencia a la página $i$. Esto se puede representar como:
```{margin}
En este grafo se tiene $|O_i| = 3$.
```
```{margin}
Para generar este grafo utilicé networkx.
```
<img src="https://dl.dropboxusercontent.com/s/sezz4empv27dtd4/inlinks_outlinks_pagerank.png?dl=0" heigth="400" width="400">
El pagerank de una página $i$ es una suma ponderada de aquellas páginas que hacen referencia a $i$. La ponderación se realiza por el número de outlinks que tienen tales páginas, esto es:
$$r_i = \displaystyle \sum_{j \in I_i} \frac{r_j}{N_j}$$
donde: $r_i$ es el pagerank de la página $i$ y $N_i = |O_i|$ es el el número de outlinks de la página $i$.
```{admonition} Comentarios
La forma que se utiliza aquí para calcular el pagerank de una página es una versión simplificada. Además de la suma ponderada se considera un parámetro para ajustar los pageranks nombrado damping factor y representa la distribución de las páginas a las que se pueden llegar desde una página (outlinks).
La ponderación anterior se realiza para hacer más "difícil" elevar el pagerank de una página $i$ pues podrían crearse muchas páginas "artificiales" que apunten a $i$ de forma sencilla. Además, la suma se realiza considerando los pageranks de las páginas $j$ por lo que páginas con pagerank bajo no tendrán una contribución alta en la suma. 
```
Para comprender cómo la definición anterior del pagerank de $i$ se relaciona con el problema de cálculo de eigenvalores y eigenvectores de una matriz, se define la matriz $Q$ de tamaño $n$ con entradas dadas por:
$$Q_{ij} = \begin{cases}
\frac{1}{N_j} & \text{ si existe un link de } j \text{ a } i \
0 & \text{ en otro caso}
\end{cases}
$$
Ejemplo
Considérese el siguiente grafo de un conjunto de páginas ${1, 2, 3, 4, 5, 6}$:
<img src="https://dl.dropboxusercontent.com/s/850i1rw3optocbt/graph_example_1_pagerank.png?dl=0" heigth="400" width="400">
Entonces la matriz asociada a este grafo es:
$$Q = \left [\begin{array}{cccccc}
0 & \frac{1}{3} & 0 & 0 & 0 & 0\
\frac{1}{3} & 0 & 0 & 0 & 0 & 0 \
0 & \frac{1}{3} & 0 & 0 & \frac{1}{3} & \frac{1}{2} \
\frac{1}{3} & 0 & 0 & 0 & \frac{1}{3} & 0 \
\frac{1}{3} & \frac{1}{3} & 0 & 0 & 0 & \frac{1}{2} \
0 & 0 & 1& 0 & \frac{1}{3} & 0
\end{array}
\right ]
$$
La matriz anterior podemos leerla como sigue:
(perspectiva de los renglones) La página $3$ es referenciada o tiene inlinks de las páginas ${2,5,6}$.
(perspectiva de las columnas) La página $3$ hace referencia o tiene outlink hacia la página ${6}$.
```{admonition} Observación
:class: tip
Obsérvese que la suma por columna es igual a $1$.
```
En el ejemplo anterior asúmase que inicialmente $r_2 = r_5 = r_6 = \frac{1}{6}$, entonces el pagerank de la página $3$ es: $r_3 = \frac{1}{6}\frac{1}{N_2} + \frac{1}{6}\frac{1}{N_5} + \frac{1}{6}\frac{1}{N_6} = \frac{1}{6}\frac{1}{3} + \frac{1}{6}\frac{1}{3} + \frac{1}{6}\frac{1}{2} = 0.19\bar{4}$
Continuando con la descripción del problema en el cálculo del pagerank, si se define el vector $r$ con los pageranks de cada página se tiene:
$$r = \left [\begin{array}{c}
r_1 \
r_2 \
\vdots \
r_n
\end{array}
\right ]
= \left [ \begin{array}{c}
\displaystyle \sum_{j \in I_1} \frac{r_j}{N_j} \
\displaystyle \sum_{j \in I_2} \frac{r_j}{N_j} \
\vdots \
\displaystyle \sum_{j \in I_n} \frac{r_j}{N_j} \
\end{array}
\right ]
$$
por lo que la suma $r_i = \displaystyle \sum_{j \in I_i} \frac{r_j}{N_j}$ se puede reescribir en una forma matricial como:
$$\lambda r = Qr, \quad \lambda = 1.$$
Esto es, $r$ es un eigenvector de $Q$ con eigenvalor $\lambda = 1$. Para calcular tal eigenvector se utiliza el {ref}método de la potencia&lt;MPOTMATSIM&gt;.
El modelo random surfer
En el método de pagerank se asume un modelo de caminata aleatoria en el grafo de la web, en el que una random surfer da clicks de forma sucesiva y aleatoria con igual probabilidad entre las páginas de los outlinks de cada página. Tal modelo induce una cadena de Markov. El eigenvector $r$ corresponde a la distribución de probabilidad estacionaria de la cadena. Los estados de la cadena son las páginas y las transiciones son los links entre las páginas.
```{admonition} Definición
Una matriz es estocástica por columnas si tiene elementos no negativos y la suma de cada columna es $1$. Ver stochastic matrix.
```
```{admonition} Definición
Un vector $x \in \mathbb{R}^n$ se nombra un vector de probabilidad si todas sus entradas son no negativas y suman $1$.
```
```{admonition} Definición
(por consistencia de esta nota se considerará a la matriz de transición como aquella que es estocástica por columnas)
Si $P$ es una matriz estocástica por columnas, $x_0$ es un vector de probabilidad o estado inicial, se define la siguiente sucesión de vectores de probabilidad:
$$x_{t+1} = Px_t, \quad t=0,1,2,\dots$$
A esta sucesión se le nombra cadena de Markov generada por $P$ y $x_0$.
```
```{admonition} Comentario
Una cadena de Markov es un proceso aleatorio en el que se describen posibles estados (eventos) que ocurren en puntos del tiempo y cada estado se encuentra completamente determinado (depende) del estado presente. A tales cadenas se les asocian matrices de transición cuyas entradas representan probabilidades de pasar de un estado a otro. En la literatura sobre cadenas de Markov se consideran matrices estocásticas por renglones, en este caso para el ejemplo anterior $Q^T$ es matriz de transición de la cadena de Markov.
```
¿Qué interpretación tiene la posición $i$ en el vector $r$?
El elemento en la posición $i$ de $r$, $r_i$, es la probabilidad que después de un número grande de clicks a páginas Web de forma aleatoria, la surfer llegue a la página $i$ en un tiempo determinado. El modelo por tanto asume que la surfer da clicks a ligas de forma aleatoria hasta llegar a su sitio destino.
Adición de un factor para obtener matrices estocásticas por columnas
```{margin}
Este supuesto es factible pues es poco probable que la surfer no pueda salir de una página por carecer de outlinks.
```
El cálculo del pagerank de $i$ se modifica asumiendo que siempre existe al menos un outlink de la página $i$ por un damping factor. Esto ayuda en la situación en el que una página no apunte a ninguna otra y corrige la no propagación de su pagerank.
```{margin}
Esta modificación en la literatura se conoce como teleportation.
```
En el ejemplo anterior este supuesto no se cumple pues la página $4$ no tiene outlinks (tiene una columna de ceros). Para modificar la situación anterior de forma matemática se realiza:
$$P = Q + \frac{1}{n}ed^T$$
donde: $e = (1, \dots, 1)^T$ y $d$ es un vector con entradas $d_j = \begin{cases} 1 & \text{ si } N_j = 0 \ 0 & \text{ en otro caso } \end{cases}$ para $j=1, \dots, n$.
```{admonition} Observación
:class: tip
La interpretación que tiene tal adición del término $\frac{1}{n}ed^T$ por llegar a una página que no tiene outlinks se realiza como sigue: la surfer necesariamente elegirá una página de forma aleatoria de acuerdo a una distribución de probabilidad. Una de las distribuciones que se eligen para modelar tal elección es la uniforme, esto es, elige cualquiera de las páginas con la misma probabilidad. 
```
```{margin}
Observa que si $P$ es una matriz estocástica por columna entonces $e^TP = e^T$ por lo que tiene eigenvalor $1$. Análogamente si $P$ es una matriz con entradas no negativas y $e^TP = e^T$ entonces $P$ es una matriz estocástica por columna.
```
```{admonition} Comentario
La matriz $P$ anterior es una matriz estocástica por columna. En la literatura sobre cadenas de Markov se consideran matrices estocásticas por renglones, en este caso para el ejemplo anterior $P^T$ es matriz de transición de la cadena de Markov.
```
Ejemplo
Para el ejemplo anterior se tiene:
End of explanation
print(1/n*np.outer(e,d))
Explanation: ```{margin}
Matriz que se le añadirá a $Q$, tiene en sus columnas el valor igual a $\frac{1}{6}$ para el ejemplo dado.
```
End of explanation
P = Q + 1/n*np.outer(e,d) #definition of P matrix
print(P)
Explanation: ```{margin}
Se obtiene $P$ a partir de $Q + \frac{1}{n}ed^T$.
```
End of explanation
k = 2
P_2 = np.linalg.matrix_power(P,k)
print(P_2)
print(np.sum(P_2, axis=0))
Explanation: ```{admonition} Comentario
Una propiedad que tienen las matrices estocásticas por columnas es $P^k$ es estocástica por columnas, $k \in \mathbb{N}$.
```
Ejemplo
End of explanation
k = 10
P_10 = np.linalg.matrix_power(P,k)
print(P_10)
np.sum(P_10, axis=0)
Explanation: ```{margin}
Ver sección stochastic matrices and Markov chains del libro Matrix Analysis and Applied Linear Algebra para caso general de matrices que convergen a un vector de distribución estacionario y que tienen mismos valores en sus entradas por renglones o columnas.
```
End of explanation
x_0 = np.array([.2, .2, .3, .2, .05, .05])
print(np.sum(x_0))
print(np.sum(P@x_0))
evalues, evectors = np.linalg.eig(P)
print(evalues)
evalues, evectors = np.linalg.eig(P_2)
print(evalues)
evalues, evectors = np.linalg.eig(P_10)
print(evalues)
print(np.sum(P_10@x_0))
Explanation: ```{admonition} Comentarios
Otras propiedades que tienen las matrices estocásticas por columnas son:
Dado un vector de probabilidad $x$ entonces $Px$ es un vector de probabilidad.
Si $\lambda$ es eigenvalor de $P$ entonces $|\lambda| \leq 1$.
```
Ejemplo
End of explanation
evalues, evectors = np.linalg.eig(P)
evector_1 = evectors[:,0]
print(evector_1)
print(P@evector_1)
Explanation: Una pregunta natural es si dada una matriz estocástica por columnas y un estado inicial $x_0$, los vectores de la cadena de Markov se estabilizan, esto es, la sucesión ${x_0, x_1, x_2, \dots }$ converge a algún vector $x$.
```{admonition} Definición
Dada una matriz estocástica por columna $P$ se le nombra $x$ un vector estacionario de P si y sólo si  $x$ es un vector de probabilidad con la propiedad:
$$Px = x$$
y $x$ es eigenvector de $P$ con eigenvalor $1$ asociado.
```
Ejemplo
En el ejemplo anterior.
End of explanation
l = [[0, 1/2, 1/2, 1/2, 0, 0],
     [1/2, 0, 1/2, 0, 0, 0],
     [1/2, 1/2, 0, 0, 0, 0],
     [0, 0, 0, 0, 0, 0],
     [0, 0, 0, 1/2, 0, 1],
     [0, 0, 0, 0, 1, 0]]
Q = np.array(l)
print(Q)
Explanation: Matriz irreducible y estocástica por columnas
Para asegurar que la matriz $P$ anterior tiene un único eigenvalor $\lambda_1$ igual a $1$ se utiliza la idea de teleportation y además se busca que la matriz sea irreducible.
```{admonition} Definición
Nombramos a una matriz irreducible como aquella para la que su grafo asociado siempre existe una ruta dirigida que une dos nodos distintos. Esto se nombra grafo fuertemente conectado.
```
```{margin}
En general la modificación es de la forma: $P = \alpha Q + (1-\alpha)ve^T$ con $v$ vector no negativo y $||v||_1 = 1$.
```
Una posible modificación es:
$$P = \alpha Q + (1-\alpha) \frac{1}{n}ee^T$$
para $\alpha \in [0,1]$. Se cumple que $P$ es estocástica por columnas, es irreducible y tiene entradas no negativas.
```{admonition} Comentarios
Es posible probar que si una matriz $P$ es estocástica por columnas e irreducible entonces $\lambda_1 = 1$ y existe un único eigenvector, $r$ asociado a $\lambda_1$ que tiene entradas positivas y $||r||_1 = 1$. Si además $P$ tiene entradas no negativas entonces el módulo de sus eigenvalores restantes son menores a $1$. En el contexto de cadenas de Markov tal vector $r$ es el vector estacionario.
La interpretación de la expresión anterior en la caminata aleatoria sobre el grafo de la web, es como sigue: en cada tiempo la surfer al visitar una página, si:
la página no tiene outlinks, elige una página aleatoria entre todas las páginas con probabilidad $\frac{1}{n}$.
la página tiene outlinks, elige un outlink con probabilidad $1 - \alpha$ y con probabilidad $\alpha$ elige cualquier página del grafo.
Típicamente $\alpha$ es $0.85$.
```
Ejemplo
Considérese el grafo
<img src="https://dl.dropboxusercontent.com/s/u2m8af7hfi83wq9/graph_example_2_pagerank.png?dl=0" heigth="400" width="400">
Entonces una surfer que haya entrado a la parte izquierda o derecha del grafo nunca saldrá. La matriz $Q$ asociada al grafo anterior es:
End of explanation
_,n = Q.shape
alpha = 0.85
P = alpha*Q + (1-alpha)*1/n*np.outer(e,e)
print(P)
Explanation: Tomando un valor de $\alpha = 0.85$ se tiene:
```{margin}
$P = \alpha Q + (1-\alpha) \frac{1}{n}ee^T$.
```
End of explanation
print(np.linalg.eigvals(Q))
eigvalues, eigvectors = np.linalg.eig(P)
print(eigvalues)
Explanation: Y sus eigenvalores son:
```{margin}
Compárese con los eigenvalores de $Q$.
```
End of explanation
print(-1*eigvectors)
Explanation: sus eigenvectores son:
```{margin}
Como se describió, el eigenvector asociado a $\lambda_1$ tiene entradas positivas.
```
End of explanation
x_0 = np.array([0.4, 0.2, 0.2, 0.1, 0.05, 0.05])
q_k = x_0.copy()
max_iter = 30
q_k_iter = np.zeros((n, max_iter))
lambda_k_iter = np.zeros(max_iter)
for k in range(max_iter):
    z_k = P@q_k
    q_k = z_k #changing a little the power method to have a probability vector
    q_k_iter[:,k] = q_k
    lambda_k = q_k.T@P@q_k/np.linalg.norm(z_k)**2
    lambda_k_iter[k] = lambda_k
print(lambda_k)
print(eigvalues)
Explanation: Verifiquemos que si definimos un vector de probabilidad $x_0$ con norma Euclidiana unitaria entonces al aplicar el método de la potencia
End of explanation
print(-1*eigvectors)
Explanation: Obsérvese que la aproximación es un múltiplo de la primer columna de la matriz siguiente
End of explanation
print(z_k)
division_to_get_constant = [(eigvectors[:,0][k]/z_k[k]).round(1) \
                            for k in range(len(z_k))]
print(division_to_get_constant)
cte = np.unique(division_to_get_constant)[0]
print(z_k*cte)
Explanation: Aproximación
End of explanation
np.sum(z_k)
Explanation: z_k es un vector de probabilidad:
End of explanation
n = 4
Q = np.array([[0,   1/3, 1, 1/2],
              [1/2, 0,   0, 0],
              [1/2, 1/3, 0, 1/2],
              [0,   1/3, 0, 0]])
Explanation: Ejemplo
Imaginemos tenemos $4$ páginas $web$ conectadas como se muestra a continuación:
<img src="https://dl.dropboxusercontent.com/s/ryb515n4y2273h9/graph_example_simple_pagerank.png?dl=0" heigth="700" width="700">
Entonces la matriz $Q$ asociada al grafo anterior es:
End of explanation
x_0 = 1/n*np.ones(n)
print(x_0)
Explanation: Definamos el vector de probabilidad inicial para aplicar el método de la potencia. Para ello, supongamos que la surfer visita aleatoriamente alguna de las páginas con la misma probabilidad, esto es:
End of explanation
e = np.ones(n)
_,n = Q.shape
alpha = 0.85
P = alpha*Q + (1-alpha)*1/n*np.outer(e,e)
print(P)
Explanation: $Q$ en una matriz irreducible y estocástica con la que podríamos trabajar, lo siguiente se realiza sólo para mostrar la aplicación de la transformación descrita.
End of explanation
q_k = x_0.copy()
max_iter = 30
q_k_iter = np.zeros((n, max_iter))
lambda_k_iter = np.zeros(max_iter)
for k in range(max_iter):
    z_k = P@q_k
    q_k = z_k #changing a little the power method to have a probability vector
    q_k_iter[:,k] = q_k
    lambda_k = q_k.T@P@q_k/(np.linalg.norm(z_k)**2)
    lambda_k_iter[k] = lambda_k
print(z_k)
print(lambda_k)
Explanation: Aplicamos el método de la potencia
End of explanation
np.sum(z_k)
Explanation: z_k es un vector de probabilidad:
End of explanation
evalues, evectors = np.linalg.eig(P)
print(evalues)
print(evectors)
Explanation: Verificamos que se aproxima al eigenvector asociado al eigenvalor igual a $1$
End of explanation
division_to_get_constant = [(evectors[:,0][k]/z_k[k]).round(1) \
                            for k in range(len(z_k))]
print(division_to_get_constant)
cte = np.unique(division_to_get_constant)[0]
print(z_k*cte)
Explanation: Obsérvese que la aproximación es un múltiplo de la primer columna de la matriz anterior:
End of explanation
print(z_k)
Explanation: Interpretación
End of explanation
<END_TASK>