Exercise: 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.

構造A：

利用內建函式numpy.random. normal(loc=0.0, scale=1.0, size=None)

numpy.random.normal(loc=0.0, scale=1.0, size=None)

Draw random samples from a normal(Gaussian) distribution.

The probability density function of thenormal distribution, first derived by De Moivre and 200 years later by bothGauss and Laplace independently [R500500], is often called the bell curvebecause of its characteristic shape (see the example below).

The normal distributions occurs often innature. For example, it describes the commonly occurring distribution ofsamples influenced by a large number of tiny, random disturbances, each withits own unique distribution [R500500].

Parameters:

loc : float or array_like of floats Mean (“centre”) of the distribution. scale : float or array_like of floats Standard deviation (spread or “width”) of the distribution. size : int or tuple of ints, optional Output shape. If the given shape is, e.g., (m, n, k), then m * n * k samples are drawn. If size is None (default), a single value is returned if loc and scale are both scalars. Otherwise, np.broadcast(loc, scale).size samples are drawn.

Returns:

out : ndarray or scalar Drawn samples from the parameterized normal distribution.

scipy.linalg.toeplitz(c, r=None)[source]

Construct a Toeplitz matrix.

The Toeplitz matrix has constant diagonals,with c as its first column and r as its first row. If r is not given, r == conjugate(c)is assumed.

Parameters:

c : array_like First column of the matrix. Whatever the actual shape of c, it will be converted to a 1-D array. r : array_like, optional First row of the matrix. If None, r = conjugate(c) is assumed; in this case, if c[0] is real, the result is a Hermitian matrix. r[0] is ignored; the first row of the returned matrix is [c[0], r[1:]]. Whatever the actual shape of r, it will be converted to a 1-D array.

Returns:

A : (len(c), len(r)) ndarray The Toeplitz matrix. Dtype is the same as (c[0] + r[0]).dtype.

程式碼：

import numpy as np
from scipy.linalg import toeplitz
import time
mu, sigma = 0, 1.0 # mean and standard deviation
n, m = 200, 500
A = np.random.normal(loc=mu, scale=sigma, size=(n, m))

#scipy.linalg.toeplitz Used to construct the convolution operator
c = [a for a in range(1, m+1)]
B = toeplitz(c, c)

Exercise 9.1: Matrix op erations
Calculate A + A, AA >, A >A and AB. Write a function that computes A(B − λI) for any λ.

#Exercise 9.1
def Exercise_9_1(A, B, n, m):
	#Calculate A + A
	print("A + A:")
	C = A + A
	print(C)
	
	#這一句支援將結果矩陣儲存到txt檔案功能
	#按需要選擇是否啟用
	#C.tofile("A+A.txt", sep=" ", format="%s")
	
	#Calculate AA >
	print("AA^:")
	print(np.dot(A, A.T))
	
	#Calculate A >A
	print("A^A:")
	print(np.dot(A.T, A))
	
	#Calculate AB
	print("AB:")
	print(np.dot(A, B))
	
	#Write a function that computes A(B − λI) for any λ
	fun_lambda(A, B, 2.0);

def fun_lambda(A, B, lamda):
	print("A(B − λI):")
	C = B - lamda * (np.eye(m))
	print(np.dot(A, C))

Exercise 9.2: Solving a linear system

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

#exercise 9.2
def Exercise_9_2(A, B, n, m):
	#Generate a vector b with m entries and solve Bx = b.
	b = np.ones((m, 1))
	x = np.linalg.solve(B, b)
	print(x)

Exercise 9.3: Norms
Compute the Frobenius norm of A: k Ak F and the infinity norm of B: k Bk ∞. Also find the largest and

smallest singular values of B.

#exercise 9.3
def Exercise_9_3(A, B, n, m):
	#Compute the Frobenius norm of A: k Ak F
	A_F = np.linalg.norm(A, 'fro')
	print("the Frobenius norm:", A_F)
	
	#the infinity norm of B: k Bk ∞.
	B_F = np.linalg.norm(B, np.inf)
	print("the infinity norm:", B_F)
		
	#find the largest and smallest singular values of B
	lar_sin = np.linalg.norm(B, 2)
	smal_sin = np.linalg.norm(B, -2)
	print("the largest singular:", lar_sin)
	print("the smallest singular:", smal_sin)

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

#exercise 9.4
def Exercise_9_4(A, B, n, m):
	#Generate a matrix Z, n × n, with Gaussian entries
	Z = np.random.standard_normal((n, n))
	#use the power iteration to find the largest eigenvalue and corresponding eigenvector of Z
	num = 0
	u_k = np.ones(n)
	v_k_norm = 0
	v_k = np.zeros(n)
	
	begin = time.clock()
	while(True):
		# calculate the matrix-by-vector product Ab
		v_k = np.dot(Z, u_k)
		# calculate the norm
		v_k_norm_temp = v_k_norm
		v_k_norm = np.linalg.norm(v_k)
		# re normalize the vector
		u_k = v_k / v_k_norm
		num += 1
		if(abs(v_k_norm_temp - v_k_norm) < 0.0005):
			break;
	end = time.clock()
	
	print("the largest eigenvalue:", v_k_norm)
	print("the corresponding eigenvector:", u_k)
	#How many iterations are needed till convergence
	print("The number of iterations:", num)
	#Optional: use the time.clock() method to compare computation time when varying n.
	print("computation time when varying n:", end-begin)

Exercise 9.5: Singular values
Generate an n × n matrix, denoted by C, where each entry is 1 withprobability 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?

#exercise 9.5
def Exercise_9_5(A, B, n, m):
	#Generate an n × n matrix, denoted by C, where each entry is 1 with probability p and 0 otherwise
	p = 0.5
	C = np.random. binomial(1, p, (n, n))
	
	#Use the linear algebra library of Scipy to compute the singular values of C
	lar_sin = np.linalg.norm(C, 2)
	smal_sin = np.linalg.norm(C, -2)
	print("the smallest singular:", smal_sin)
	print("the largest singular:", lar_sin)
	
	#What can you say about the relationship between n, p and the largest singular value?
	print("n * p:", n*p)
	print("the largest singular is closed with n * p \nso that we can say they are equal!")

Exercise 9.6: Nearest neighb or
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 findthis value manually. In

particular, use brackets and argmin.

#exercise 9.6
def Exercise_9_6(A, B, n, m):
	#Write a function that takes a value z and an array A
	z = -5
	closest = fun_closest(A, z)
	
	#The function should return the closest value, not index.
	print("the closest value:", closest)
	
def fun_closest(A, z):
	#finds the element in A that is closest to z
	B, C = A[A>z], A[A<=z]
	ceil, floor = 0, 0
	
	#Hint: Use the built-in functionality of Numpy rather than writing code to find this value manually
	if(len(B)):
		ceil = np.argmin(B)
	else:
		return C[np.argmax(C)]
	
	#	   In particular, use brackets and argmin.	
	if(len(C)):
		floor = np.argmax(C)
	else:
		return B[ceil]
		
	if(abs(B[ceil]-z) < abs(C[floor]-z)):
		return B[ceil]
	else:
		return C[floor]

完整程式碼：

import numpy as np
from scipy.linalg import toeplitz
import time

#Exercise 9.1
def Exercise_9_1(A, B, n, m):
	#Calculate A + A
	print("A + A:")
	C = A + A
	print(C)
	
	#這一句支援將結果矩陣儲存到txt檔案功能
	#按需要選擇是否啟用
	#C.tofile("A+A.txt", sep=" ", format="%s")
	
	#Calculate AA >
	print("AA^:")
	print(np.dot(A, A.T))
	
	#Calculate A >A
	print("A^A:")
	print(np.dot(A.T, A))
	
	#Calculate AB
	print("AB:")
	print(np.dot(A, B))
	
	#Write a function that computes A(B − λI) for any λ
	fun_lambda(A, B, 2.0);

def fun_lambda(A, B, lamda):
	print("A(B − λI):")
	C = B - lamda * (np.eye(m))
	print(np.dot(A, C))
	
#exercise 9.2
def Exercise_9_2(A, B, n, m):
	#Generate a vector b with m entries and solve Bx = b.
	b = np.ones((m, 1))
	x = np.linalg.solve(B, b)
	print(x)

#exercise 9.3
def Exercise_9_3(A, B, n, m):
	#Compute the Frobenius norm of A: k Ak F
	A_F = np.linalg.norm(A, 'fro')
	print("the Frobenius norm:", A_F)
	
	#the infinity norm of B: k Bk ∞.
	B_F = np.linalg.norm(B, np.inf)
	print("the infinity norm:", B_F)
		
	#find the largest and smallest singular values of B
	lar_sin = np.linalg.norm(B, 2)
	smal_sin = np.linalg.norm(B, -2)
	print("the largest singular:", lar_sin)
	print("the smallest singular:", smal_sin)
	
#exercise 9.4
def Exercise_9_4(A, B, n, m):
	#Generate a matrix Z, n × n, with Gaussian entries
	Z = np.random.standard_normal((n, n))
	#use the power iteration to find the largest eigenvalue and corresponding eigenvector of Z
	num = 0
	u_k = np.ones(n)
	v_k_norm = 0
	v_k = np.zeros(n)
	
	begin = time.clock()
	while(True):
		# calculate the matrix-by-vector product Ab
		v_k = np.dot(Z, u_k)
		# calculate the norm
		v_k_norm_temp = v_k_norm
		v_k_norm = np.linalg.norm(v_k)
		# re normalize the vector
		u_k = v_k / v_k_norm
		num += 1
		if(abs(v_k_norm_temp - v_k_norm) < 0.0005):
			break;
	end = time.clock()
	
	print("the largest eigenvalue:", v_k_norm)
	print("the corresponding eigenvector:", u_k)
	#How many iterations are needed till convergence
	print("The number of iterations:", num)
	#Optional: use the time.clock() method to compare computation time when varying n.
	print("computation time when varying n:", end-begin)

#exercise 9.5
def Exercise_9_5(A, B, n, m):
	#Generate an n × n matrix, denoted by C, where each entry is 1 with probability p and 0 otherwise
	p = 0.5
	C = np.random. binomial(1, p, (n, n))
	
	#Use the linear algebra library of Scipy to compute the singular values of C
	lar_sin = np.linalg.norm(C, 2)
	smal_sin = np.linalg.norm(C, -2)
	print("the smallest singular:", smal_sin)
	print("the largest singular:", lar_sin)
	
	#What can you say about the relationship between n, p and the largest singular value?
	print("n * p:", n*p)
	print("the largest singular is closed with n * p \nso that we can say they are equal!")

#exercise 9.6
def Exercise_9_6(A, B, n, m):
	#Write a function that takes a value z and an array A
	z = -5
	closest = fun_closest(A, z)
	
	#The function should return the closest value, not index.
	print("the closest value:", closest)
	
def fun_closest(A, z):
	#finds the element in A that is closest to z
	B, C = A[A>z], A[A<=z]
	ceil, floor = 0, 0
	
	#Hint: Use the built-in functionality of Numpy rather than writing code to find this value manually
	if(len(B)):
		ceil = np.argmin(B)
	else:
		return C[np.argmax(C)]
	
	#	   In particular, use brackets and argmin.	
	if(len(C)):
		floor = np.argmax(C)
	else:
		return B[ceil]
		
	if(abs(B[ceil]-z) < abs(C[floor]-z)):
		return B[ceil]
	else:
		return C[floor]

mu, sigma = 0, 1.0 # mean and standard deviation
n, m = 200, 500
A = np.random.normal(loc=mu, scale=sigma, size=(n, m))

#scipy.linalg.toeplitz Used to construct the convolution operator
c = [a for a in range(1, m+1)]
B = toeplitz(c, c)

Exercise_9_1(A, B, n, m)
Exercise_9_2(A, B, n, m)
Exercise_9_3(A, B, n, m)
Exercise_9_4(A, B, n, m)
Exercise_9_5(A, B, n, m)
Exercise_9_6(A, B, n, m)