import numpy as np
import matplotlib.pyplot as plt

def get_central_diff_matrix(n):
    '''
    n = number of interior points, h = 1/(n+1).
    '''
    ###ADD YOUR CODE HERE
    A = np.zeros((n,n))
    return A 


## ALL THE CODE BELOW IS PROVIDED FOR YOU.

# Computes the grid [ 1/(n+1), 2/(n+1), ..., n/(n+1)]
def get_grid(n):
    return np.linspace(start = 0, stop = 1, num = n+2)[1:-1]


# Computes the grid [0, 1/(n+1), 2/(n+1), ..., n+1/(n+1)]
# A test case, u = sin(k pi x ), with known second derivative 
k = 3
u_fn = lambda x : np.sin(k * np.pi * x)
f_fn = lambda x : (k*np.pi)**2 * np.sin(k * np.pi*x)
fine_n = int(1e6)
fine_grid = get_grid(fine_n)


# Plots successive finite difference discretization for increasing number of points
plt.figure()
plt.plot(fine_grid, f_fn(fine_grid), label = 'true', color = 'black')
for n in [5, 10,  20]:
    A = get_central_diff_matrix(n)
    grid = get_grid(n)
    v = u_fn(grid)
    b = A @ v
    plt.plot(grid, b, '-o', alpha = 0.5, label = 'n ='+ str(n))
    
plt.legend()
plt.ylabel('$-d^2 f(x)$')
plt.xlabel('x')
plt.title('-$d^2f$ and its approximations')


# Plot error as a function of n
ns = np.logspace(2, 4, 10).astype(int)
errors = np.zeros(10)
for idx, n in enumerate(ns):
    h = 1 / (n + 1)
    A = get_central_diff_matrix(n) 


    grid = get_grid(n)
    u = u_fn(grid)
    b = A @ u 
    exact_df2 =  f_fn(grid)
    errors[idx] = np.max(np.abs(exact_df2 - b))
    
hs = 1 / ( ns + 1)
plt.figure()
plt.loglog( hs, errors, label = 'error')

# Make a line with slope O(h^2), that intersect the timing plot at last index
# The error should decay at the same rate at this line.
line = hs**2 
scale = errors[0] / line[0] 
plt.loglog( hs, scale * line, '--', label = 'O(h^2)')
plt.legend()
plt.ylabel('$e(h) = |u'' - D_h^2 u|_{inf}$')
plt.xlabel('h')
plt.gca().invert_xaxis()

n_tests = 10
ns = np.logspace(2, 4, n_tests).astype(int)
errors = np.zeros(n_tests)
### ADD YOUR CODE HERE
for idx, n in enumerate(ns):
    # Solve ODE of with grid size n. 
    

    # Compute relative error in the solution 
    errors[idx] = 0 

# produce a log-log scaling plot
# As above, plot a line O(h^2)

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

def fast_truncated_svd_matvec(Uk, Sk, Vk, x):
    # Takes in:
    # Uk: a nxk matrix 
    # Sk: a vector of length k
    # Vk: a nxk matrix
    # x: a vector of length n

    ###ADD YOUR CODE HERE

    return np.zeros_like(x)


## ALL THE CODE BELOW IS PROVIDED.
def slow_truncated_svd_matvec(Uk, Sk, Vk, x):
    # Takes in:
    # Uk: a nxk matrix 
    # Sk: a vector of length k
    # Vk: a nxk matrix
    # x: a vector of length n
    # Computes A_k x the slow way:
    A_k = (Uk @ np.diag(Sk)) @ Vk.T 
    return A_k @ x


# Testing for correctness
n = 100
A = np.random.randn(n, n)
[U,S,Vh] = np.linalg.svd(A)
k = 10 #rank 10 truncated SVD
Uk = U[:,:k]; Sk = S[:k]; Vk = Vh[:k,:].T

x = np.random.randn(n)
error = np.linalg.norm(slow_truncated_svd_matvec(Uk, Sk, Vk, x) - fast_truncated_svd_matvec(Uk, Sk, Vk, x) ) / np.linalg.norm(slow_truncated_svd_matvec(Uk, Sk, Vk, x) )
# Should be close to machine epsilon
print('error in truncated SVD mat-vec:', error)

# Testing complexity
ns = np.logspace(3,5.5, 30).astype(int)
timing = np.zeros(30)
for n_idx,n in enumerate(ns):
    # "Fake rank 10 SVD"

    # These matrices may not have orthonormal columns, but for checking scaling, it doesn't matter
    Uk = np.random.randn(n, 10) 
    Vk = np.random.randn(n, 10)
    Sk = np.random.randn(10)
    x = np.random.randn(n)

    # Repeat 10 times for more accurate timing
    st_time = time.time()
    for _ in range(10):
        y = fast_truncated_svd_matvec(Uk, Sk, Vk, x)
    timing[n_idx] =  (time.time() - st_time) / 10

# Scaling plot
plt.figure()
plt.loglog(ns, timing, label = 'timing')
scale =   timing[-1] / ns[-1]
plt.loglog(ns, ns * scale, label = '$O(n)$')
plt.xlabel('n')
plt.legend()

error in truncated SVD mat-vec: 1.0

<matplotlib.legend.Legend at 0x136291970>

import scipy.linalg

def circulant_matvec(c,x):
    ###ADD YOUR CODE HERE
    return np.zeros_like(x)


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

# Testing for correctness
n = 10
c = np.random.randn(n)
C = scipy.linalg.circulant(c) # Forms the n \times n circulant matrix
x = np.random.randn(n)
# C@x is an O(n^2) algorithm!
error = np.linalg.norm(C@x - circulant_matvec(c,x) )
print('error in circulant mat-vec:', error)


# Testing complexity
ns = np.logspace(3,5, 30).astype(int)
timing = np.zeros(30)
for n_idx,n in enumerate(ns):
    c = np.random.randn(n)
    x = np.random.randn(n)
    st_time = time.time()
    for _ in range(10):
        y = circulant_matvec(c,x)
    timing[n_idx] =  (time.time() - st_time) / 10

# Scaling plot
plt.figure()
plt.loglog(ns, timing, label = 'timing')
scale =   timing[-1] / (ns[-1] * np.log(ns[-1]))
# If your implementation is correct, the two lines should be roughly parallel (up to some noise)
plt.loglog(ns, ns * np.log(ns) * scale, label = '$O(n \log(n))$')
plt.legend()

error in circulant mat-vec: 5.885244950849171

<>:36: SyntaxWarning: invalid escape sequence '\l'
<>:36: SyntaxWarning: invalid escape sequence '\l'
/var/folders/0b/bgs7xkxj0cj4jc50c5v9hvwc0000gq/T/ipykernel_11650/2541319750.py:36: SyntaxWarning: invalid escape sequence '\l'
  plt.loglog(ns, ns * np.log(ns) * scale, label = '$O(n \log(n))$')

<matplotlib.legend.Legend at 0x1364c01a0>

def toeplitz_matvec(c,r,x):
    ###ADD YOUR CODE HERE
    return np.zeros_like(x)

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

# Testing for correctness
n = 10
c = np.random.randn(n) # First column of T
r = np.random.randn(n) # First row of T. Must have r[0] == c[0] (since they are the same entry)
r[0] = c[0]
T = scipy.linalg.toeplitz(c,r) # Forms the n \times n Toeplitz matrix
x = np.random.randn(n)
# Note that T@x is a O(n^2) algorithm once more
error = np.linalg.norm(T@x - toeplitz_matvec(c,r,x) )
print('error in Toeplitz mat-vec:', error)


# Testing complexity 
ns = np.logspace(3,5, 30).astype(int)
timing = np.zeros(30)
for n_idx,n in enumerate(ns):
    c = np.random.randn(n)
    r = np.random.randn(n)
    r[0] = c[0]
    x = np.random.randn(n)
    st_time = time.time()
    for _ in range(10):
        y = toeplitz_matvec(c,r,x)
    timing[n_idx] =  (time.time() - st_time) / 10

# Scaling plot
plt.figure()
plt.loglog(ns, timing, label = 'timing')
scale =   timing[-1] / (ns[-1] * np.log(ns[-1]))
plt.loglog(ns, ns * np.log(ns) * scale, label = '$O(n log(n))$')
plt.legend()

error in Toeplitz mat-vec: 5.191177014108935

<matplotlib.legend.Legend at 0x1571e69f0>

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

def hankel_matvec(c, r, x):
    """
    Return y = H x where H = hankel(c, r), with r[0] == c[-1].
    Target complexity: O(n log n).
    """
    ### ADD YOUR CODE HERE
    return np.zeros_like(x)


# Testing for correctness
n = 10
c = np.random.randn(n)           # first column of H
r = np.random.randn(n)           # last row of H
r[0] = c[-1]                     # Hankel consistency: overlap
H = scipy.linalg.hankel(c, r)    # forms the n x n Hankel matrix
x = np.random.randn(n)

err = np.linalg.norm(H @ x - hankel_matvec(c, r, x))
print('error in Hankel mat-vec:', err)

# Testing complexity
ns = np.logspace(3, 5, 30).astype(int)
timing = np.zeros_like(ns, dtype=float)
for k, n in enumerate(ns):
    c = np.random.randn(n)
    r = np.random.randn(n); r[0] = c[-1]
    x = np.random.randn(n)
    t0 = time.time()
    for _ in range(10):
        y = hankel_matvec(c, r, x)
    timing[k] = (time.time() - t0) / 10

# Scaling plot
plt.figure()
plt.loglog(ns, timing, label='timing')
scale = timing[-1] / (ns[-1] * np.log(ns[-1]))
plt.loglog(ns, ns * np.log(ns) * scale, label=r'$O(n \log n)$')
plt.xlabel('n'); plt.ylabel('time (s)'); plt.legend(); plt.title('Hankel mat-vec scaling')
plt.show()

error in Hankel mat-vec: 7.314123156776811

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

def fast_pairwise_distance(X,Y):
    # Computes the Euclidean pairwise distance between the columns of X and Y 
    # using high level BLAS operations

    ###ADD YOUR CODE HERE
    return np.zeros((X.shape[1], Y.shape[1]))



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

def naive_pairwise_distance(X,Y):
    # Computes the Euclidean pairwise distance between the columns of X and Y the slow way
    B = np.zeros((X.shape[1], Y.shape[1]))
    for i in range(X.shape[1]):
        for j in range(Y.shape[1]):
            B[i,j] = np.linalg.norm(X[:,i] - Y[:,j])**2
    return B


# Testing for correctness
n, l, m = 11, 12, 13
X = np.random.randn(m,n)
Y = np.random.randn(m,l)

error = np.linalg.norm(fast_pairwise_distance(X,Y) - naive_pairwise_distance(X,Y) ) / np.linalg.norm(naive_pairwise_distance(X,Y))
print('relative error in pairwise distance:', error)


# Testing for speed
n, l, m = 1000, 500, 900
X = np.random.randn(m,n)
Y = np.random.randn(m,l)

st_time = time.time()
naive_pairwise_distance(X,Y) 
naive_time = time.time() - st_time
print('time for naive computation:' , naive_time)

st_time = time.time()
fast_pairwise_distance(X,Y) 
fast_time = time.time() - st_time
print('time for fast computation:', fast_time)

print('Fast computation is :', naive_time / fast_time, ' faster!')

relative error in pairwise distance: 1.0

time for naive computation: 1.3081071376800537
time for fast computation: 0.0004088878631591797
Fast computation is : 3199.1830903790087  faster!

def fast_translation_computation(y1, y2):
    # Implement a O(n log n) method to compute 
    # the norm squared difference of the shifts of y1 and y2.
    
    ###ADD YOUR CODE HERE

    return np.zeros(y1.size)


### Code below is provided for you. No need to change it
def naive_translation_computation(y1, y2):
    m = y1.size
    B = np.zeros(m)
    for k in range(m):
        B[k] = np.linalg.norm(np.roll(y1, k) - y2)**2
    return B


# Testing for correctness
m = 10
y1 = np.random.randn(m)
y2 = np.random.randn(m)
error = np.linalg.norm(naive_translation_computation(y1,y2) - fast_translation_computation(y1,y2) ) / np.linalg.norm(naive_translation_computation(y1,y2))
print('absolute error in shifted distances:', error)


# Testing complexity
ns = np.logspace(3,5, 30).astype(int)
timing = np.zeros(30)
for n_idx,n in enumerate(ns):
    y1 = np.random.randn(n)
    y2 = np.random.randn(n)
    st_time = time.time()
    for _ in range(10):
        fast_translation_computation(y1,y2)
    timing[n_idx] =  time.time() - st_time

# Scaling plot
plt.figure()
plt.loglog(ns, timing, label = 'timing')
scale =   timing[-1] / (ns[-1] * np.log(ns[-1]))
plt.loglog(ns, ns * np.log(ns) * scale, label = '$O(n \log(n))$')
plt.legend()

absolute error in shifted distances: 1.0

<>:42: SyntaxWarning: invalid escape sequence '\l'
<>:42: SyntaxWarning: invalid escape sequence '\l'
/var/folders/0b/bgs7xkxj0cj4jc50c5v9hvwc0000gq/T/ipykernel_11650/90526894.py:42: SyntaxWarning: invalid escape sequence '\l'
  plt.loglog(ns, ns * np.log(ns) * scale, label = '$O(n \log(n))$')

<matplotlib.legend.Legend at 0x16a1ea0f0>

import numpy as np
## These is no code for you to write in this cell. Just run it.
def givens_rotation(x):
    assert len(x) == 2, "Input vector must be of length 2"
    r = np.hypot(x[0], x[1])
    # np.hypot is more stable version of the following line
    # r = np.sqrt(x[0]**2 + x[1]**2)

    if r == 0:
        return np.eye(2)

    cos_theta = x[0] / r
    sin_theta = -x[1] / r
    G = np.zeros((2,2))
    G[0,0] = cos_theta
    G[0,1] = -sin_theta
    G[1,0] = sin_theta
    G[1,1] = cos_theta
    return G

# Example usage:
x = np.array([3, 4])
G = givens_rotation(x)
print(G)
print(G @ x)

[[ 0.6  0.8]
 [-0.8  0.6]]
[ 5.0000000e+00 -4.4408921e-16]

import numpy as np
import time
import matplotlib.pyplot as plt
def hessenberg_qr(H):
    ## ADD YOUR CODE HERE
    n = H.shape[0]
    Q = np.eye(n)
    R = H.copy()

    for i in range(n-1):
        
        ## ADD YOUR CODE HERE
        # placeholder
        n = n

    return Q, R


# Test for correctness
H = np.triu(np.random.rand(20), -1)
Q, R = hessenberg_qr(H)
print("Is Q orthogonal?:", np.allclose(Q.T @ Q, np.eye(Q.shape[0])))
print("Is R upper triangular?:", np.allclose(R, np.triu(R)))
print("Does Q @ R equal H?:", np.allclose(Q @ R, H))


# Test for complexity and scaling plot
# But for debugging, you can start with smaller ns
ns = [50, 100, 200, 400]

# Once you think you got the right complexity, try on this large example to make sure it has the correct scaling (if your algorithm is in fact n^3, it will take a very long time to run)
# This takes ~1s on my machine. If it takes much longer, you might not have the right complexity.
# ns = [50, 100, 200, 400, 800, 1600, 3200]
times = []

for n in ns:
    H_large = np.triu(np.random.randn(n, n), -1)  # Generate a random Hessenberg matrix
    start_time = time.time()
    Q_large, R_large = hessenberg_qr(H_large)
    end_time = time.time()
    times.append(end_time - start_time)
    print(f"Time taken for n={n}: {end_time - start_time:.4f} seconds")

# 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.xlabel('Matrix size (n)')
plt.ylabel('Time (seconds)')
plt.title('Scaling of Hessenberg QR Decomposition')
plt.legend()
plt.grid(True)
plt.show()


## Also observe the interesting sparsity pattern of the Q: matrix:
# Can you explain it?
plt.imshow(np.abs(Q)>1e-14)
plt.title("Sparsity of factor Q")

Is Q orthogonal?: True
Is R upper triangular?: False
Does Q @ R equal H?: True
Time taken for n=50: 0.0000 seconds
Time taken for n=100: 0.0000 seconds
Time taken for n=200: 0.0003 seconds
Time taken for n=400: 0.0001 seconds

Text(0.5, 1.0, 'Sparsity of factor Q')

HW2: MAT 321 - Numerical Analysis and Scientific Computing.¶

ADD YOUR CODE HERE¶

Problem 1: Poisson's equation¶

1.1¶

Write your solutions here¶

1.2¶

1.3¶

Problem 2: Fast matvecs¶

2.1: truncated SVD.¶

2.2: Circulant matrices¶

2.3: Toeplitz matrices¶

2.4: Hankel matrices¶

Problem 3: Always a BLAS¶

3.1:¶

Write your answer to part 3.1 in this cell (use math and words, NOT code!)¶

3.2:¶

3.3: Fast image alignment¶

Q4: Givens rotations¶

Q4.1 Give a brief explanation of your algorithm and justify its complexity below.¶

Write your answer here¶

4.2 Implement your algorithm in the cell below¶

Problem 5: SVD and Eigenvalue decomposition¶

Write your answer for problem 5 in this cell. You may also turn in hand-written answers for this problem.¶