【高级编程技术】【作业】【第十一周】【1】

numpy练习

Generate matrices $A$, with random Gaussian entries, $B$, a Toeplitz matrix, where $A \in \mathbb R^{n \times m}$ and $B \in \mathbb R^{m \times m}$, for $n=200$, $m=500$.

import numpy
import scipy
from scipy.linalg import toeplitz

n, m, mu, sigma = 200, 500, numpy.random.randint(10), numpy.random.randint(10)
A = mu+sigma*numpy.random.randn(n, m)
c_and_r = list(range(1, 2*m))
numpy.random.shuffle(c_and_r)
B = toeplitz(c_and_r[:m], c_and_r[m-1:])

Exercise 9.1: Matrix operations

Calculate $A+A$, $AA^\top$, $A^\top A$ and $AB$. Write a function that computes $A(B-\lambda I)$ for any $\lambda$.

ans1 = A+A
ans2 = numpy.dot(A, A.T)
ans3 = numpy.dot(A.T, A)
ans4 = numpy.dot(A, B)
def f(A, B, Lambda):
    return numpy.dot(A, B-Lambda*numpy.identity(m))

Exercise 9.2: Solving a linear system

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

b = numpy.random.randint(m, size=m)
x = numpy.linalg.solve(B, b)

Exercise 9.3: Norms

Compute the Frobenius norm of $A$: $\Vert A \Vert_F$ and the infinity norm of $B$: $\Vert B \Vert_\infty$. Also find the largest and smallest singular values of $B$.

矩阵范数 - 维基百科，自由的百科全书

forbenius_norm = numpy.linalg.norm(A, 'fro')
infinity_norm = numpy.linalg.norm(B, numpy.inf)
singular_values = numpy.linalg.svd(B)[1]
largest_singular_value = numpy.max(singular_values)
smallest_singualr_value = numpy.min(singular_values)

Exercise 9.4: Power iteration

Generate a matrix $Z$, $n \times 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$.

import time

def Rayleigh_quotient(M, x):
    return numpy.dot(x.T, numpy.dot(M, x))/numpy.dot(x.T, x)

def power_iteration(Z: numpy.ndarray, epsilon):
    b = Zb = numpy.random.rand(Z.shape[0])
    last_b = numpy.zeros(Z.shape[0])
    num_iteration = 0

    while num_iteration==0 or numpy.linalg.norm(b-last_b)>=epsilon:
        last_b = b
        Zb = numpy.dot(Z, b)
        b = Zb/numpy.linalg.norm(Zb)
        num_iteration += 1

    return Rayleigh_quotient(Z, b), b, num_iteration

for i in [2**i for i in range(1, 15)]:
    Z = mu+sigma*numpy.random.randn(i, i)
    time.clock()
    ans = power_iteration(Z, 0.0001)
    end_time = time.clock()
    print('n = %d, num_iteration = %d, computation_time = %.5fs' % (i, ans[-1], end_time))

结果如下：

n = 2, num_iteration = 4, computation_time = 0.00016s
n = 4, num_iteration = 6, computation_time = 0.00045s
n = 8, num_iteration = 6, computation_time = 0.00337s
n = 16, num_iteration = 6, computation_time = 0.00374s
n = 32, num_iteration = 5, computation_time = 0.00420s
n = 64, num_iteration = 5, computation_time = 0.00470s
n = 128, num_iteration = 4, computation_time = 0.00616s
n = 256, num_iteration = 4, computation_time = 0.00917s
n = 512, num_iteration = 4, computation_time = 0.02576s
n = 1024, num_iteration = 4, computation_time = 0.09117s
n = 2048, num_iteration = 3, computation_time = 0.35692s
n = 4096, num_iteration = 3, computation_time = 1.36184s
n = 8192, num_iteration = 3, computation_time = 5.33402s
n = 16384, num_iteration = 3, computation_time = 22.03539s

注：time.clock()函数在Unix和Windows上的逻辑不同。本代码运行在Windows下。

Exercise 9.5: Singular values

Generate an $n \times 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?

for i in range(1, 21):
    p = numpy.random.random()
    C = numpy.random.rand(i, i)
    for x in numpy.nditer(C, op_flags=['readwrite']):
        x[...] = int(x>p)
    largest_singular_value = scipy.linalg.svd(C)[1].max()
    print(i, p, largest_singular_value/(i*(1-p)))

结果如下：

1 0.7062732095948546 1.0 3.404524315336278
2 0.3092046755501813 1.6180339887498951 1.1711384917367327
3 0.6980664733541699 1.4142135623730951 1.561285777130647
4 0.42589199638974373 2.6294927554336645 1.1450340088006958
5 0.6448744457574925 2.0 1.1263621984433392
6 0.34815120316043613 4.541265371562903 1.161124428850387
7 0.9184178628185704 1.4142135623730951 2.4764061826567354
8 0.28297257641663554 6.001927119649966 1.0463210545098764
9 0.1577365464789231 7.314497480172666 0.9649260440350683
10 0.8376846303366561 2.031442113282388 1.2515402068798378
11 0.0029423337276608397 11.0 1.0029510166033437
12 0.05532093632253121 10.927195150044199 0.9639248193902098
13 0.7840148431713382 4.026942985937468 1.4341950605330034
14 0.5128857973897747 6.879554650947621 1.0087916922742128
15 0.830317372713335 3.714399938573341 1.4593518884712393
16 0.6176825252651831 6.824246467296105 1.1156053081325834
17 0.2235203787788148 13.201967143664898 1.0001373910459754
18 0.0377629061862087 16.882941260965882 0.9747505964951693
19 0.24950047002444076 14.838343682548814 1.0405941986421847
20 0.1794485471751679 17.01989927285184 1.0371012819646532

猜测$n$，$p$和最大奇异值$S$的关系是$S \approx n(1-p)$.

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.

def closest_value(A, z):
    return A[numpy.abs(A-z).argmin()]