An eigendecomposition that ensures eigenvectors are perpendicular.

    numpy.linalg.eig doesn't guarantee that eigenvectors from the same
    eigenspace will be perpendicular. This method uses Gram-Schmidt to recover
    a perpendicular set. It further checks that all eigenvectors are
    perpendicular and raises an ArithmeticError otherwise.

    Args:
        matrix: The matrix to decompose.
        rtol: Relative threshold for determining whether eigenvalues are from
              the same eigenspace and whether eigenvectors are perpendicular.
        atol: Absolute threshold for determining whether eigenvalues are from
              the same eigenspace and whether eigenvectors are perpendicular.

    Returns:
        The eigenvalues and column eigenvectors. The i'th eigenvalue is
        associated with the i'th column eigenvector.

    Raises:
        ArithmeticError: Failed to find perpendicular eigenvectors.

def _perp_eigendecompose(matrix: np.ndarray,
                         rtol: float = 1e-5,
                         atol: float = 1e-8,
                         ) -> Tuple[np.array, List[np.ndarray]]:
    """An eigendecomposition that ensures eigenvectors are perpendicular.

    numpy.linalg.eig doesn't guarantee that eigenvectors from the same
    eigenspace will be perpendicular. This method uses Gram-Schmidt to recover
    a perpendicular set. It further checks that all eigenvectors are
    perpendicular and raises an ArithmeticError otherwise.

    Args:
        matrix: The matrix to decompose.
        rtol: Relative threshold for determining whether eigenvalues are from
              the same eigenspace and whether eigenvectors are perpendicular.
        atol: Absolute threshold for determining whether eigenvalues are from
              the same eigenspace and whether eigenvectors are perpendicular.

    Returns:
        The eigenvalues and column eigenvectors. The i'th eigenvalue is
        associated with the i'th column eigenvector.

    Raises:
        ArithmeticError: Failed to find perpendicular eigenvectors.
    """
    vals, cols = np.linalg.eig(matrix)
    vecs = [cols[:, i] for i in range(len(cols))]

    # Convert list of row arrays to list of column arrays.
    for i in range(len(vecs)):
        vecs[i] = np.reshape(vecs[i], (len(vecs[i]), vecs[i].ndim))

    # Group by similar eigenvalue.
    n = len(vecs)
    groups = _group_similar(
        list(range(n)),
        lambda k1, k2: np.allclose(vals[k1], vals[k2], rtol=rtol))

    # Remove overlap between eigenvectors with the same eigenvalue.
    for g in groups:
        q, _ = np.linalg.qr(np.hstack([vecs[i] for i in g]))
        for i in range(len(g)):
            vecs[g[i]] = q[:, i]

    return vals, vecs