Eritque Arcus Math Notes

Krylov subspace

2 min read

shifting does not improve Lanczos / Arnoldi

Concept

Lanczos iteration for symmetric

import numpy as np
 
def lanczos_iter(A, tol=1e-8):
    n = A.shape[0]
    
    # q0 = np.zeros(n)  # q0 = 0
    # b0 = 0            # b0 = 0
    x0 = np.random.rand(n)
    q1 = x0 / np.linalg.norm(x0)
	
    Q = [q1]  # start with k = 1
    alphas = []
    betas = [b0]
	
    qk = q1
    qk_prev = q0
    bk_prev = b0
    
    for k in range(1, n+1):
        # u_k = A*q_k
        uk = A @ qk
		
		# .T same as .H cause it's real
        # a_k = q_k^T * u_k
        ak = qk.T @ uk
        alphas.append(ak)
        
        # u_k = u_k - b_{k-1}*q_{k-1} - a_k*q_k
        uk = uk - bk_prev * qk_prev - ak * qk
        
        # b_k = ||u_k||
        bk = np.linalg.norm(uk)
	    
        if abs(bk) < tol:
            break
            
        betas.append(bk)
        
        # q_{k+1} = u_k/b_k
        qk_next = uk / bk
        Q.append(qk_next)
        
	    
        qk_prev = qk
        qk = qk_next
        bk_prev = bk
    
    
    m = len(alphas)
    T = np.zeros((m, m))
    np.fill_diagonal(T, alphas)
    for i in range(m-1):
        T[i, i+1] = betas[i+1]
        T[i+1, i] = betas[i+1]
	    
    return np.column_stack(Q), T

Optimization

turn to Symmetric positive definite matrix.

Arnoldi iteration for nonsymmetric

import numpy as np
 
def arnoldi_iter(A, tol=1e-8):
    n = A.shape[0]
    x0 = np.random.rand(n)
    q1 = x0 / np.linalg.norm(x0)
    
    # Initialize properly
    Q = np.zeros((n, n+1))
    Q[:, 0] = q1
    H = np.zeros((n+1, n))
    
    for k in range(n):
        v = A @ Q[:, k]
        
        # Orthogonalization
        for j in range(k+1):
            H[j, k] = Q[:, j].T @ v
            v = v - H[j, k] * Q[:, j]
        
        H[k+1, k] = np.linalg.norm(v)
        
        if H[k+1, k] < tol:
            break
            
        Q[:, k+1] = v / H[k+1, k]
    
    return Q[:, :k+1], H[:k+2, :k+1]

Graph View