Generate matrices A, with random Gaussian entries, B, a Toeplitz matrix, where A ∈ R
n×m and B ∈ Rm×m,for n = 200, m = 500.
import numpy as np
from scipy.linalg import toeplitz
n_size = 200
m_size = 500
mu = 0.0
sigma = 0.1
np.random.seed(0)
A = np.random.normal(mu, sigma, (n_size, m_size))
print(A)
print(A.shape)
base = list(range(1, m_size + 1))
# print(base)
B = toeplitz(base, base)
print(B)
print(B.shape)
输出:
[[ 0.17640523 0.04001572 0.0978738 ..., -0.0514234 -0.10180419
-0.00778548]
[ 0.03827324 -0.00342423 0.10963468 ..., 0.00941923 -0.11476109
-0.03581141]
[ 0.05559627 0.08924739 -0.04223148 ..., -0.02467936 -0.10793432
-0.01142256]
...,
[-0.16058274 -0.035732 -0.09369679 ..., -0.1123346 -0.06918121
0.05208148]
[-0.00682909 -0.09105953 0.01875499 ..., -0.0513064 0.00768273
-0.09786463]
[ 0.02218512 -0.09420463 -0.0277842 ..., 0.0726182 -0.11678305
-0.12852076]]
(200, 500)
[[ 1 2 3 ..., 498 499 500]
[ 2 1 2 ..., 497 498 499]
[ 3 2 1 ..., 496 497 498]
...,
[498 497 496 ..., 1 2 3]
[499 498 497 ..., 2 1 2]
[500 499 498 ..., 3 2 1]]
(500, 500)
Exercise 9.1: Matrix operations
Calculate A + A, AA^T, A^TA and AB. Write a function that computes A(B − λI) for any λ.
# 计算A + A
sum = A + A
print('sum:\n', sum)
# 计算A * (A^T)
res1 = np.dot(A, A.T)
print('res1:\n', res1)
# 计算(A^T) * A
res2 = np.dot(A.T, A)
# 计算A * B
res3 = np.dot(A, B)
print('res3:\n', res3)
def func(A, B, s, size): # 计算A(B − λI)
return np.dot(A, (B - s * np.identity(size)))
result = func(A, B, 0.5, m_size)
print('result:\n', result)
输出:
sum:
[[ 0.35281047 0.08003144 0.1957476 ..., -0.10284679 -0.20360838
-0.01557095]
[ 0.07654649 -0.00684846 0.21926937 ..., 0.01883846 -0.22952219
-0.07162282]
[ 0.11119254 0.17849478 -0.08446296 ..., -0.04935871 -0.21586863
-0.02284511]
...,
[-0.32116548 -0.07146399 -0.18739357 ..., -0.22466921 -0.13836242
0.10416295]
[-0.01365817 -0.18211906 0.03750997 ..., -0.1026128 0.01536546
-0.19572927]
[ 0.04437024 -0.18840927 -0.05556839 ..., 0.1452364 -0.2335661
-0.25704153]]
res1:
[[ 4.984799 -0.19033257 -0.22672115 ..., 0.10982337 0.33523647
0.2211457 ]
[-0.19033257 4.77802726 -0.05226853 ..., 0.1097388 0.19826363
-0.20766612]
[-0.22672115 -0.05226853 4.66202048 ..., -0.24424636 0.10341813
0.07524698]
...,
[ 0.10982337 0.1097388 -0.24424636 ..., 4.78065551 -0.21885776
-0.1114218 ]
[ 0.33523647 0.19826363 0.10341813 ..., -0.21885776 5.18029085
-0.29804843]
[ 0.2211457 -0.20766612 0.07524698 ..., -0.1114218 -0.29804843
4.28037337]]
res3:
[[ -770.16766139 -768.54712896 -766.84656508 ..., 137.33964977
136.29110713 135.03895612]
[ -21.2917197 -17.95722443 -14.62957762 ..., -1604.79749053
-1607.75429431 -1610.94062028]
[ 107.89619671 106.54753873 105.37737552 ..., 620.30765296
622.00621723 623.48891285]
...,
[ 270.00991989 267.0670624 264.05274091 ..., 1038.2843583
1040.94024978 1043.45777884]
[ 1299.92652164 1295.17208978 1290.23553886 ..., 1065.34345608
1070.26459357 1075.20109653]
[ -558.78922472 -557.96853555 -557.33625563 ..., 170.65842728
170.37271597 169.85343857]]
result:
[[ -770.25586401 -768.56713682 -766.89550198 ..., 137.36536147
136.34200923 135.04284885]
[ -21.31085632 -17.95551232 -14.68439496 ..., -1604.80220015
-1607.69691377 -1610.92271458]
[ 107.86839858 106.50291503 105.39849126 ..., 620.31999264
622.06018438 623.49462413]
...,
[ 270.09021126 267.08492839 264.0995893 ..., 1038.3405256
1040.97484038 1043.4317381 ]
[ 1299.92993618 1295.21761955 1290.22616136 ..., 1065.36910928
1070.26075221 1075.25002885]
[ -558.80031728 -557.92143323 -557.32236354 ..., 170.62211818
170.4311075 169.91769895]]
Exercise 9.2: Solving a linear system
Generate a vector b with m entries and solve Bx = b.
vector_b = np.random.random((500, 1))
# print('vector_b:\n', vector_b)
sloution = np.linalg.solve(B, vector_b)
print('sloution:\n', sloution)
输出:
sloution:
[[-0.19933145]
[ 0.46112091]
[-0.28651194]
[-0.36007313]
[ 0.45699715]
[-0.1136803 ]
[ 0.06858123]
[-0.02439619]
[ 0.01277973]
[ 0.07296403]
[ 0.10677244]
[-0.53755885]
[ 0.70831459]
[-0.50730103]
[-0.12370938]
...
Exercise 9.3: Norms
Compute the Frobenius norm of A: ||A||F and the infinity norm of B: ||B||∞. Also find the largest and smallest singular values of B.
Frobenius_norm_A = np.linalg.norm(A, ord='fro')
print('Frobenius_norm_A:\n', Frobenius_norm_A)
infinity_norm_B = np.linalg.norm(B, np.inf)
print("infinity_norm_B:\n", infinity_norm_B)
U, sigma_singular, VT = np.linalg.svd(B)
max_singular_values = max(sigma_singular)
min_singular_values = min(sigma_singular)
print('max_singular_values:', max_singular_values)
print('min_singular_values:', min_singular_values)
输出:
Frobenius_norm_A:
31.5386895558
infinity_norm_B:
125250.0
max_singular_values: 87334.5204564
min_singular_values: 0.500004934835
Exercise 9.4: Power iteration
Generate a matrix Z, n × 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 numpy as np
from scipy.linalg import toeplitz
import time
def power_iteration(A, n_size):
u = np.random.randint(0, 10, size=n_size)
num_iter = 0
start = time.clock()
next_lam = 1
pre_lam = 0
while abs(next_lam - pre_lam) > 0.00001:
pre_lam = next_lam
num_iter += 1
v = np.dot(A, u)
next_lam = max(abs(v))
u = v/next_lam
end = time.clock()
return u, next_lam, num_iter, end - start
n_size = 200
m_size = 500
mu = 0.0
sigma = 0.1
np.random.seed(0)
Z = np.random.normal(mu, sigma, (n_size, n_size))
eigenvector, lam, num_iters, time = power_iteration(Z, n_size)
print('the largest eigenvalue: ', lam)
print('the corresponding eigenvector: ', eigenvector)
print('the number of iterations: ', num_iters)
print('the time of iterations: ', time)
这里取n = 200输出:
the largest eigenvalue: 1.38374435975
the corresponding eigenvector: [-0.23448063 -0.23857682 -0.45558981 -0.22375566 0.49945187 -0.04983987
0.22586259 0.27561021 -0.25905579 0.27833055 0.10392899 0.10268637
0.27160415 0.21686688 -0.41957299 0.17116327 -0.17859623 0.16854047
-0.74401348 -0.41213827 0.42046522 0.34716834 0.05462962 0.68213181
-0.26967417 0.22156731 -0.05764227 0.15715657 -0.34579995 0.22889528
-0.43437992 -0.01709412 0.27750828 -0.32423244 -0.29912263 -0.33907246
-0.13649076 -0.33804399 -0.12328962 0.07250919 -0.0715865 0.2687078
0.12319253 0.07354623 -0.09890419 0.1051427 -0.37148188 -0.67210815
0.21043618 0.04296881 0.35635652 0.28666593 0.05953488 -0.52401179
0.10210595 0.58416896 0.49189862 0.01899168 0.05358757 0.28626702
-0.31010217 -0.75048177 0.12335981 -0.80049872 0.65508207 0.00988996
0.92204464 -0.41418936 -0.16885457 -0.16748709 -0.11395704 0.16557866
-0.22648951 -0.09870179 0.57416871 -0.07121615 -0.62246095 -0.20429331
-0.49580195 0.30371582 -0.75810448 -0.00355056 0.41848717 0.07971469
0.28304665 0.48246567 -0.28948918 0.57802151 0.17559049 -0.14606612
-0.00850067 0.43262172 0.04759938 0.52080095 0.06836337 -0.14321346
0.20378579 -0.05444155 -0.49864496 -0.11220402 0.03517885 0.53004052
0.15553566 0.25167114 0.00179509 -0.62592081 -0.31307944 -0.12592512
0.05764811 0.10770587 -0.26141666 0.32830107 -0.23000649 0.32579796
-0.42460092 -0.54333748 -0.07546899 -0.1202275 -0.22148577 0.46604526
-0.37403106 -0.19443784 -0.02189891 0.3196339 0.45715056 -0.95319163
-0.11534788 0.3458295 -0.07180602 -0.18107832 -0.4378873 -0.70818384
0.06307144 -0.75049885 0.15457617 -0.65519509 -0.24892478 0.04379539
-0.05674392 -0.45369933 -0.60605137 0.14269135 -0.15789244 0.71103203
0.10365828 -0.50814752 -0.3803296 -0.3921604 0.34154827 0.03067514
0.36761996 -0.01638594 -0.47892651 -0.62457798 0.62175791 0.60974948
-0.34869693 -0.15907353 0.47225017 0.357265 -0.01573852 -0.36644187
1. 0.11222906 0.14525716 -0.24414313 0.84366931 -0.77558439
0.61776754 0.08202249 -0.27274454 0.05338878 -0.58895049 0.43344831
0.32846749 -0.06703431 0.06770655 0.06771246 0.28469391 -0.04991912
-0.03054763 -0.07803063 -0.23599849 0.26720273 0.05086394 -0.05258738
-0.48809871 0.40606355 -0.34869432 0.07958705 0.04265011 -0.32885641
-0.35524689 0.68314129 -0.14964527 -0.222974 0.09612197 -0.42326915
0.11732494 0.15071393]
the number of iterations: 1813
the time of iterations: 0.13679134993804806
n = 100时:
the largest eigenvalue: 1.00512510268
the corresponding eigenvector: [ 0.25042316 0.45906092 -0.43155535 0.17084169 -0.31571496 -0.32947339
0.09475778 -0.07051941 -0.47606022 -0.44952443 1. -0.2734709
0.83039778 -0.07502614 -0.30195033 0.43850008 -0.28705981 -0.34490665
-0.8302928 -0.14746054 -0.20792767 0.46113797 0.11940885 -0.43873384
-0.4026181 0.18259404 -0.29167808 -0.57137324 -0.26113824 0.0447006
0.6175856 -0.20932656 0.1497286 0.40798432 -0.48649966 0.05361968
-0.15613001 -0.63377699 -0.19008996 -0.50509828 0.20765657 -0.67042046
-0.04948563 0.1925169 0.30199847 -0.05921492 0.08473466 0.40901562
0.56484163 0.60907853 0.52669863 0.50567972 -0.04177342 0.231218
-0.87164059 -0.18402144 0.25972169 -0.06419379 0.64224447 -0.09154909
-0.26810684 0.01138896 -0.40211867 -0.11467527 0.52627956 0.08511445
0.06288489 -0.06932288 0.19850499 -0.0042502 -0.44447338 0.00977573
-0.26873374 0.02494805 0.4358977 -0.4499939 0.42533057 -0.00221055
-0.3709814 -0.44295519 0.0393151 0.24926072 0.26211296 0.29795194
0.06267313 0.02638521 -0.34356405 0.40541127 -0.04932551 -0.53103175
-0.07837982 0.02486878 -0.0533286 -0.31930375 -0.20309002 -0.00346277
0.37197201 0.34570481 0.2211697 0.48039722]
the number of iterations: 54074
the time of iterations: 1.5376810345121839
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
def one_or_zero(p):
if np.random.random() > 1 - p:
return 1
else:
return 0
def largest_singular_value(A):
U, S, V = np.linalg.svd(A)
return max(S)
for n in range(2, 200, 10):
C = np.zeros((n, n))
for k in range(1, 10):
p = k / 10
for i in range(n):
for j in range(n):
C[i][j] = one_or_zero(0.5)
value = largest_singular_value(C)
print('the largest singular value: ', value, 'p: ', p, 'n: ', n)
输出:
the largest singular value: 1.41421356237 p: 0.1 n: 2
the largest singular value: 1.0 p: 0.2 n: 2
the largest singular value: 1.41421356237 p: 0.3 n: 2
the largest singular value: 1.61803398875 p: 0.4 n: 2
the largest singular value: 1.0 p: 0.5 n: 2
the largest singular value: 1.0 p: 0.6 n: 2
the largest singular value: 1.0 p: 0.7 n: 2
the largest singular value: 1.0 p: 0.8 n: 2
the largest singular value: 1.61803398875 p: 0.9 n: 2
the largest singular value: 6.09803062927 p: 0.1 n: 12
the largest singular value: 5.36168216768 p: 0.2 n: 12
the largest singular value: 6.63331122699 p: 0.3 n: 12
the largest singular value: 6.13206436997 p: 0.4 n: 12
the largest singular value: 6.71450698537 p: 0.5 n: 12
the largest singular value: 6.43079837273 p: 0.6 n: 12
the largest singular value: 6.57675764838 p: 0.7 n: 12
the largest singular value: 6.72227602896 p: 0.8 n: 12
the largest singular value: 6.43344014997 p: 0.9 n: 12
the largest singular value: 11.5611945294 p: 0.1 n: 22
the largest singular value: 10.861829775 p: 0.2 n: 22
the largest singular value: 11.2604419494 p: 0.3 n: 22
the largest singular value: 10.7399483624 p: 0.4 n: 22
the largest singular value: 11.427691543 p: 0.5 n: 22
the largest singular value: 11.3774277544 p: 0.6 n: 22
the largest singular value: 11.1564436536 p: 0.7 n: 22
the largest singular value: 12.663423727 p: 0.8 n: 22
the largest singular value: 10.861248278 p: 0.9 n: 22
...
the largest singular value: 87.0336062627 p: 0.1 n: 172
the largest singular value: 86.8825137409 p: 0.2 n: 172
the largest singular value: 86.1686471393 p: 0.3 n: 172
the largest singular value: 86.7796838102 p: 0.4 n: 172
the largest singular value: 87.3335809673 p: 0.5 n: 172
the largest singular value: 86.6432805456 p: 0.6 n: 172
the largest singular value: 86.3703671617 p: 0.7 n: 172
the largest singular value: 87.0565486864 p: 0.8 n: 172
the largest singular value: 86.0200230692 p: 0.9 n: 172
the largest singular value: 90.7639134334 p: 0.1 n: 182
the largest singular value: 91.5999226806 p: 0.2 n: 182
the largest singular value: 91.913131263 p: 0.3 n: 182
the largest singular value: 92.0368133674 p: 0.4 n: 182
the largest singular value: 91.29022559 p: 0.5 n: 182
the largest singular value: 91.089694189 p: 0.6 n: 182
the largest singular value: 91.2739383084 p: 0.7 n: 182
the largest singular value: 91.5252808991 p: 0.8 n: 182
the largest singular value: 91.9294664348 p: 0.9 n: 182
the largest singular value: 95.9877956052 p: 0.1 n: 192
the largest singular value: 96.3236948637 p: 0.2 n: 192
the largest singular value: 96.1700156123 p: 0.3 n: 192
the largest singular value: 97.0419055764 p: 0.4 n: 192
the largest singular value: 96.3372471599 p: 0.5 n: 192
the largest singular value: 96.6318842389 p: 0.6 n: 192
the largest singular value: 96.2941907085 p: 0.7 n: 192
the largest singular value: 96.8051842424 p: 0.8 n: 192
the largest singular value: 96.2303337938 p: 0.9 n: 192
可以看到the largest singular value在n比较小时其值受到p的影响还是比较大,但在当n比较大时,其值并没有太大的变化,而且n比较时的数值并没有太大的代表性,所以可以只看n比较大时the largest singular value的变化,而对于p相同时,n增加10,the largest singular value的变化还是很明显的,并且随着n的增大而增大。
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.
Hint: Use the built-in functionality of Numpy rather than writing code to find this value manually. In particular, use brackets and argmin.
import numpy as np
def get_nearest_neighbor(z, A):
B = abs(A - z)
return A[np.argmin(B)]
z = np.random.rand() * 10
A = np.random.rand(100) * 10
print('z: ', z)
print('A: ', A)
print('the nearest neighbor: ', get_nearest_neighbor(z, A))
输出:
z: 5.140756589959544
A: [ 2.38548496 9.49087376 1.28884664 0.63702924 6.67374531 9.51206904
5.84934177 8.86422373 4.18638771 4.99405837 3.96835962 7.17309549
9.3484444 1.04547294 1.88207947 4.10086577 5.2077526 3.7403351
0.49729221 2.69092716 4.49822518 2.48918472 9.79170684 1.99704282
6.90493967 1.35927367 0.58062983 3.58336034 3.34607267 1.63150145
7.13385139 5.0682374 2.87122973 7.39368304 6.20506582 7.67409263
2.50880758 2.15953104 3.44241519 5.26641213 7.80353381 9.44048273
8.34729033 1.8021764 6.61065728 6.02415552 2.81593119 2.72401768
5.28152751 3.18988028 3.57815635 8.04498665 6.87776589 6.5423706
4.53766417 4.68319582 7.91227892 2.09053791 6.13504116 5.27837593
3.90897699 6.74414173 4.44737472 4.3569735 9.65642109 6.08779349
5.45481335 3.04881375 2.23628713 0.22330483 4.68740936 1.07713896
0.3546105 1.53122082 7.82766266 8.33500204 8.45982052 5.94027026
5.80184327 8.2297121 3.49990858 8.6759586 5.70309204 7.70432551
8.30971419 9.86307761 6.39693714 1.44703552 9.81276596 4.73014594
0.61624376 3.24323494 7.41945733 8.60669928 9.49122329 6.63813434
5.43821111 3.24123646 4.52327842 0.33085725]
the nearest neighbor: 5.20775260306