import numpy as np
def one_term_matrix_equation_solver(A, B, F):
    ## ADD YOUR CODE HERE
    return np.zeros_like(F)

## CODE BELOW IS PROVIDED FOR YOU. NO NEED TO CHANGE IT.

n = 10
A = np.random.randn(n,n)
X = np.random.randn(n,n)
B = np.random.randn(n,n)
F = B @ X @ A.T

Xsolve = one_term_matrix_equation_solver(A, B, F)
# Should be close to machine epsilon ~1e-13 is fine
print('error: ', np.linalg.norm(Xsolve - X )/np.linalg.norm(X))

error:  1.0

import numpy as np

def kronecker_solver(Bs, As, F):
    ## ADD YOUR CODE HERE
    return np.zeros_like(F)

## CODE BELOW IS PROVIDED FOR YOU. NO NEED TO CHANGE IT.

## the vec(X) function implemented in python.
def vec(X):
    # Numpy is row-ordered, so we need to transpose before we flatten.
    return X.T.reshape(-1)

## Inverse of vec function.
def unvec(x):
    n = np.sqrt(x.size).astype(int)
    return x.reshape(n,n).T


n = 10
k = 5
X = np.random.randn(n,n)
As = np.random.randn(k,n,n) # k (nxn) matrices
Bs = np.random.randn(k, n,n)
F = np.zeros((n,n))
for i in range(k):
    F += Bs[i] @ X @ As[i].T

Xsolve = kronecker_solver(Bs, As, F)

# Should be close to machine epsilon ~1e-13 is fine
print('error: ', np.linalg.norm(Xsolve - X )/np.linalg.norm(X))

error:  1.0

import numpy as np

def fast_kron_matvec(A,B, x):
    ## ADD YOUR CODE HERE
    return np.zeros_like(x)



## CODE BELOW IS PROVIDED FOR YOU. NO NEED TO CHANGE IT.

def naive_kron_matvec(A,B,x):
    # A O(n^2) algorithm
    return np.kron(A,B) @ x

n_root = 30
n = n_root**2
A = np.random.randn(n_root,n_root) 
B = np.random.randn(n_root,n_root)
x = np.random.randn(n)

b = naive_kron_matvec(A,B,x)
b_fast = fast_kron_matvec(A,B, x)

# Should be close to machine epsilon ~1e-13 is fine
print('error: ', np.linalg.norm(b - b_fast )/np.linalg.norm(b))

error:  1.0

import scipy.linalg
import numpy as np
import time
import matplotlib.pyplot as plt

def fast_toeplitz_cholesky(c):
    n = c.size
    R = np.zeros((n,n))
    ## ADD YOUR CODE HERE

    return R

## CODE BELOW IS PROVIDED FOR YOU. NO NEED TO CHANGE IT.


# A helper function you may use to compute the Q matrix from part3.
def get_Q(gamma):
    if np.abs(gamma) >= 1:
        print('Something went wrong! |gamma| > 1')
    Q  = np.array( [[1, -gamma], [ -gamma, 1]] ) / ( np.sqrt( 1- gamma**2))
    return Q


# Form a random Toeplitz SPD matrix using its circulant embedding.
n = 10
fx = np.random.randn(10)**2
x = np.fft.ifft(fx)
C = scipy.linalg.circulant(x)
T = C[:n//2, :n//2].real
s,_ = np.linalg.eig(T)
print('T symmetric?: ', np.allclose(T, T.T) )
print('T positive definite?: ', np.all(s > 0) )

c = T[:,0]

R = fast_toeplitz_cholesky(c)
print('Backward error: ', np.linalg.norm(R.T @ R - T) / np.linalg.norm(T) )


# Running with smallish ns first to make sure the complexity is right. Running this whole cell takes ~1s on my machine. If it takes a lot longer, your solution might be O(n^3)
ns = 2**(np.arange(4,10))

# Once you have debugged, uncomment this line. Careful: if your algorithm is O(n^3), it will take very long.
# It takes ~ 4s on my machine to run with the O(n^2) alg.
# ns = 2**(np.arange(4,13))
times = []

for n in ns:
    fx = np.random.randn(2*n)**2
    x = np.fft.ifft(fx)
    C = scipy.linalg.circulant(x)
    T = C[:n, :n].real
    c = T[:,0]

    start_time = time.time()
    for _ in range(10):
        R = fast_toeplitz_cholesky(c)
    end_time = time.time()
    times.append(end_time - start_time)

    print(f"Time taken for n={n}: {end_time - start_time:.4f} seconds")

# On my machine, the plot looks even faster than O(n^2), but blows up after ~ 8000.
# The scaling is not always observed for smallish problem, where the real bottleneck is moving data to and from the CPU. 
# Plotting the results
plt.figure(figsize=(10, 6))
plt.loglog(ns, times, marker='o', label='Measured time')

plt.loglog(ns, [times[2] * (n / ns[2])**2 for n in ns], linestyle='--', label='O(n^2) reference')
plt.loglog(ns, [times[2] * (n / ns[2])**3 for n in ns], linestyle='--', label='O(n^3) reference')

plt.xlabel('Matrix size (n)')
plt.ylabel('Time (seconds)')
plt.title('Scaling of Toeplitz Cholesky Decomposition')
plt.legend()
plt.grid(True)
plt.show()

T symmetric?:  True
T positive definite?:  True
Backward error:  1.0
Time taken for n=16: 0.0000 seconds
Time taken for n=32: 0.0000 seconds
Time taken for n=64: 0.0000 seconds
Time taken for n=128: 0.0000 seconds
Time taken for n=256: 0.0000 seconds
Time taken for n=512: 0.0009 seconds

#If you don't have skimage, uncomment the line below and run it
# !pip install scikit-image

from skimage import data
import numpy as np
import matplotlib.pyplot as plt
## THIS WHOLE CELL IS PROVIDED FOR YOU. NO NEED TO CHANGE ANY CODE HERE.

def get_PSF_diagonal(defocus,n):
    X = np.mgrid[-n//2:n//2, -n//2:n//2]
    R = np.sqrt(X[0]**2 + X[1]**2)
    return np.fft.ifftshift(np.sin(- np.pi*defocus * R + 0.1))

def blur_image(image, kernel):
    return np.fft.ifft2(np.fft.fft2(image, norm = 'ortho') * kernel, norm = 'ortho')

def display_subimage(img, ax, title):
    im = ax.imshow(img.reshape(n,n))
    ax.set_title(title)
    ax.set_axis_off()
    plt.colorbar(im, ax=ax)

sigma = 1
tau = 1
clean_image = data.shepp_logan_phantom()[::4,::4]

n = 100
convolved_image = blur_image(clean_image, get_PSF_diagonal(0.1,100))
convolved_noisy_image = convolved_image + sigma * np.random.randn(n,n)

fig, axs = plt.subplots(1, 3,figsize = (20,5))
display_subimage(clean_image, axs[0], 'clean image x')
display_subimage(convolved_image.real, axs[1], 'convolved image C_i x')
display_subimage(convolved_noisy_image.real, axs[2], 'convolved image + noise y_i')

# ( You need to have run the previous cell before you run this one)

def average_and_deblur(images, ks, lambda_val = 1):
    ## ADD YOUR CODE HERE
    return np.zeros_like(images[0])


## CODE BELOW IS PROVIDED FOR YOU. NO NEED TO CHANGE IT.

n = 100 # Images of size n x n
m = 1000 # Number of images
noisy_images = np.zeros((m,n,n))
ks = np.zeros((m,n,n)) 
# diagonal of matrix D_i, C_i = F_2^{-1} D_i F_2
# The diagonal vectors of size (n^2) are stored as matrices of size (n x n) 

# Simulate some defocus values, which defines the point spread function.
defocus = np.random.rand(m)*0.1 + 0.05

tau = 1
lambda_val = sigma / tau
for i in range(m):
    # Get kernel
    ks[i] = get_PSF_diagonal(defocus[i],n)
    # Blur image
    noisy_images[i] = blur_image(clean_image, ks[i]).real
    # Add noise
    noisy_images[i] += np.random.randn(n,n) * sigma


averaged_deblurred_image = average_and_deblur(noisy_images, ks, lambda_val).real

fig, axs = plt.subplots(1, 3,figsize = (20,5))
display_subimage(clean_image, axs[0], 'clean image x')
display_subimage(noisy_images[0], axs[1], 'one  noisy image')
display_subimage(averaged_deblurred_image, axs[2], 'deblurred, denoised image')

# relative error should be around ~0.2 (there is noise, not just floating point error so we don't expect to get to O(machine epsilon)), and the recovered image should look pretty good! (still a bit noisy)
print('relative error:', np.linalg.norm(clean_image - averaged_deblurred_image) / np.linalg.norm(clean_image))

relative error: 1.0

import numpy as np
from scipy.linalg import cholesky

def block_cholesky_extend(R11, A12, A22):
    """Compute the Cholesky factorization of a block matrix A given R11, A12, and A22.
    
    Given a symmetric positive definite matrix A partitioned as:
        A = [[A11, A12],
             [A12.T, A22]]
    where A11 = R11.T @ R11 is the Cholesky factorization of A11,
    computes the Cholesky factorization R of the full matrix A.

    Args:
        R11 (ndarray): Upper triangular Cholesky factor of A11
        A12 (ndarray): Off-diagonal block of A
        A22 (ndarray): Lower right diagonal block of A
        
    Returns:
        ndarray: Upper triangular Cholesky factor R of the full matrix A
    """
    n, k = R11.shape[0], A22.shape[0]
    R = np.zeros((n + k, n + k), dtype=R11.dtype)

    ## ADD YOUR CODE HERE
    return R

rng = np.random.default_rng(0)
n1, n2 = 80, 5
M = rng.standard_normal((n1 + n2, n1 + n2))
A = M @ M.T + np.eye(n1 + n2)
A11 = A[:n1, :n1]
A12 = A[:n1, n1:]
A22 = A[n1:, n1:]
R11 = cholesky(A11, lower=False)
R_block = block_cholesky_extend(R11, A12, A22)
R_full = cholesky(A, lower=False)
print(f"relative factor error (Q4.2): {np.linalg.norm(R_block - R_full) / np.linalg.norm(R_full):.2e}")

relative factor error (Q4.2): 1.00e+00

## Note: This code is provided for you. No need to change it.

import numpy as np

def phi(r, ell=0.2):
    # Single fixed-bandwidth Gaussian kernel
    return np.exp(-(r*r)/(ell*ell))

def pairwise_dists(X, Y):
    """
    X: (d,n), Y: (d,m)
    returns: (n,m) with ||x_i - y_j||_2
    """
    X2 = np.sum(X*X, axis=0)[:,None]  # (n,1)
    Y2 = np.sum(Y*Y, axis=0)[None,:]  # (1,m)
    XY = X.T @ Y                      # (n,m)
    sq = X2 + Y2 - 2*XY
    np.maximum(sq, 0, out=sq)         # numerical safety
    return np.sqrt(sq)

def Phi_block(X, Y, ell=0.2):
    D = pairwise_dists(X, Y)
    return phi(D, ell=ell)

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import cholesky

def fit_rbf(X, f, ell=0.2):
    """
    Fits a Radial Basis Function interpolant by solving the linear system Phi_XX * c = f
    
    Args:
        X: (d,n) array of training points
        f: (n,) array of function values at training points
        ell: bandwidth parameter for the Gaussian kernel
        
    Returns:
        coeffs: (n,) array of RBF coefficients
        R: Upper triangular Cholesky factor of Phi_XX
    """
    n = X.shape[1]
    coeffs = np.zeros(n)
    R =  np.zeros((n,n))
    # ADD YOUR CODE HERE
    return coeffs, R

def evaluate_rbf(X, coeffs, X_query, ell=0.2):
    """
    Evaluates the RBF interpolant at query points
    
    Args:
        X: (d,n) array of training points
        coeffs: (n,) array of RBF coefficients
        X_query: (d,m) array of query points
        ell: bandwidth parameter for the Gaussian kernel
        
    Returns:
        (m,) array of interpolated values at query points
    """
    n = X.shape[1]
    f_at_query = np.zeros(X_query.shape[1])
    return f_at_query

rng = np.random.default_rng(1)
X_train = np.linspace(0.0, 1.0, 6)[None, :]
f_train = np.sin(2 * np.pi * X_train)[0]
coeffs, R = fit_rbf(X_train, f_train, ell=0.18)
Phi_train = Phi_block(X_train, X_train, ell=0.18)
residual = np.linalg.norm(evaluate_rbf(X_train, coeffs, X_train, ell=0.18) - f_train) / np.linalg.norm(f_train)
print(f"relative training residual (Q4.3): {residual:.2e}")

x_dense = np.linspace(0.0, 1.0, 400)[None, :]
s_dense = evaluate_rbf(X_train, coeffs, x_dense, ell=0.18)
fig, ax = plt.subplots(figsize=(6.0, 3.0))
ax.plot(x_dense.ravel(), np.sin(2 * np.pi * x_dense).ravel(), label='ground truth', lw=2)
ax.plot(x_dense.ravel(), s_dense, label='rbf interpolant', lw=2, linestyle='--')
ax.scatter(X_train.ravel(), f_train, color='black', s=20, label='data points')
ax.set_title('Gaussian RBF interpolation (Question 4.3)')
ax.set_xlabel('x')
ax.set_ylabel('f(x)')
ax.legend(loc='upper right')
ax.grid(alpha=0.25)
fig.tight_layout()
plt.show()

relative training residual (Q4.3): 1.00e+00

import numpy as np
import matplotlib.pyplot as plt


def update_rbf_interpolant(X_old, f_old, R_old, X_new, f_new, ell=0.2):
    """
    Updates an existing RBF interpolant with new data points efficiently.

    Args:
        X_old (ndarray): Matrix of existing points, shape (d, n)
        f_old (ndarray): Vector of existing function values, shape (n,)
        R_old (ndarray): Upper triangular Cholesky factor of existing RBF matrix
        X_new (ndarray): Matrix of new points to add, shape (d, p) 
        f_new (ndarray): Vector of new function values, shape (p,)
        ell (float): Length scale parameter for RBF kernel

    Returns:
        tuple:
            - coeffs_ext (ndarray): Updated RBF coefficients for all points
            - R_ext (ndarray): Updated Cholesky factor
            - X_ext (ndarray): Combined matrix of all points
    """
    # ADD YOUR CODE HERE
    n_old = X_old.shape[1]
    k_new = X_new.shape[1]
    coeffs_ext, R_ext, X_ext =  np.zeros(n_old + k_new), np.zeros((n_old + k_new, n_old + k_new)), np.zeros((d, n_old + k_new))
    return coeffs_ext, R_ext, X_ext



## The code below is provided for you. No need to change it.
ell = 0.2
rng = np.random.default_rng(4)
d, n_old, k_new = 2, 6, 3
X_old = rng.random((d, n_old))
X_new = rng.random((d, k_new))
ffun = lambda x: np.cos(2 * np.pi * x[0]) * np.sin(2 * np.pi * x[1])
f_old = ffun(X_old)
f_new = ffun(X_new)
R_old = np.linalg.cholesky(Phi_block(X_old, X_old, ell=ell)).T
coeffs_ext, R_ext, X_ext = update_rbf_interpolant(X_old, f_old, R_old, X_new, f_new, ell=ell)

# Estimate the naive way to do it
X_ext_2 = np.concatenate([X_old, X_new], axis=1)
Phi_full = Phi_block(X_ext_2, X_ext_2, ell=ell)
coeffs_ref = np.linalg.solve(Phi_full, np.concatenate([f_old, f_new]))
print(f"relative coefficient error (Q4.4): {np.linalg.norm(coeffs_ext - coeffs_ref) / np.linalg.norm(coeffs_ref):.2e}")

# A visualization of the update
x_plot = np.linspace(0.0, 1.0, 1000)[None, :]
fun1d = lambda x: np.sin(6 * np.pi * x[0]) + 0.3 * np.cos(2 * np.pi * x[0])
X_old_1d = np.linspace(0.0, 1.0, 7)[None, :]
X_new_1d = np.array([[0.18, 0.74]])
f_old_1d = fun1d(X_old_1d)
f_new_1d = fun1d(X_new_1d)
coeffs_old_1d, R_old_1d = fit_rbf(X_old_1d, f_old_1d, ell=ell)
coeffs_new_1d, R_new_1d, X_ext_1d = update_rbf_interpolant(X_old_1d, f_old_1d, R_old_1d, X_new_1d, f_new_1d, ell=ell)
interp_old = evaluate_rbf(X_old_1d, coeffs_old_1d, x_plot, ell=ell)
interp_new = evaluate_rbf(X_ext_1d, coeffs_new_1d, x_plot, ell=ell)

fig, ax = plt.subplots(figsize=(6.2, 3.0))
ax.plot(x_plot.ravel(), fun1d(x_plot), label='ground truth', lw=2)
ax.plot(x_plot.ravel(), interp_old, label='before update', linestyle='--', lw=2)
ax.plot(x_plot.ravel(), interp_new, label='after update', linestyle='-.', lw=2)
ax.scatter(X_old_1d.ravel(), f_old_1d, color='black', s=28, label='old data', zorder=3)
ax.scatter(X_new_1d.ravel(), f_new_1d, color='tab:orange', s=38, label='new data', zorder=3)
ax.set_title('Effect of batch update (Question 4.4)')
ax.set_xlabel('x')
ax.set_ylabel('f(x)')
ax.legend(loc='best')
ax.grid(alpha=0.25)
fig.tight_layout()
plt.show()

relative coefficient error (Q4.4): 1.00e+00

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import cholesky
plt.style.use('seaborn-v0_8-colorblind')

def fill_single_missing_value(R, f_values, missing_index):
    """
    Assuming that entry missing_index of the measurement vector fX is missing,
    fill it with an approximate value from interpolating the remaining points
    using block elimination.
    
    Args:
        R: Upper triangular Cholesky factor of the RBF kernel matrix
        f_values: Array of function values, including corrupted value
        missing_index: Index of the corrupted/missing value to repair
        
    Returns:
        repaired_value: The repaired/interpolated value at the missing index
        coeffs: RBF coefficients with the missing point zeroed out
    """
    n = R.shape[0]
    repaired_value, coeffs = 1, np.zeros(n)
    return repaired_value, coeffs

## The code below is provided for you. No need to change it.

# Set random seed and generate test data
rng = np.random.default_rng(5)
d = 2
n = 45
X = rng.random((d, n))
true_fun = lambda x: np.exp(x[0]) * np.cos(3 * x[1])
f_full = true_fun(X)

# Test 2D case
Phi_full = Phi_block(X, X, ell=0.27)
R_full = cholesky(Phi_full, lower=False)
missing_index = 6

# Add corruption and repair
f_corrupted = f_full.copy()
f_corrupted[missing_index] += 0.25
repaired_value, coeffs_hat = fill_single_missing_value(R_full, f_corrupted, missing_index)

# Compare with reference solution
mask = np.ones(n, dtype=bool)
mask[missing_index] = False
Phi_reduced = Phi_full[np.ix_(mask, mask)]
f_reduced = f_full[mask]
coeffs_reduced = fit_rbf(X[:,mask], f_reduced, ell=0.27)[0]
predicted_ref = evaluate_rbf(X[:,mask], coeffs_reduced, X[:,[missing_index]], ell=0.27)
print(f"Relative error in repaired value: {np.linalg.norm(repaired_value - predicted_ref.item()) / np.linalg.norm(predicted_ref.item()):.6e}")

# Test 1D case for visualization
x_plot = np.linspace(0.0, 1.0, 400)[None, :]
true_fun_1d = lambda x: np.sin(4 * np.pi * x)
X_1d = np.linspace(0, 1, 20)[None, :]
f_true_1d = true_fun_1d(X_1d).ravel()

# Set up 1D problem
Phi_1d = Phi_block(X_1d, X_1d, ell=0.18)
R_1d = cholesky(Phi_1d, lower=False)
corrupt_idx = 6

# Add corruption and repair for 1D case
f_corrupted_1d = f_true_1d.copy()
f_corrupted_1d[corrupt_idx] += 0.3
coeffs_bad = fit_rbf(X_1d, f_corrupted_1d, ell=0.18)[0]
repaired_value_1d, coeffs_fixed = fill_single_missing_value(R_1d, f_corrupted_1d, corrupt_idx)

# Compute interpolations
interp_bad = evaluate_rbf(X_1d, coeffs_bad, x_plot, ell=0.18)
interp_fixed = evaluate_rbf(X_1d, coeffs_fixed, x_plot, ell=0.18)

# Create visualization
fig, axes = plt.subplots(1, 2, figsize=(10.0, 3.2), sharey=True)

# Left plot - Before repair
axes[0].plot(x_plot.ravel(), true_fun_1d(x_plot).ravel(), label='ground truth', lw=2)
axes[0].plot(x_plot.ravel(), interp_bad, label='using corrupted data', linestyle='--', lw=2)
axes[0].scatter(X_1d.ravel(), f_corrupted_1d, color='tab:red', s=32, alpha=0.6, zorder=3)
axes[0].scatter(X_1d[0, corrupt_idx], f_corrupted_1d[corrupt_idx], 
                s=140, facecolors='none', edgecolors='red', linewidths=2.2, 
                label='corrupted point', zorder=4)
axes[0].annotate('corrupted', 
                xy=(X_1d[0, corrupt_idx], f_corrupted_1d[corrupt_idx]),
                xytext=(X_1d[0, corrupt_idx] + 0.07, f_corrupted_1d[corrupt_idx] + 0.6),
                arrowprops=dict(arrowstyle='->', color='red', linewidth=1.5),
                color='red', fontsize=9)
axes[0].set_title('Before repair')
axes[0].set_xlabel('x')
axes[0].set_ylabel('f(x)')
axes[0].legend(loc='upper right')
axes[0].grid(alpha=0.25)

# Right plot - After repair
f_repaired = f_corrupted_1d.copy()
f_repaired[corrupt_idx] = repaired_value_1d
axes[1].plot(x_plot.ravel(), true_fun_1d(x_plot).ravel(), label='ground truth', lw=2)
axes[1].plot(x_plot.ravel(), interp_fixed, label='after repair', linestyle='-.', lw=2)
axes[1].scatter(X_1d.ravel(), f_repaired, color='tab:green', s=32, alpha=0.6, zorder=3)
axes[1].scatter(X_1d[0, corrupt_idx], repaired_value_1d,
                s=140, facecolors='none', edgecolors='black', linewidths=2.2,
                label='repaired value', zorder=4)
axes[1].set_title('After repair')
axes[1].set_xlabel('x')
axes[1].legend(loc='upper right')
axes[1].grid(alpha=0.25)

# Finalize plot
fig.suptitle('Effect of single-point repair (Question 4.5)', fontsize=13)
fig.tight_layout()
plt.show()

Relative error in repaired value: inf

/var/folders/0b/bgs7xkxj0cj4jc50c5v9hvwc0000gq/T/ipykernel_11682/370252130.py:52: RuntimeWarning: divide by zero encountered in scalar divide
  print(f"Relative error in repaired value: {np.linalg.norm(repaired_value - predicted_ref.item()) / np.linalg.norm(predicted_ref.item()):.6e}")

MAT 321 - PS3¶

ADD YOUR CODE HERE¶

Q1 Matrix equations¶

Q1.1 Not-Kronecker solve¶

Q 1.2 Kronecker solve¶

Answer: 1.2 write the complexity here¶

Q1.3 Fast Kronecker matvec¶

Q2: Cholesky Decomposition of Toeplitz matrices¶

Q2.1-Q2.4 Add your answers to parts 1-4 here¶

Q3: 2D averaging and deblurring¶

Q3.1¶

3.1 Add your answer here¶

Q3.2¶

3.2 Add your answer here¶

Q3.3¶

Q4 Cholesky update and RBF¶

4.1 - Schur complement of symmetric positive definite matrices¶

Add your solution here¶

4.2 Update Cholesky factorization¶

4.3 RBF Interpolation¶

4.4¶

4.5: Corrupted data¶

¶

Extra credit¶