Numpy Exercise

Numpy Exercise

Generate matrices A, with random Gaussian entries, B, a Toeplitz matrix, where A ∈ Rn×m and B ∈ Rm×m, for n = 200, m = 500.

Exercise 9.1: Matrix operations

Calculate A + A, AA⊤, A⊤A and AB. Write a function that computes A(B − λI) for any λ.

from scipy.linalg import toeplitz
import numpy as np
import time
n = 2
m = 4

A = np.random.normal(size=(n, m))
B = toeplitz(range(m))

print(A+A)
print(np.matmul(A, A.T))
print(np.matmul(A.T, A))
print(np.matmul(A, B))


def compute_with_constant(A, B, constant):
    return np.matmul(A, (B - constant*np.eye(m)))

Exercise 9.2: Solving a linear system

Generate a vector b with m entries and solve Bx = b.

b = np.random.random(m)


def solve_linear(B, b):
    return np.linalg.solve(B, b)

Exercise 9.3: Norms

Compute the Frobenius norm of A: ∥A∥F and the infinity norm of B: ∥B∥∞. Also find the largest and smallest singular values of B.

print(np.linalg.norm(A, 'fro'))
print(np.linalg.norm(B, np.inf))
U, S, VT = np.linalg.svd(B)
print(max(S))
print(min(S))

Exercise 9.4: Power iteration

Generate a matrix Z, n × n, with Gaussian entries, and use the power iteration to find the largest eigenvalue and corresponding eigenvector of Z. How many iterations are needed till convergence?
Optional: use the time.clock() method to compare computation time when varying n.

Z = np.random.normal(size=(n, n))
# eigvals, eigvts = np.linalg.eig(Z)
# print(max(eigvals))
stop = 1e-8


def power_iteration(Z, stop):
    u = np.random.normal(size=[len(Z)]).transpose()
    u = u/np.max(u)
    v = u.copy()
    eigen_change = np.inf
    eigen_val = np.max(v)

    iteration = 0
    s_time = time.clock()
    while abs(eigen_change) > stop:
        v = np.matmul(Z, u)
        max_v = np.max(v)
        eigen_change = eigen_val - max_v
        eigen_val = max_v
        u = v/max_v
        iteration += 1
    print("The result is {}".format(np.allclose(np.matmul(Z, u), eigen_val*u)))
    print("Total iteration: {}".format(iteration))
    print("Total time: {}s".format(time.clock() - s_time))

Exercise 9.5: Singular values

Generate an n × n matrix, denoted by C, where each entry is 1 with probability p and 0 otherwise. Use the linear algebra library of Scipy to compute the singular values of C. What can you say about the relationship between n, p and the largest singular value?

p_list = np.arange(0, 1, 0.1)

for p in p_list:
    C = np.random.random((n, n))
    for i in range(n):
        for j in range(n):
            C[i][j] = 1 if C[i][j] > p else 0

    U, S, VT = np.linalg.svd(C)
    max_singular = np.max(S)
    print("p = {},   max_singular = {}".format(p, max_singular))

Exercise 9.6: Nearest neighbor

Write a function that takes a value z and an array A and finds the element in A that is closest to z. The function should return the closest value, not index.
Hint: Use the built-in functionality of Numpy rather than writing code to find this value manually. In particular, use brackets and argmin.

def find_closest(A, z):
    index = np.argmin(np.abs(A-z))
    return A[index]