Davidson Method

import numpy as np

def gram_schmidt_orthogonalize(w, V):
    """
    Orthogonalize vector w against the columns of V using Gram-Schmidt.
    
    Parameters:
        w: Vector to orthogonalize
        V: Matrix whose columns form the current subspace
    
    Returns:
        Orthogonalized and normalized vector
    """
    w = w.copy().reshape(-1)
    
    # Subtract projections onto each vector in V
    for i in range(V.shape[1]):
        v = V[:, i].reshape(-1)
        projection = np.dot(w, v) / np.dot(v, v) * v
        w -= projection
    
    # Normalize
    norm = np.linalg.norm(w)
    if norm > 1e-12:
        w /= norm
    else:
        # If vector is nearly zero, return a random vector
        w = np.random.randn(len(w))
        w /= np.linalg.norm(w)
    
    return w.reshape(-1, 1)

def davidson_method(A, num_eigenvalues=1, max_subspace_dim=50, tol=1e-6, max_iter=100):
    """
    Davidson method for computing the smallest eigenvalues of a symmetric matrix.
    
    Parameters:
        A: Symmetric matrix (n x n)
        num_eigenvalues: Number of eigenvalues to compute
        max_subspace_dim: Maximum dimension of the search subspace
        tol: Convergence tolerance
        max_iter: Maximum number of iterations
    
    Returns:
        eigenvalues: Converged eigenvalues
        eigenvectors: Corresponding eigenvectors
        iterations: Number of iterations performed
    """
    n = A.shape[0]
    assert A.shape[0] == A.shape[1], "Matrix must be square"
    assert np.allclose(A, A.T), "Matrix must be symmetric"
    
    # Diagonal preconditioner
    D = np.diag(np.diag(A))
    
    # Initial subspace: random normalized vector
    V = np.random.randn(n, 1)
    V /= np.linalg.norm(V)
    
    eigenvalues = []
    eigenvectors = []
    
    for iteration in range(max_iter):
        # Project matrix onto subspace
        H = V.T @ A @ V
        
        # Solve eigenvalue problem in subspace
        eigvals, eigvecs = np.linalg.eigh(H)
        
        # Sort eigenvalues (ascending)
        idx = np.argsort(eigvals)
        eigvals = eigvals[idx]
        eigvecs = eigvecs[:, idx]
        
        # Check convergence for each desired eigenvalue
        converged_count = 0
        for i in range(min(num_eigenvalues, len(eigvals))):
            # Ritz value and vector
            mu = eigvals[i]
            y = eigvecs[:, i]
            x = V @ y  # Expand to full space
            
            # Compute residual
            residual = (mu * np.eye(n) - A) @ x
            residual_norm = np.linalg.norm(residual)
            
            if residual_norm < tol:
                # Converged eigenvalue
                if i >= len(eigenvalues):
                    eigenvalues.append(mu)
                    eigenvectors.append(x)
                converged_count += 1
            else:
                # Add correction vector to subspace
                if V.shape[1] < max_subspace_dim:
                    # Correction vector: (mu*I - D)^(-1) * r
                    correction = residual / (mu - np.diag(D) + 1e-10)
                    
                    # Orthogonalize against existing subspace
                    correction_ortho = gram_schmidt_orthogonalize(correction, V)
                    
                    # Add to subspace
                    V = np.hstack([V, correction_ortho])
        
        # If all eigenvalues converged, stop
        if converged_count >= num_eigenvalues:
            break
        
        # Limit subspace size (restart if needed)
        if V.shape[1] >= max_subspace_dim:
            # Keep only the best vectors
            V = V[:, :num_eigenvalues]
    
    return np.array(eigenvalues), np.array(eigenvectors).T, iteration + 1

# Example: Diagonally dominant matrix
def generate_diagonally_dominant_matrix(n, sparsity=1e-4, seed=None):
    """Generate a diagonally dominant symmetric matrix."""
    if seed is not None:
        np.random.seed(seed)
    
    # Create matrix with diagonal elements 1, 2, ..., n
    A = np.zeros((n, n))
    for i in range(n):
        A[i, i] = i + 1
    
    # Add small random off-diagonal elements
    A = A + sparsity * np.random.randn(n, n)
    A = (A + A.T) / 2  # Symmetrize
    
    return A

# Test the method
np.random.seed(42)
n = 200
A = generate_diagonally_dominant_matrix(n, sparsity=1e-3)

print("Computing smallest eigenvalues using Davidson method...")
eigenvalues, eigenvectors, iterations = davidson_method(A, num_eigenvalues=4, tol=1e-8)

print(f"\nConverged in {iterations} iterations")
print(f"\nSmallest eigenvalues:")
for i, ev in enumerate(eigenvalues):
    print(f"  λ_{i+1} = {ev:.8f}")

# Compare with exact eigenvalues
exact_eigenvalues = np.sort(np.linalg.eigvalsh(A))[:4]
print(f"\nExact eigenvalues (from full diagonalization):")
for i, ev in enumerate(exact_eigenvalues):
    print(f"  λ_{i+1} = {ev:.8f}")

print(f"\nErrors:")
for i in range(len(eigenvalues)):
    error = abs(eigenvalues[i] - exact_eigenvalues[i])
    print(f"  |λ_{i+1} - λ_{i+1}^exact| = {error:.2e}")