numpy

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

import numpy as np
from scipy.linalg import toeplitz
A = np.random.normal(2, 1, [200, 500])
B = toeplitz(np.random.normal(2, 1, 500), np.random.normal(2, 1, 500))

Exercise 9.1

Matrix operations Calculate A + A, AAT, ATA and AB. Write a function that computes A(B−λI) for any λ.

print("A:", A)
print("B:", B)
print("A + A:", A + A)
print("AAT:", np.dot(A, A.T))
print("ATA:", np.dot(A.T, A))
print("AB:", np.dot(A, B))

def foo(λ):
return np.dot(A, B - λ)

A: [[1.65428653 1.29250801 0.12205231 0.74224551 2.11893814]
[1.31702596 1.57987211 1.76251052 1.09167503 0.72512259]]
B: [[1.75756406 2.38712337 0.4477556  3.09135636 1.92525955]
[2.42381625 1.75756406 2.38712337 0.4477556  3.09135636]
[2.31593935 2.42381625 1.75756406 2.38712337 0.4477556 ]
[2.61638874 2.31593935 2.42381625 1.75756406 2.38712337]
[2.29753941 2.61638874 2.31593935 2.42381625 1.75756406]]
A + A: [[3.30857306 2.58501602 0.24410461 1.48449102 4.23787628]
[2.63405192 3.15974422 3.52502104 2.18335006 1.45024519]]
AAT: [[9.46296488 6.78263494]
[6.78263494 9.05455374]]
ATA: [[4.47122129 4.21891117 2.52318159 2.6656511  4.4603361 ]
[4.21891117 4.16657284 2.94229479 2.6840652  3.88434548]
[2.52318159 2.94229479 3.12134009 2.0146815  1.53665748]
[2.6656511  2.6840652  2.0146815  1.74268277 2.36437055]
[4.4603361  3.88434548 1.53665748 2.36437055 5.01570162]]
AB: [[13.13332887 13.77944548 10.74700582 12.42453151 12.73118452]
[14.74818863 14.61808725 11.78414015 12.66237593 12.08915085]]

Exercise 9.2

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

print("B:")
print(B)

b = np.random.random(5)
print("b:", b)
x = np.linalg.solve(B, b)
print("x:", x)
print("np.dot(B, x):\n", np.dot(B, x))

B:
[[ 2.5121108   1.84558016  1.69360915  1.14977797  3.51985693]
[ 2.89274206  2.5121108   1.84558016  1.69360915  1.14977797]
[ 3.59327592  2.89274206  2.5121108   1.84558016  1.69360915]
[ 0.87142587  3.59327592  2.89274206  2.5121108   1.84558016]
[-0.09840123  0.87142587  3.59327592  2.89274206  2.5121108 ]]
b: [0.63848187 0.88284346 0.00751235 0.26442866 0.93755242]
x: [ 0.23215661 -0.41407927 -2.29703059  2.99870062  0.35851666]
np.dot(B, x):
[0.63848187 0.88284346 0.00751235 0.26442866 0.93755242]

Exercise 9.3

Norms Compute the Frobenius norm of A: ||A||F and the inﬁnity norm of B: ||B||∞. Also ﬁnd the largest and smallest singular values of B.

print("A:")
print(A)
print("B:")
print(B)

print("the Frobenius norm of A:", np.linalg.norm(A))
print("the infinity norm of B:", np.linalg.norm(B, np.inf))
print("the largest singular value of B:", np.linalg.norm(B, 2))
print("the smallest singular value of B:", np.linalg.norm(B, -2))

A:
[[ 0.89763757  3.90870489  2.24489208  1.06493451 -0.38395911]
[ 1.9064736  -0.60809743  2.44124628  3.29251693  2.95935474]]
B:
[[3.26704861 1.56438763 1.90880029 3.5782664  2.64795173]
[2.55789806 3.26704861 1.56438763 1.90880029 3.5782664 ]
[3.60588272 2.55789806 3.26704861 1.56438763 1.90880029]
[1.30917891 3.60588272 2.55789806 3.26704861 1.56438763]
[0.96130538 1.30917891 3.60588272 2.55789806 3.26704861]]
the Frobenius norm of A: 7.208837168639553
the infinity norm of B: 12.966454652614956
the largest singular value of B: 12.56650752516994
the smallest singular value of B: 1.8246744060618487

Exercise 9.4

Power iteration Generate a matrix Z, n × n, with Gaussian entries, and use the power iteration to ﬁnd 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.

import numpy as np
import time
def foo(n):
Z = np.random.normal(0, 1, [n, n])
b_k = np.random.rand(Z.shape[0])
b_k1 = np.random.rand(Z.shape[0])
while b_k.all() != b_k1.all():
b_k1 = np.dot(Z, b_k)
b_k1_norm = np.linalg.norm(b_k1)
b_k = b_k1 / b_k1_norm
# print(b_k)
# print(np.dot(Z, b_k))
print("n:", n, "time:", time.clock())
t = time.clock()
foo(5)
foo(10)
foo(100)
foo(1000)
foo(10000)

n: 5 time: 0.00014279107049498438
n: 10 time: 0.0002451910413678816
n: 100 time: 0.001004088603281464
n: 1000 time: 0.07183528623351858
n: 10000 time: 5.034103741410491

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?

import numpy as np
import scipy as sp
# print(help(sp.linalg.eigvals))
def foo(p, n):
C = np.random.binomial(1, p, [n, n])
# print(sp.linalg.eig(C))
print("p:", p, "n:", n, max(sp.linalg.eigvals(C)))

foo(.5, 10)
foo(.5, 20)
foo(.5, 100)
foo(.3, 10)
foo(.3, 20)
foo(.3, 100)
foo(.9, 10)
foo(.9, 20)
foo(.9, 100)

p: 0.5 n: 10 (5.641844752967631+0j)
p: 0.5 n: 20 (9.242385892004453+0j)
p: 0.5 n: 100 (50.14720097026784+0j)
p: 0.3 n: 10 (3.456889014392734+0j)
p: 0.3 n: 20 (6.498898769869614+0j)
p: 0.3 n: 100 (30.05949038152027+0j)
p: 0.9 n: 10 (9.010646243839297+0j)
p: 0.9 n: 20 (17.458772002141625+0j)
p: 0.9 n: 100 (90.72610340531895+0j)

$\underset{n\to \infty }{lim}max\left(eigvals\right)=np$

Exercise 9.6

Nearest neighbor Write a function that takes a value z and an array A and ﬁnds 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 ﬁnd this value manually. In particular, use brackets and argmin.

import numpy as np

print(help(np.argmin))
def foo(A, z):
B = np.abs(A - z)
return A[np.argmin(B)]

a = np.random.randint(0, 100, 6)
print(a)
b = np.random.randint(100)
print("b:", b, "Nearest neighbor:", foo(a,b))

[17 93 87 77 20 26]
b: 9 Nearest neighbor: 17