1. Conceptos básicos

Un autovalor \(\lambda\) de una matriz \(A\) satisface la ecuación \(A\,v = \lambda\,v\), donde \(v\) es el correspondiente autovector. En problemas de gran escala (matrices de varios miles o millones de filas) los métodos directos (descomposición de Schur, QR) son prohibitivamente costosos, de ahí la necesidad de métodos iterativos que convergen a los valores más relevantes sin requerir la matriz completa en memoria.

2. Principales métodos iterativos

v_{k+1} = A v_k / \|A v_k\|

\mu_k = (v_k^T A v_k) / (v_k^T v_k)

3. Implementación práctica en Python

Usaremos numpy y scipy.sparse para operar con matrices densas y esparcidas. Cada algoritmo se presenta como una función reutilizable.

3.1 Power Method

import numpy as np
def power_method(A, num_iter=1000, tol=1e-10):
    n = A.shape[0]
    b_k = np.random.rand(n)
    b_k /= np.linalg.norm(b_k)
    lambda_old = 0.0
    for _ in range(num_iter):
        b_k1 = A @ b_k
        b_k1_norm = np.linalg.norm(b_k1)
        b_k = b_k1 / b_k1_norm
        lambda_new = b_k.T @ A @ b_k
        if np.abs(lambda_new - lambda_old) < tol:
            break
        lambda_old = lambda_new
    return lambda_new, b_k

Ejemplo con una matriz aleatoria:

np.random.seed(0)
A = np.random.rand(5,5)
val, vec = power_method(A)
print(f"Autovalor dominante: {val:.6f}")

3.2 Inverse Power con shift

import scipy.sparse.linalg as spla
def inverse_power_method(A, shift=0.0, num_iter=100, tol=1e-10):
    n = A.shape[0]
    I = np.eye(n)
    M = A - shift * I
    # Factorizamos una sola vez (LU) para acelerar iteraciones
    lu = spla.splu(spla.csc_matrix(M))
    b_k = np.random.rand(n)
    b_k /= np.linalg.norm(b_k)
    lambda_old = 0.0
    for _ in range(num_iter):
        # Resolver M y_k = b_k
        y = lu.solve(b_k)
        b_k = y / np.linalg.norm(y)
        lambda_new = b_k.T @ A @ b_k
        if np.abs(lambda_new - lambda_old) < tol:
            break
        lambda_old = lambda_new
    return lambda_new, b_k

Buscar el autovalor más cercano a 2.5:

val, vec = inverse_power_method(A, shift=2.5)
print(f"Autovalor cercano a 2.5: {val:.6f}")

3.3 Rayleigh Quotient Iteration (RQI)

def rqi(A, x0=None, max_iter=30, tol=1e-12):
    n = A.shape[0]
    if x0 is None:
        x = np.random.rand(n)
    else:
        x = x0
    x = x / np.linalg.norm(x)
    mu = x.T @ A @ x
    for _ in range(max_iter):
        # (A - mu I) y = x
        try:
            y = np.linalg.solve(A - mu * np.eye(n), x)
        except np.linalg.LinAlgError:
            # fallback a solver iterativo
            y = spla.cg(A - mu * np.eye(n), x)[0]
        x = y / np.linalg.norm(y)
        mu_new = x.T @ A @ x
        if np.abs(mu_new - mu) < tol:
            break
        mu = mu_new
    return mu, x

Ejemplo con una matriz simétrica:

A_sym = np.array([[4,1,0],[1,3,1],[0,1,2]], dtype=float)
val, vec = rqi(A_sym)
print(f"Autovalor con RQI: {val:.12f}")

3.4 Lanczos (usando scipy.sparse.linalg.eigsh)

from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigsh
# Matriz esparcida grande (ejemplo 10000x10000)
N = 10000
A_sparse = csr_matrix(np.diag(np.arange(1, N+1)))
# Autovalores extremos (k=5)
vals, vecs = eigsh(A_sparse, k=5, which='LM')
print("5 mayores autovalores:", vals)

3.5 Arnoldi (usando scipy.sparse.linalg.eigs)

from scipy.sparse.linalg import eigs
# Matriz no simétrica aleatoria esparcida
A_non_sym = csr_matrix(np.random.rand(N, N))
vals, vecs = eigs(A_non_sym, k=6, which='LR')
print("6 autovalores con mayor parte real:", vals)

4. Buenas prácticas, rendimiento y escalabilidad

• Pre‑condicionamiento: Para Power Inverse y RQI, usar pre‑condicionadores (ILU, Jacobi) reduce drásticamente el número de iteraciones.
• Elección de la semilla: Un vector inicial aleatorio suele ser suficiente, pero si se conoce alguna dirección privilegiada, usarla acelera la convergencia.
• Uso de precisión doble vs. simple: En GPUs o aceleradores, el cálculo en precisión simple es más rápido; sin embargo, los métodos iterativos pueden ser sensibles a la pérdida de ortogonalidad. Aplicar re‑ortogonalización periódica (Gram‑Schmidt) cuando sea necesario.
• Paralelismo: Operaciones de producto matriz‑vector (A·v) son triviales de paralelizar con numpy.dot multihilo o con cupy en CUDA. Las factorizaciones LU para Inverse Power pueden usar scipy.sparse.linalg.splu con back‑ends multihilo.
• Escalabilidad: En matrices con más de 10⁶ filas, prefiera métodos Krylov (Lanczos/Arnoldi) que sólo almacenan una base de dimensión k (usualmente < 100) en lugar de mantener la matriz completa en memoria.

5. Troubleshooting frecuente

1. Convergencia lenta o estancada
   • Verifique que el autovalor objetivo sea aislado; si hay valores cercanos, el método puede oscilar.
   • Intente cambiar la semilla o aplicar un shift diferente en Inverse Power.
   • Para Lanczos, controle la pérdida de ortogonalidad usando reorthogonalize=True (en paquetes como scikit‑sparse).
2. Desbordamiento numérico
   • Normalice el vector en cada iteración (ya incluido en los ejemplos).
   • Use tipos de dato float64 o mayor precisión si el rango de valores es amplio.
3. Fallo al resolver (A‑μI)·y = x
   • El shift puede coincidir exactamente con un autovalor, haciendo singular la matriz. Ajuste ligeramente el valor de \(μ\).
   • Emplee un solver iterativo tolerante a indefinidos (GMRES) con pre‑condicionador.

6. Conclusiones

Los métodos iterativos ofrecen una vía eficaz para estimar autovalores en contextos de alta dimensionalidad y matrices esparcidas. La elección del algoritmo depende de la naturaleza de la matriz (simétrica vs. general), del autovalor de interés (extremo, cercano a un shift) y de los recursos computacionales disponibles. Con las implementaciones mostradas y las buenas prácticas descritas, puedes integrar rápidamente estimaciones de eigen‑valores en pipelines de análisis, modelado de sistemas dinámicos o algoritmos de aprendizaje automático.