solve_lyapunov.py

def lyapunov_spectral(A, B, cond=None):
    u.check_symmetric(A)
    u.check_symmetric(B)

    s, U = torch.symeig(A, eigenvectors=True)
    if cond is None:
        cond = get_condition(s.dtype)
    cutoff = cond * max(s)
    s = torch.where(s > cutoff, s, torch.tensor(0.).to(s.device))

    C = U.t() @ B @ U    # TODO(y): throw away eigenvectors corresponding to discarded evals. U=U[:num_eigs]
    s = s.unsqueeze(1) + s.unsqueeze(0)
    si = torch.where(s > 0, 1 / s, s)
    Y = C * si
    X = U @ Y @ U.t()

    # cancel small asymetries introduces by multiplication by small numbers
    X = (X + X.t()) / 2
    return X