Numpy

Numpy Exercises

注：为方便测试，在测试过程中，以下程序生成的矩阵均采用较小的维度

Exercise 9.1: Matrix operations

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

n = 200
m = 500

def MyPrint(ans):
    print (ans)
    print ("")

def CalculateAnswer(A, B):
    """Do calculation"""
    ans1 = A + A
    print ("A + A")
    MyPrint(ans1)

    ans2 = A.dot(A.T)
    print ("A * A.T")
    MyPrint(ans2)

    ans3 = (A.T).dot(A)
    print ("A.T * A")
    MyPrint(ans3)

    ans4 = A.dot(B)
    print ("A * B")
    MyPrint(ans4)

    I = np.identity(m)
    k = int(input("Input any λ for A(B − λI)\n"))
    ans5 = A.dot(B - k * I)
    print ("A * (B - kI)")
    MyPrint(ans5)

# Generate matrix A with random Gaussian entries
A = np.random.normal(0, 1, n * m).reshape(n, m)
print ("Matrix A")
MyPrint(A)

# Generate Toeplitz matrix B
B_data_col = []
B_data_row = []
for i in range(0, m):
    rand = random.randint(1, m)
    B_data_col.append(rand)
    rand = random.randint(1, m)
    B_data_row.append(rand)
B = toeplitz(B_data_col, B_data_row)
print ("Matrix B")
MyPrint(B)

# Another method to generate toeplitz matrix
# B = np.zeros((m, m))
# i, j = np.indices(B.shape)
# for index in range(-m + 1, m):
#   B[i == j + index] = random.randint(1, m)
# print ("Matrix B")
# MyPrint(B)

CalculateAnswer(A, B)

Test Exercise 9.1

Exercise 9.2: Solving a linear system

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

m = 500

def MyPrint(ans):
    print (ans)
    print ("")

def CalculateAnswer(B):
    """Bx = b"""
    arr = []
    for i in range(0, m):
        rand = random.randint(1, m)
        arr.append(rand)
    b = np.array(arr)
    print ("b")
    MyPrint (b)

    x = np.linalg.solve(B, b)
    print ("x")
    MyPrint (x)

# Generate Toeplitz matrix B
B_data_col = []
B_data_row = []
for i in range(0, m):
    rand = random.randint(1, m)
    B_data_col.append(rand)
    rand = random.randint(1, m)
    B_data_row.append(rand)
B = toeplitz(B_data_col, B_data_row)
print ("Matrix B")
MyPrint(B)

CalculateAnswer(B)

Test Exercise 9.2

Exercise 9.3: Norms

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

n = 200
m = 500

def MyPrint(ans):
    print (ans)
    print ("")

def CalculateAnswer(A ,B):
    """Compute the Frobenius norm of A and the infinity norm of B"""
    print ("the Frobenius norm of A")
    MyPrint(np.linalg.norm(A))

    print ("the Infinity norm of B")
    MyPrint(np.linalg.norm(B, np.inf))

    print ("largest and smallest singular values of B")
    U, sigma, VT = np.linalg.svd(B)
    print ("largest value", max(sigma))
    print ("smallest value", min(sigma))

# Generate matrix A with random Gaussian entries
A = np.random.normal(0, 1, n * m).reshape(n, m)
print ("Matrix A")
MyPrint(A)

# Generate Toeplitz matrix B
B_data_col = []
B_data_row = []
for i in range(0, m):
    rand = random.randint(1, m)
    B_data_col.append(rand)
    rand = random.randint(1, m)
    B_data_row.append(rand)
B = toeplitz(B_data_col, B_data_row)
print ("Matrix B")
MyPrint(B)

CalculateAnswer(A, B)

Test Exercise 9.3

Exercise 9.4: Power iteration

import numpy as np
import time

n = 3

def eigenvalue(Z, v):
    Zv = Z.dot(v)
    return v.dot(Zv)

def power_iteration(Z):
    n, d = Z.shape
    v = np.ones(d) / np.sqrt(d)
    ev = eigenvalue(Z, v)
    iteration_count = 0

    while True:
        iteration_count += 1
        Zv = Z.dot(v)
        v_new = Zv / np.linalg.norm(Zv)
        ev_new = eigenvalue(Z, v_new)
        if np.abs(ev - ev_new) < 0.01:
            break
        v = v_new
        ev = ev_new

    return ev_new, v_new, iteration_count

Z = np.random.normal(10, 10, n * n).reshape(n, n)
print ("Matrix Z")
print (Z)
print ("")

t_prev = time.clock()
eigenvalue, eigenvector, iteration_count = power_iteration(Z)
t = time.clock()

print ("Eigenvalue:", eigenvalue)
print ("Eigenvector:", eigenvector)
print ("Iteration count:", iteration_count)
print ("Time cost:", t - t_prev)

Test Exercise 9.4

Exercise 9.5: Singular values

import numpy as np
from scipy import linalg
import random

n = 200
p = 1 / 3

def MyPrint(ans):
    print (ans)
    print ("")

def random_pick():
    """Generate a random number dependent on p"""
    x = random.uniform(0, 1)
    cumulative_probability = 0.0
    for item, item_probability in zip([1, 0], [p, 1 - p]):
        cumulative_probability += item_probability
        if x < cumulative_probability: break
    return item

def CalculateAnswer(C):
    """Compute the singular values of C"""
    U, s, Vh = linalg.svd(C)
    print ("singular values")
    MyPrint(s)
    print ("max singular value", max(s))
    print ("Conclusion: max singular value nearly equals to n * p")

# Generate matrix C with entry 1 and 0
C_data = []
for i in range(0, n * n):
    C_data.append(random_pick())
C = np.reshape(C_data, (n, n))
print ("Matric C")
MyPrint(C)

CalculateAnswer(C)

Test Exercise 9.5

Exercise 9.6: Nearest neighbor

import numpy as np
import random

n = 100         # 数组大小100
m = 500         # 数组元素范围1-500

def Myprint(ans):
    print (ans)
    print ("")

A = []
for i in range(0, n):
    A_data = random.randint(1, m)
    A.append(A_data)
A = np.array(A)
print ("Array A")
Myprint(A)

z = random.randint(1, m)
print ("Value z =", z)
Z = np.array([z for i in range(0, n)])
A_temp = abs(A - Z)

print ("The element in A that is closet to z:", A[np.argmin(A_temp)])

Test Exercise 9.6