popipopi

import numpy as np

def outerProductCholesky(A, tol=1e-10):
    """
    Computes the outer product Cholesky decomposition of a symmetric and 
    positive definite matrix A.
    
    Parameters:
        A : numpy.ndarray
            A symmetric matrix.
        tol : float
            Tolerance for symmetry check and for deciding positivity.
            
    Returns:
        L : numpy.ndarray
            Lower-triangular matrix such that A = L @ L.T
            
    Raises:
        ValueError: if A is not symmetric positive definite.
    """
    A = np.array(A, dtype=float)
    n, m = A.shape
    if n != m:
        raise ValueError("Matrix A must be square.")
    
    # Check for symmetry
    if not np.allclose(A, A.T, atol=tol):
        raise ValueError("Matrix A is not symmetric.")
    
    L = np.zeros_like(A)
    
    # Outer product algorithm:
    for k in range(n):
        # Compute the pivot value
        temp = A[k, k] - np.sum(L[k, :k]**2)
        if temp <= tol:
            raise ValueError("Matrix is not positive definite. "
                             "Encountered nonpositive pivot at index {}.".format(k))
        L[k, k] = np.sqrt(temp)
        
        # Compute the k-th column below the diagonal.
        for i in range(k+1, n):
            # Using the formula: L[i,k] = (A[i,k] - sum_{j=0}^{k-1} L[i,j] * L[k,j]) / L[k,k]
            L[i, k] = (A[i, k] - np.sum(L[i, :k] * L[k, :k])) / L[k, k]
            
    return L

# Define the matrices from the assignment:

# Matrix A from Problem 1
A = np.array([[4, 0, -2],
              [0, 1,  4],
              [-2, 4, 21]], dtype=float)

# Matrix B from point 5
B = np.array([[-1, 2, 3],
              [ 2,-2, 5],
              [ 3, 5,-1]], dtype=float)

# Matrix C from point 5
C = np.array([[ 4, -2, 6],
              [-2, 17, 1],
              [ 6,  1,11]], dtype=float)

# --- Testing the function ---

print("Cholesky decomposition for matrix A:")
try:
    L_A = outerProductCholesky(A)
    print("L_A =\n", L_A)
    # To verify, we check that L_A @ L_A.T equals A
    print("Reconstruction A from L_A @ L_A.T:\n", np.dot(L_A, L_A.T))
except ValueError as e:
    print("Error:", e)

print("\nCholesky decomposition for matrix B:")
try:
    L_B = outerProductCholesky(B)
    print("L_B =\n", L_B)
except ValueError as e:
    print("Error:", e)

print("\nCholesky decomposition for matrix C:")
try:
    L_C = outerProductCholesky(C)
    print("L_C =\n", L_C)
    print("Reconstruction C from L_C @ L_C.T:\n", np.dot(L_C, L_C.T))
except ValueError as e:
    print("Error:", e)