<jupyter_start><jupyter_text>GMRES Algunos Ejemplos<jupyter_code>import numpy as np
import scipy.linalg as spla
from scipy.sparse.linalg import gmres, LinearOperator
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Error Para analizar convergencia.<jupyter_code>error = lambda x, x_a: np.linalg.norm(x - x_a, ord=2, axis=1) / np.linalg.norm(x)
back_error = lambda x, A, b: np.linalg.norm(b - np.dot(A, x), ord=2) / np.linalg.norm(b)<jupyter_output><empty_output><jupyter_text>Además, para graficar la convergencia.<jupyter_code>def plot(err, b_e):
    Ni_e = np.arange(err.shape[-1])
    Ni_b = np.arange(b_e.shape[-1])
    plt.figure(figsize=(12, 6))
    plt.plot(Ni_e, err, 'bo', label=r'$||\mathbf{x} - \mathbf{x}_a||_2/||\mathbf{x}||_2$')
    plt.plot(Ni_b, b_e, 'rx', label=r'$||\mathbf{b} - A\,\mathbf{x}||_2/||\mathbf{b}||_2$')
    plt.yscale('log')
    plt.xlabel('# Iterations')
    plt.ylabel('Error')
    plt.legend()
    plt.grid(True)
    plt.show()<jupyter_output><empty_output><jupyter_text>Implementación $QR$ reducidoEl algoritmo de Gram-Schmidt para obtener $\hat{Q}$ y $\hat{R}$.<jupyter_code>def QR(A, modified=True):
    m, n = A.shape
    Q = np.zeros((m, n))
    R = np.zeros((n, n))
    # Gram-Schmidt orthogonalization
    for j in range(n):
        y = A[:, j]
        if modified:
            for i in range(j):
                R[i, j] = np.dot(Q[:, i], y)
                y = y - R[i, j] * Q[:, i]
        else:
            for i in range(j):
                R[i, j] = np.dot(Q[:, i], A[:, j])
                y = y - R[i, j] * Q[:, i]
        R[j, j] = np.linalg.norm(y)
        Q[:, j] = y / R[j, j]
    return Q, R<jupyter_output><empty_output><jupyter_text>Mínimos Cuadrados- Utilizando $QR$    1. Obtener $QR$ reducido, es decir, $A=\hat{Q}\,\hat{R}$    2. Resolver $\hat{R}\,\hat{\mathbf{x}} = \hat{Q}^*\mathbf{b}$        - Utilizando Ecuaciones Normales    1. Resolver $A^*A\hat{\mathbf{x}} = A^*\mathbf{b}$<jupyter_code>def lstsq(A, b, qr=True):
    if qr:
        Q, R = QR(A) # QR reduced
        x = spla.solve_triangular(R, np.dot(Q.T, b)) # Solve Rx=Q^*b
    else:
        x = np.linalg.solve(np.dot(A.T, A), np.dot(A.T, b)) # Solve A^*Ax=A^*b
    return x<jupyter_output><empty_output><jupyter_text>Implementación GMRESVersión modificada de [1].<jupyter_code>def GMRES(A, b, x_0=None, tol=1e-10, max_iter=100, show=False, np_ls=False):
    m = b.shape[0]
    
    M = np.min([m, max_iter])
    
    X = np.zeros((M + 1, m))
    if x_0 is not None:
        X[0] = x_0
    
    r_0 = b - np.dot(A, X[0])
    n_r0 = np.linalg.norm(r_0)
    
    Q = np.zeros((m, m))
    H = np.zeros((m, m))
    
    Q[:, 0] = r_0 / n_r0

    for k in range(M):
        
        # Arnoldi Iteration
        y = np.dot(A, Q[:, k])
        for j in np.arange(k+1):
            H[j, k] = np.dot(Q[:,j], y)
            y = y - H[j, k] * Q[:,j]
        
        # Build H matrix and fill Q_k
        if k + 1 < m: 
            H[k+1, k] = np.linalg.norm(y)
            
            if np.abs(H[k+1, k]) > 1e-16:
                Q[:, k+1] = y / H[k+1, k]  

            e_1 = np.zeros((k + 1) + 1)        
            e_1[0] = 1
            H_tilde = H[:(k + 1) + 1, :k + 1]
            s = k + 1
        else: # When we use all Krylov subspace
            H_tilde = H[:k + 1, :k + 1]
            s = m + 1
            
        # Solving the 'SMALL' least squares problem. 
        if np_ls: 
            c_k, _, _, _ = np.linalg.lstsq(H_tilde, n_r0 * e_1, rcond=None) # Using Numpy Implementation (faster)
        else:
            c_k = lstsq(H_tilde, n_r0 * e_1) # Using algorithms seen in classes (slower)
        
        
        if show:
            print("k:", k)
            print("H:\n", H)
            print("H_t:\n", H_tilde)
            print("e_1:\n", e_1)
            print("c_k:\n", c_k)
            print("Q_k:\n", Q[:,:s])
            print()
            
        X[k+1] = X[0] + np.dot(Q[:, :s], c_k)  
        
        if np.linalg.norm(np.dot(A, X[k+1]) - b) < tol:
            X = X[:k + 2] 
            break
        
    return X<jupyter_output><empty_output><jupyter_text>Ejemplo 1Matriz $A$ aleatoria.<jupyter_code>np.random.seed(1234567890)
m = 50
A_1 = np.random.rand(m, m)
b_1 = np.dot(A_1, np.ones(m))<jupyter_output><empty_output><jupyter_text>Solver de ```numpy```<jupyter_code>x_1 = np.linalg.solve(A_1, b_1)
x_1<jupyter_output><empty_output><jupyter_text>Implementación<jupyter_code>x_g_1 = GMRES(A_1, b_1, tol=1e-10)
x_g_1[-1]<jupyter_output><empty_output><jupyter_text>GMRES ```scipy```<jupyter_code>x_s_1, _ = gmres(A_1, b_1)
x_s_1<jupyter_output><empty_output><jupyter_text>Convergencia<jupyter_code>b_e_1 = np.array([back_error(x_i, A_1, b_1) for x_i in x_g_1])
err_1 = error(x_1, x_g_1)
plot(err_1, b_e_1)<jupyter_output><empty_output><jupyter_text>¿Podemos Mejorar la convergencia? Intentemos utilizando *GMRES* con reinicio.<jupyter_code>p = 20
M = 5
X_r = np.zeros(((M * p), m))
x_0_ = None
for r in range(M):
    x_g_1_a = GMRES(A_1, b_1, x_0=x_0_, max_iter=p, tol=1e-16)
    X_r[r * p: (r + 1) * p] = x_g_1_a[1:]
    x_0_ = x_g_1_a[-1]
b_e_1_a = np.array([back_error(x_i, A_1, b_1) for x_i in X_r])
err_1_a = error(x_1, X_r)
plot(err_1_a, b_e_1_a)<jupyter_output><empty_output><jupyter_text>¿Problema? Quizás la matriz $A$ es mal condicionada.<jupyter_code>np.linalg.cond(A_1)<jupyter_output><empty_output><jupyter_text>Ejemplo 2Modifiquemos levemente nuestra matriz $A$<jupyter_code>A_2 = A_1 + m * np.eye(m) 
b_2 = np.dot(A_2, np.ones(m))<jupyter_output><empty_output><jupyter_text>Solver de ```numpy```<jupyter_code>x_2 = np.linalg.solve(A_2, b_2)
x_2<jupyter_output><empty_output><jupyter_text>Implementación<jupyter_code>x_g_2 = GMRES(A_2, b_2, tol=1e-10)
x_g_2[-1]<jupyter_output><empty_output><jupyter_text>GMRES ```scipy```<jupyter_code>x_s_2, _ = gmres(A_2, b_2)
x_s_2
b_e_2 = np.array([back_error(x_i, A_2, b_2) for x_i in x_g_2])
err_2 = error(x_2, x_g_2)
plot(b_e_2, err_2)<jupyter_output><empty_output><jupyter_text>Ejemplo 3: Evitar construir $A$ explícitamenteSupongamos necesitamos resolver el sistema $A\mathbf{x}=\mathbf{b}$, con $A\in\mathbb{R}^{n\times n}$ tridiagonal, y $\mathbf{b}\in\mathbb{R}^n$. ¿Necesitamos calcular el producto matriz vector $A\mathbf{v}$ para $\mathbf{v}\in\mathbb{R}^n$ cualquiera?La estructura de $A$ tridiagonal es:\begin{equation}    A=    \begin{bmatrix}        a_{1} & b_{1}\\        c_{1}& a_{2} & b_{2}\\        &c_{2}&\ddots &\ddots \\        & & \ddots & \ddots &b_{n-1}\\        & & & c_{n-1} & a_{n}    \end{bmatrix}\end{equation}  Construcción $A$ tridiagonal<jupyter_code>def buildA(n):
    np.random.seed(1234)
    # Random arrays to fill diagonals
    a_ = np.random.rand(n)
    b_ = np.random.rand(n-1)
    c_ = np.random.rand(n-1)
    # Temporal matrices to add diagonals
    a_tmp = np.zeros((n, n))
    b_tmp = np.zeros_like(a_tmp)
    c_tmp = np.zeros_like(a_tmp)
    # Diagonals
    a_tmp = a_ * np.eye(n) # Main diagonal
    b_tmp[:-1,1:] = b_ * np.eye(n-1) # Above diagonal
    c_tmp[1:,:-1] = c_ * np.eye(n-1) # Below diagonal 

    # Build A
    A = a_tmp + b_tmp + c_tmp
    
    del a_tmp, b_tmp, c_tmp
    
    return A
n = 5
A = buildA(n)
# Just show when A is small...
if n < 7:
    print(A)<jupyter_output>[[0.19151945 0.27259261 0.         0.         0.        ]
 [0.87593263 0.62210877 0.27646426 0.         0.        ]
 [0.         0.35781727 0.43772774 0.80187218 0.        ]
 [0.         0.         0.50099513 0.78535858 0.95813935]
 [0.         0.         0.         0.68346294 0.77997581]]<jupyter_text>Construcción ```afun(v)```Con esta función queremos calcular el producto matriz vector $A\mathbf{v}$ sin necesidad de ocupar todos los coeficientes de $A$. Dado que la matriz tiene muchas entradas nulas evitamos aplicar el algoritmo completo de multiplicación matriz-vector.<jupyter_code>def buildAfun(A, v):
    # Get matrix diagonals
    a = np.diag(A)
    b = np.diag(A, 1)
    c = np.diag(A, -1)
    
    # Vector size
    n = v.shape[0]
    u = np.zeros(n) # Output vector
    
    u[0] = a[0] * v[0] + b[0] * v[1] # First element of u
    u[-1] = c[-1] * v[-2] + a[-1] * v[-1] # Last element of u
    for i in range(1, n-1):
        u[i] = c[i-1] * v[i-1] + a[i] * v[i] + b[i] * v[i+1] # Inside u...
    return u<jupyter_output><empty_output><jupyter_text>Probemos el tiempo de cómputo multiplicando por un vector aleatorio.<jupyter_code>n_test = 10000 # For large n we can see the difference in computational time
A_test = buildA(n_test)
v = np.random.rand(n_test)
afun_test = lambda v: buildAfun(A_test, v)
%timeit np.dot(A_test, v)
%timeit afun_test(v)
Av_1 = np.dot(A_test, v)
Av_2 = afun_test(v)
print(np.linalg.norm(Av_1- Av_2))<jupyter_output>4.8266229755959866e-15<jupyter_text>Notamos que calcular el producto punto es mucho más costoso que utilizar nuestra función. Luego vemos que el resultado de ambos $A\mathbf{v}$ es similar. Importante que la diferencia de tiempo computacional se aprecia para valores de $n$ grandes, ya que ```numpy``` ya tiene considerablemente optimizado el producto matriz-vector.  "Matrix-free" GMRESModificamos nuestro algoritmo para que no calcule el producto matriz-vector y utilizamos nuestra ```afun(v)```.<jupyter_code>def GMRES2(fun, b, x_0=None, tol=1e-10, max_iter=100, show=False, np_ls=False):
    m = b.shape[0]
    M = np.min([m, max_iter])
    X = np.zeros((M + 1, m))
    if x_0 is not None:
        X[0] = x_0
    
    r_0 = b - fun(X[0])
    n_r0 = np.linalg.norm(r_0)
    Q = np.zeros((m, m))
    H = np.zeros((m, m))
    Q[:, 0] = r_0 / n_r0

    for k in range(M):
        
        # Arnoldi Iteration
        y = fun(Q[:, k])
        for j in np.arange(k+1):
            H[j, k] = np.dot(Q[:,j], y)
            y = y - H[j, k] * Q[:,j]
        
        # Build H matrix and fill Q_k
        if k + 1 < m: 
            H[k+1, k] = np.linalg.norm(y)
            
            if np.abs(H[k+1, k]) > 1e-16:
                Q[:, k+1] = y / H[k+1, k]  

            e_1 = np.zeros((k + 1) + 1)        
            e_1[0] = 1
            H_tilde = H[:(k + 1) + 1, :k + 1]
            s = k + 1
        else: # When we use all Krylov subspace
            H_tilde = H[:k + 1, :k + 1]
            s = m + 1
            
        # Solving the 'SMALL' least squares problem. 
        if np_ls: 
            c_k, _, _, _ = np.linalg.lstsq(H_tilde, n_r0 * e_1, rcond=None) # Using Numpy Implementation
        else:
            c_k = lstsq(H_tilde, n_r0 * e_1) # Using algorithms seen in classes
        
        
        if show:
            print("k:", k)
            print("H:\n", H)
            print("H_t:\n", H_tilde)
            print("e_1:\n", e_1)
            print("c_k:\n", c_k)
            print("Q_k:\n", Q[:,:s])
            print()
            
        X[k+1] = X[0] + np.dot(Q[:, :s], c_k)  
        
        if np.linalg.norm(fun(X[k+1]) - b) < tol:
            X = X[:k + 2] 
            break
        
    return X
n = 100
A_3 = buildA(n)
b_3 = np.dot(A_3, np.arange(1, n+1))
afun = lambda v: buildAfun(A_3, v)
x_3 = GMRES(A_3, b_3, np_ls=True, max_iter=n)
x_3[-1]
x_3 = GMRES2(afun, b_3, np_ls=True, max_iter=n)
x_3[-1]<jupyter_output><empty_output><jupyter_text>Notar que en ambos casos obtenemos la solución, pero el impacto de ```afun(v)``` se ve cuando tenemos sistemas de ecuaciones muy grandes.  Comentarios Recordar que en cada iteración la complejidad de nuestro algoritmo es la complejidad del producto matriz-vector $O(m^2)$ más la resolución del problema de mínimos cuadrados $O(k^3)$ si utilizamos $QR$ por ejemplo. Se omite la obtención de $h_{j,k}$ ya que es un producto punto, es decir, $O(m)$. De acuerdo a [3] podemos reducir la complejidad del problema de mínimos cuadrados aprovechando la estructura de $\tilde{H}_k$ de manera que la mayor parte del costo computacional se encuentra en resolver $\hat{R}\,\mathbf{x}=\hat{Q}^*\mathbf{b}$, que sería $O(k^2)$. Ahora, nuestro algoritmo tiene como complejidad $k\left(O(m^2) + O(k^2)\right)$, donde el peor caso es $k=m$, o sea obtenemos una complejidad de $O(m^3)$ similar a $PA=LU$ por ejemplo. El objetivo es resolver el problema cuando $k<m$, es decir, con complejidad $O(k^3)$. Notar que el uso de ```afun(v)``` nos ayuda a evitar el costo de la multiplicación matriz-vector. ---  Ejercicio 4: Ecuación de SylvesterLa ecuación de Sylvester es una ecuación matricial de la forma:\begin{equation}    A\,X + X\,B = C,\end{equation}donde $A, B, C \in \mathbb{R}^{n \times n}$ son matrices conocidas y $X\in \mathbb{R}^{n \times n}$ nuestra matriz incógnita. Utilicemos un método iterativo para resolverlo. Algoritmo obtenido de [4]. Idea 1\begin{equation}    \begin{split}        A\,X + X\,B & = C \\        A\,X & = C - X\,B     \end{split}\end{equation}Luego iteramos $A\,X_{i+1} = C - X_i\,B$, con $X_0$ una estimación inicial. Idea 2\begin{equation}    \begin{split}        A\,X + X\,B & = C \\        X\,B & = C - A\,X \\        (X\,B)^T & = (C - A\,X)^T  \\        B^TX^T &= (C - A\,X)^T     \end{split}\end{equation}Luego iteramos $B^T\,X^T_{i+1} = (C - A\,X_i)^T$, con $X_0$ una estimación inicial.<jupyter_code>np.random.seed(2468)
n = 10
A_4 = np.random.rand(n, n) + n * np.eye(n)
B_4 = np.random.rand(n, n)
C_4 = np.random.rand(n, n)
def iterativeSylvester(A, B, C, m, alg=1, tol=1e-10):
    X0 = np.zeros_like(A)
    X1 = np.zeros_like(A)
    if alg == 1: # Algorithm 1
        for i in range(m):
            X1 = np.linalg.solve(A, C - np.dot(X0, B))
            X0 = X1
            if np.linalg.norm(np.dot(A, X1) + np.dot(X1, B) - C) < tol:
                break
        return X1
    elif alg == 2: # Algorithm 2
        for i in range(m):
            X1 = np.linalg.solve(B.T, (C - np.dot(A, X0)).T).T
            X0 = X1
            if np.linalg.norm(np.dot(A, X1) + np.dot(X1, B) - C) < tol:
                break
        return X1<jupyter_output><empty_output><jupyter_text>Propuesta 1<jupyter_code>X_1 = iterativeSylvester(A_4, B_4, C_4, 70)
np.linalg.norm(np.dot(A_4, X_1) + np.dot(X_1, B_4) - C_4)<jupyter_output><empty_output><jupyter_text>Propuesta 2<jupyter_code>X_2 = iterativeSylvester(A_4, B_4, C_4, 10, 2)
np.linalg.norm(np.dot(A_4, X_2) + np.dot(X_2, B_4) - C_4)<jupyter_output><empty_output><jupyter_text>Propuesta 3Utilizando GMRES. En este enfoque, vectorizamos la matriz $C$ para resolver un problema $A\,\mathbf{x}=\mathbf{b}$ donde $\mathbf{b}=\operatorname{vec}(C)$. Además, nuestro "producto matriz vector", será $AX+XB$.<jupyter_code>def AXBX(x, A, B):
    X = np.reshape(x, A.shape) # Get flattened version of x and reshape to get X
    out = np.dot(A, X) + np.dot(X, B) # Out "matrix-vector" multiplication
    return out.flatten()

Ax = lambda x: AXBX(x, A_4, B_4)<jupyter_output><empty_output><jupyter_text>Resolviendo con GMRES implementado<jupyter_code>x_u = GMRES2(Ax, C_4.flatten(), np_ls=True, max_iter=n ** 2) # C.flatten() is vec(C)...
X_3_u = np.reshape(x_u[-1], A_4.shape)
np.linalg.norm(np.dot(A_4, X_3_u) + np.dot(X_3_u, B_4) - C_4)<jupyter_output><empty_output><jupyter_text>Resolviendo con GMRES de ```scipy```<jupyter_code># This is required by Scipy GMRES... It's just a funcion to compute matrix-vector multiplication
afun2 = LinearOperator((n ** 2, n ** 2), matvec = Ax) 
x_s, exitCode = gmres(afun2, C_4.flatten(), tol=1e-10)
X_3_s = np.reshape(x_s, A_4.shape)
np.linalg.norm(np.dot(A_4, X_3_s) + np.dot(X_3_s, B_4) - C_4)<jupyter_output><empty_output><jupyter_text>Solución del problema<jupyter_code>X = spla.solve_sylvester(A_4, B_4, C_4)
print("Propuesta 1: ", np.linalg.norm(X_1 - X))
print("Propuesta 2: ", np.linalg.norm(X_2 - X)) 
print("Propuesta 3a:", np.linalg.norm(X_3_s - X))
print("Propuesta 3b:", np.linalg.norm(X_3_u - X))<jupyter_output>Propuesta 1:  4.689803856203944e-12
Propuesta 2:  1.690995795460254e+21
Propuesta 3a: 1.6368741557609354e-11
Propuesta 3b: 2.5798177172434354e-12
{'fraction_non_alphanumeric': 0.16093257718966603, 'fraction_numerical': 0.03938248267170762, 'mean_word_length': 2.780609814197237, 'pattern_counts': {'":': 0, '<': 118, '<?xml version=': 0, '>': 112, 'https://': 0, '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': '28', 'dir': 'jupyter-structured-clean-dedup', 'id': '390670', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7165, 'n_tokens_neox': 6865, 'n_words': 1512, 'path': 'material/07_gmres/gmres.ipynb', 'repo_name': 'DoVirgo/INF-285', 'size': '15870'}