Ogita Aishima

ogita_aishima.py

import numpy as np


def estimate_eigeinvalues(A, X, extra=False):
    """From estimated eigenvectors X of A, compute estimated eigenvalues of A."""
    n, _ = A.shape
    R = np.eye(n) - X.T @ X
    S = X.T @ A @ X
    w = np.diag(S) / (1 - np.diag(R))
    if extra:
        return w, R, S
    return w


def ogita_aishima_step(A, X):
    """Perform one step of the Ogita-Aishima method."""
    w, R, S = estimate_eigeinvalues(A, X, True)
    D = np.diag(w)
    delta = 2 * (
        np.linalg.norm(S - D, ord=2)
        + np.linalg.norm(A, ord=2) * np.linalg.norm(R, ord=2)
    )
    ll = w - w[:, None]
    lr = R * w
    e = np.where(np.abs(ll) > delta, (S + lr) / ll, R / 2)
    return X + X @ e