1. Fundamentos teóricos

El algoritmo de potencias (Power Iteration) es uno de los métodos iterativos más simples para obtener el eigenvalor dominante (de mayor magnitud) y su eigenvector asociado de una matriz cuadrada A\in\mathbb{R}^{n\times n}. Partiendo de un vector inicial v₀ no ortogonal al subespacio propio del eigenvalor dominante, la sucesión

v_{k+1}=\frac{A v_{k}}{\|A v_{k}\|}

converge a la dirección del eigenvector dominante u₁. El eigenvalor se puede estimar mediante el cociente de Rayleigh:

λ_{k}=\frac{v_{k}^{T} A v_{k}}{v_{k}^{T} v_{k}}

El algoritmo presenta convergencia lineal con razón |λ₂/λ₁|, donde λ₁ y λ₂ son el mayor y segundo mayor eigenvalor en módulo, respectivamente.

2. Implementación paso a paso en Python

Utilizaremos numpy para operaciones vector‑matriz y timeit para medir rendimiento. El código está estructurado para ser reutilizable y fácil de depurar.

import numpy as np
import time
def power_method(A, num_iters=1000, tol=1e-10, verbose=False):
    """Estimación del eigenvalor dominante y su eigenvector.
    Parámetros
    ----------
    A : np.ndarray
        Matriz cuadrada (n x n).
    num_iters : int
        Número máximo de iteraciones.
    tol : float
        Tolerancia para la convergencia del eigenvalor.
    verbose : bool
        Si True muestra el progreso en cada iteración.
    Returns
    -------
    lambda_est : float
        Aproximación del eigenvalor dominante.
    v : np.ndarray
        Eigenvector unitario correspondiente.
    """
    n = A.shape[0]
    # Vector inicial aleatorio (no ortogonal al eigenvector dominante)
    v = np.random.rand(n)
    v = v / np.linalg.norm(v)
    lambda_old = 0.0
    for k in range(num_iters):
        w = A @ v               # Producto matriz‑vector
        v = w / np.linalg.norm(w)
        lambda_new = v.T @ A @ v   # Cociente de Rayleigh
        if verbose:
            print(f"Iter {k+1:4d}: λ = {lambda_new:.12f}")
        if np.abs(lambda_new - lambda_old) < tol:
            break
        lambda_old = lambda_new
    return lambda_new, v
# ---- Ejemplo de uso ----------------------------------------------------
if __name__ == "__main__":
    # Matriz de prueba (simétrica positiva definida)
    A = np.array([[4, 1, 0],
                  [1, 3, 1],
                  [0, 1, 2]], dtype=float)
    start = time.time()
    lam, vec = power_method(A, num_iters=500, tol=1e-12, verbose=True)
    elapsed = time.time() - start
    print(f"\nEigenvalor dominante: {lam:.12f}")
    print(f"Eigenvector (unitario): {vec}")
    print(f"Tiempo de ejecución: {elapsed:.4f} s")

El bloque if __name__ == "__main__" permite ejecutar el script directamente o importarlo como módulo en notebooks.

3. Power Method con Shift (inverso) para el eigenvalor más pequeño

Para obtener el eigenvalor de menor magnitud, se aplica el algoritmo al inverso de A - σI. En la práctica, se resuelve el sistema lineal (A - σI) y = v en cada iteración.

from scipy.sparse.linalg import spsolve
def inverse_power_method(A, sigma=0.0, num_iters=1000, tol=1e-10, verbose=False):
    n = A.shape[0]
    I = np.eye(n)
    M = A - sigma * I
    v = np.random.rand(n)
    v = v / np.linalg.norm(v)
    lambda_old = 0.0
    for k in range(num_iters):
        # Resolver M y = v en lugar de multiplicar
        y = spsolve(M, v)
        v = y / np.linalg.norm(y)
        lambda_new = v.T @ A @ v
        if verbose:
            print(f"Iter {k+1:4d}: λ = {lambda_new:.12f}")
        if np.abs(lambda_new - lambda_old) < tol:
            break
        lambda_old = lambda_new
    return lambda_new, v

Esta variante es útil en problemas donde el eigenvalor más pequeño determina la estabilidad del sistema (por ejemplo, análisis de condición).

4. Buenas prácticas, troubleshooting y optimización

• Elección del vector inicial: Evite vectores colineales con eigenvectores asociados a eigenvalores de magnitud similar; una pequeña perturbación aleatoria suele ser suficiente.
• Escalado: Normalice el vector en cada iteración para prevenir desbordamiento/underflow numérico.
• Detección de convergencia: Use la diferencia absoluta del eigenvalor (abs(λ_k-λ_{k-1})) o la norma del residuo ||Av-λv||.
• Precisión de punto flotante: En matrices muy mal condicionadas, considere numpy.float64 o mpmath de precisión arbitraria.
• Rendimiento: Para matrices dispersas, utilice scipy.sparse y la función dot optimizada.
• Paralelismo: El producto matriz‑vector es fácilmente paralelizables con numba o cupy (GPU).
• Fallos comunes:
  • Convergencia lenta cuando |λ₂/λ₁|≈1. Solución: aplicar un shift o usar métodos de subespacio (Lanczos).
  • División por cero al normalizar v si Av=0. Verifique que A no sea singular o elija otro σ.

5. Casos de uso del mundo real