PS: 这次作业还使用到了Scipy库.
PPS: 由于作业使用的英文, 本次作业我也选择使用英文完成
Generate a n*m matrix A, with random Gaussian entries, and a m*m matrix B, a Toeplitz matrix, where n = 200, m = 500.
Preparation
Code:
import numpy as np
from scipy.linalg import toeplitz
SIZE_N = 200
SIZE_M = 500
MU = 0.0
SIGMA = 1.0
# random Gaussian entries produce sequence with expectation MU and variance SIGMA.
# Generate matrix with random Gaussian entries
matrix_a = SIGMA * np.random.randn(SIZE_N, SIZE_M) + MU
# Generate toeplitz matrix
matrix_b = toeplitz(list(range(1, SIZE_M + 1))[::-1], list(range(1, SIZE_M + 1)[::-1]))
Generated Matrix B:
array([[500, 499, 498, ..., 3, 2, 1],
[499, 500, 499, ..., 4, 3, 2],
[498, 499, 500, ..., 5, 4, 3],
...,
[ 3, 4, 5, ..., 500, 499, 498],
[ 2, 3, 4, ..., 499, 500, 499],
[ 1, 2, 3, ..., 498, 499, 500]])
Ex 9-1 Matrix operation
Quite simple. View comments in the following code:
#Ex 9-1
# Compute A+A
result = matrix_a + matrix_a
# compute A*(A)^T
result = np.dot(matrix_a, matrix_a.T)
# compute (A)^T*A
result = np.dot(matrix_a.T, matrix_a)
# compute A*B
result = np.dot(matrix_a, matrix_b)
# given lambda, compute A *(B - lambda * I)
# Attention! lambda is a keyword in python
def myfun1(A, B, p):
return np.dot(A, (B - p * np.identity(B.shape[0])))
Ex 9-2 Solving a linear system
Solve linear systems of equation by the inverse of coefficient matrix:
# Generate a random vector b
vector_b = np.random.random((matrix_b.shape[1]))
# x = B^(-1) * b
solution = np.linalg.inv(matrix_b).dot(vector_b)
Or you can use function numpy.linalg.solve directly.
solution = np.linalg.solve(matrix_b, vector_b)
Here is details for numpy.linalg.solve.
Ex 9-3 Norm
Also quite simple, view comments in the code below:
# Compute the Frobenius norm of A
# Notice that the 2nd param can be None, for it's set to get Frobenius norm for a matrix as default
frobenius_norm_A = np.linalg.norm(matrix_a, 'fro')
# Compute the infinity norm of B
# 2nd param cannot be ignored
infinity_norm_B = np.linalg.norm(matrix_b, np.inf)
# Get the maximum & minimum of B
max_element_of_B = matrix_b.max()
min_element_of_B = matrix_b.min()
Ex 9-4 Power iteration
View Here for the explanation of Power iteration.
Here I write a function power_iterationt(A, error) which computes the approximate largest eigenvalues of a given matrix A. The max error of eigenvalues will not exceed error. The function will return the largest eigenvalue, the corresponding eigenvector of A, and it's iteration times.
# function to compute eigenvalue with Power Iteration
def power_iteration(A, error):
times = 0 #Iteration count
last_eig = 1.0
new_eig = 0.0
u = np.random.rand(A.shape[0])
print(isinstance(last_eig, float))
while (times == 0 or abs(last_eig - new_eig) > error): #stopping condition
last_eig = new_eig
v = A.dot(u)
new_eig = abs(v.max())
u = v / new_eig #normalize the vector, it's a key part for power iteration.
times = times + 1
return new_eig, u, times
And I test the iteration times given different matrix sizes(10, 100, 500, 1000, 2000, 10000):
error = 0.0001
for i in [10, 100, 500, 1000, 2000, 10000]:
A = np.random.rand(i, i)
print("SIZE = %d, Iteration times = %d" % (i, power_iteration(A, error)[2]))
One of the test results are:
SIZE = 10, Iteration times = 8
SIZE = 100, Iteration times = 6
SIZE = 500, Iteration times = 6
SIZE = 1000, Iteration times = 6
SIZE = 2000, Iteration times = 5
SIZE = 10000, Iteration times = 6
No apparent differences on iteration timesare shown with different size of matrix.
And I use time.clock() to compare computation time when given different size of matrix:
error = 0.0001
for i in [10, 100, 500, 1000, 2000, 10000]:
A = np.random.rand(i, i)
start = time.clock()
power_iteration(A, error)
end = time.clock()
print("SIZE = %d, computation time = %fs" % (i, end - start))
One of the test results is:
SIZE = 10, computation time = 0.000109s
SIZE = 100, computation time = 0.000394s
SIZE = 500, computation time = 0.008850s
SIZE = 1000, computation time = 0.025489s
SIZE = 2000, computation time = 0.057395s
SIZE = 10000, computation time = 1.431594s
It seems that the running time is not proportional to the matrix size. The relation between them needs more analysis.
Ex 9-5
# Size of Matrix C
n = 10
def random_zero_or_one(p):
return (float)(np.random.random() > 1-p)
matrix_c = np.array(list(map(random_zero_or_one, [0.5] * n * n))).reshape((n, n))
# SVD factorization
U, s, Vh = scipy.linalg.svd(matrix_c)
print(s)
And I test the singular value of C given different n and p:
for n in list(range(2, 100, 5)):
average = 0
p = np.random.random()
for test_times in range(10):
matrix_c = np.array(list(map(random_zero_or_one, [p] * n * n))).reshape((n, n))
U, s, Vh = scipy.linalg.svd(matrix_c)
average = average +s.max()
average = average / 10
print("singular value: %f, p = %f, n = %d, p*n = %f" % (average, p, n, p*n))
Results:
singular value: 1.123607, p = 0.355978, n = 2, p*n = 0.711957
singular value: 2.318768, p = 0.231591, n = 7, p*n = 1.621135
singular value: 2.286219, p = 0.112528, n = 12, p*n = 1.350333
singular value: 16.922299, p = 0.995744, n = 17, p*n = 16.927640
singular value: 12.941962, p = 0.573954, n = 22, p*n = 12.626985
singular value: 17.642513, p = 0.646745, n = 27, p*n = 17.462121
singular value: 25.088633, p = 0.777588, n = 32, p*n = 24.882804
singular value: 21.929907, p = 0.581553, n = 37, p*n = 21.517446
singular value: 33.011162, p = 0.781029, n = 42, p*n = 32.803204
singular value: 13.905496, p = 0.279924, n = 47, p*n = 13.156423
singular value: 13.855356, p = 0.250213, n = 52, p*n = 13.011097
singular value: 55.627474, p = 0.974140, n = 57, p*n = 55.525976
singular value: 3.424217, p = 0.036468, n = 62, p*n = 2.261025
singular value: 64.438280, p = 0.960761, n = 67, p*n = 64.370955
singular value: 35.216698, p = 0.477921, n = 72, p*n = 34.410340
singular value: 61.532490, p = 0.794089, n = 77, p*n = 61.144878
singular value: 23.501740, p = 0.276596, n = 82, p*n = 22.680886
singular value: 26.268507, p = 0.294773, n = 87, p*n = 25.645271
singular value: 55.765074, p = 0.603855, n = 92, p*n = 55.554649
singular value: 4.250828, p = 0.031018, n = 97, p*n = 3.008710
The max singular value may be the product of p and n
Ex 9-6
Method:
First construct a array B with the same size of A, where each entries are z.
Then find the minimum absolute value of B - A with function numpy.argmin, which return the index.
Lastly index the entry in A.
def get_nearest_neighbor(z, a):
# print info of z and A
print("z = %f" % z)
print("A = ")
print(a)
# return the answer
return a[np.argmin(abs(a-z))]
# test code
# randomly generate z and A
z = np.random.rand()
a = np.random.rand(10)
print(get_nearest_neighbor(z, a))
