<jupyter_start><jupyter_text>Tema 3: Aproximación de autovalores<jupyter_code>!pip install -r https://raw.githubusercontent.com/alexmascension/ANMI/main/requirements.txt
from sympy import *
from sympy.matrices import Matrix as mat
from sympy.matrices import randMatrix
from sympy import symbols
import sympy

import numpy as np

from scipy.linalg import orth
from anmi.genericas import norma, print_verbose

from anmi.T2 import factorizacion_QR
from anmi.T3 import matriz_krylov, sucesion_krylov, potencia_iterada, metodo_autovals_QR<jupyter_output><empty_output><jupyter_text>Sucesiones de KrylovSea $A$ una matriz (aplicación lineal) y $x$ un vector. Si aplicamos la multiplicación de $A$ por $x$ de manera iterada obtenemos una serie de vectores $\{x, Ax, A^2x, A^3x, \cdots\}$. Si $x$ no es un autovector de $A$, entonces esa sucesión tendrá $n$ (dimensión de $A$) vectores independientes. Si $x$ es un autovector, con su autovalor $\lambda$, entonces la sucesión de vectores será, $\{x, \lambda x, \lambda^2x, \lambda^3x, \cdots\}$. Estas sucesiones de vectores se llaman *sucesiones de Krylov*.Por otra parte, por el teorema de Cayley-Hamilton se tiene que $A^nx$ tiene que ser una combinación lineal de los siguientes elementos de la sucesión, es decir:$$(-1)^nA^n + a_{n-1}A^{n-1} + \cdots + a_1A + a_0I = 0$$Luego si tomamos $a = \begin{bmatrix}a_0\\a_1\\ \cdots \\ a_n\end{bmatrix}$ se tiene que$$(x|Ax|\cdots|A^{n-1}x)a = (-1)^{n+1}A^nx$$Y si resolvemos $a$, entonces se tienen los coeficientes del polinómio característico $p(\lambda) = a_0 + a_1\lambda + a_2\lambda^2 + \cdots + a_n\lambda^n$$<jupyter_code>help(matriz_krylov)
# EJERCICIO 26
A = mat([[1, 1, 1], [0, 2, 2], [3, -1, 0]])
x = mat([[1, 0, 0]]).T

m_krylov = matriz_krylov(A, x)
m_krylov
sk, a = sucesion_krylov(A, x)
sk
# EJEMPLO 15
A = mat([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
x = mat([[-1, 0, 1]]).T

matriz_krylov(A, x)<jupyter_output><empty_output><jupyter_text>Método de la potencia iteradaEn el método de la potencia iterada, se aplica la matriz de krylov hasta una potencia determinada, $k$. Entonces, se tiene que $$\lim_{k \to \infty} \frac{A^kw}{A^{k-1}w} = |\lambda_1|$$Es decir, el mayor autovalor.Además, si tomamos $ B= A^{-1} $, tenemos que$$\lim_{k \to \infty} \frac{B^kw}{B^{k-1}w} = \frac{1}{|\lambda_n|}$$Donde $\lambda_n$ es el menor autovalor.<jupyter_code>help(potencia_iterada)
A.eigenvals()
x = mat([[-2, 0, 1]]).T

matriz_krylov(A, x, 17)
x = mat([[-1, 0, 0]]).T

np.array(potencia_iterada(A, x, 30, devolver_ultimo=False)[:, -3:], dtype=float)
x = mat([[-1, 0, 0]]).T

np.array(potencia_iterada(A, x, 300, devolver_ultimo=True), dtype=float)
N(2+sqrt(2))  # Vemos que converge al mayor autovalor
np.array(potencia_iterada(A**-1, x, 300, devolver_ultimo=True), dtype=float)
1/N(2-sqrt(2))  # Y lo mismo con el menor
# Si tomamos una matriz ortogonal, el metodo de la potencia no tiene validez 
# porque se requiere que haya autovalores dominantes, y en este caso los 
# autovalores tienen módulo 1.

dict_QR = factorizacion_QR(A)
Q = dict_QR['Q']

Q
Q * Q.T
Q.eigenvals()
N(-1/2 + 3*sqrt(70)/70 + 3*sqrt(14)/28 + sqrt(5)/5 + sqrt(70)*I*sqrt(6*sqrt(14) + 20*sqrt(5) + 73)/140)
matriz_krylov(Q, x, 5)
N(matriz_krylov(N(Q), x, 30), 4)
potencia_iterada(N(Q), x, 100, devolver_ultimo=False)[:, -5:]  # No hay una convergencia<jupyter_output><empty_output><jupyter_text>Método QREl método QR consiste en aplicar los siguientes pasos:$$A^{(1)} = A$$De ahí sacamos que $$A^{(1)}  = Q^{(1)}R^{(1)}$$De ahí construimos:$$A^{(2)}  = R^{(1)}Q^{(1)}$$Y se cumple que $A^{(1)}$ y $A^{(2)}$ son semejantes, luego tienen los mismos autovalores.Con ello se reitera el proceso, y se cumple que las matrices equivalentes construidas convergen a una matriz triangular superior. Los la diagonal de $A^{(k)}$ converge a los autovalores de $A$.<jupyter_code>help(metodo_autovals_QR)
dict_QR = metodo_autovals_QR(A, n_iters=10)
N(dict_QR['A'][-2], 3)
N(dict_QR['A'][-1], 3)
N(2- sqrt(2), 3), 2, N(2 + sqrt(2), 3), 
A = mat([[1, 1, 1], [0, 0, 1], [0, 1, 1]])
dict_QR = metodo_autovals_QR(A, n_iters=30, verbose=False)
N(dict_QR['A'][-15], 3)
N(dict_QR['A'][-1], 3)
[N(i, 3) for i in list(A.eigenvals().keys())]<jupyter_output><empty_output><jupyter_text>Ejercicios  Ejercicio 27Determinar las primera iteraciones generadas por el método de la potencia con normalización con norma infinito cuando se aplica a la matriz con vector inicial.<jupyter_code>A = Matrix([[0, -1, 1], [0, 1, -1], [-1, -1, 2]])
A
x0 = Matrix([1, -1, 2])
n_iters = 15
m_krylov_x, m_krylov_w = zeros(A.shape[0], n_iters), zeros(A.shape[0], n_iters)
m_krylov_x[:, 0] = x0
m_krylov_w[:, 0] = x0/max(x0)

for i in range(1, n_iters):
    kriv_i_x = A * m_krylov_w[:, i - 1]
    m_krylov_x[:, i] = kriv_i_x 
    m_krylov_w[:, i] = kriv_i_x / max(kriv_i_x) 

m_krylov_x
m_krylov_w
[N(i) for i in (np.array(m_krylov_x[0, 1:]) / np.array(m_krylov_w[0, :-1])).tolist()[0]]
# Probamos sin normalizar

n_iters = 15
m_krylov_x  = zeros(A.shape[0], n_iters)
m_krylov_x[:, 0] = x0

for i in range(1, n_iters):
    kriv_i_x = A * m_krylov_x[:, i - 1]
    m_krylov_x[:, i] = kriv_i_x 

m_krylov_x 
[N(i) for i in (np.array(m_krylov_x[0, 1:]) / np.array(m_krylov_x[0, :-1])).tolist()[0]]<jupyter_output><empty_output><jupyter_text>Ejercicio 28Realizar una iteración con el método QR para el cálculo de autovalores de la matriz $A$<jupyter_code>A = Matrix([[1, 1, 0], [2, 1, 0], [2, 0, 1]])
metodo_autovals_QR(A, n_iters=1, verbose=True)
dict_QR = metodo_autovals_QR(A, n_iters=15, verbose=False)
N(dict_QR['A'][-1])
N(dict_QR['A'][-2])
np.sort(solve(det(A - eye(3) * Symbol('lambda')), Symbol('lambda')))[::-1]<jupyter_output><empty_output><jupyter_text>Ejercicio 29Sea $A$ con $\theta$ dado. Determinar el n-ésimo términa de la sucesión de Krylov asociada a $A$ y al vector $x$. Por medio de la sucesión de Krylov determia el polinomio característico de $A$.<jupyter_code>t = Symbol('theta')
A = Matrix([[1, 0, 0], [0, cos(t), -sin(t)], [0, sin(t), cos(t)]])
x0 = Matrix([1, cos(t), -sin(t)])
ss = simplify(matriz_krylov(A, x0, n_iters=8))
ss<jupyter_output><empty_output><jupyter_text>Así a priori tiene pintas de que $A^nx = (1, \cos(n-1)\theta, sin(n-1)\theta)^t$. Como $n=3$ solo necesitamos los primeros 4 términos de la sucesión de krylov para resolver el sistema y crear el polinomio característico.<jupyter_code>M = ss[:3, :3]
am = Matrix([Symbol('a0'), Symbol('a1'), Symbol('a2')])
rhs = ss[:3, 3]
sol = simplify(M.inv() * rhs)
sol
l = Symbol('lambda')
simplify(sol[0] + sol[1] * l + sol[2] * l ** 2 - l ** 3)<jupyter_output><empty_output><jupyter_text>Ejercicio 31Si $\alpha$ es una constante, aproximar el autovalore de módulo máximo de la matriz $A$ usando el método de la potencia (sin normalización) como vector inicial $v$, probando las diferentes componentes.<jupyter_code>a = Symbol('alpha')
n = Symbol('n')
A = Matrix([[1, 0, 0], [0, -2*cos(a)**2, 2*cos(a)*sin(a)], [0, 2*cos(a)*sin(a), -2*sin(a)**2]])
A
v = Matrix([1, 1, 1])
v
pol, vecfinal = sucesion_krylov(A, v)
pol
solve(pol, l)
vecfinal
# Vamos a ir desglosando 
mk = matriz_krylov(A, v, n_iters=7)
mk
simplify(Matrix(np.array(mk[:, 1:]) / np.array(mk[:, :-1])))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17597222222222222, 'fraction_numerical': 0.03666666666666667, 'mean_word_length': 5.338908450704225, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '427543', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3304, 'n_tokens_neox': 3150, 'n_words': 818, 'path': 'notebook/T3.ipynb', 'repo_name': 'alexmascension/ANMI', 'size': '7200'}