for n = 200, m = 500.
生成n×m矩阵A,其中各元素是服从高斯分布的随机数。生成常对角矩阵B。
代码:
import numpy as np
import scipy.linalg as la
n = 200
m = 500
A = np.mat(np.random.randn(n, m))
col_1 = np.random.random(size=(m, 1))
row_1 = np.random.random(size=(1, m))
B = np.mat(la.toeplitz(col_1, row_1))
Exercise 9.1: Matrix operations
Calculate A + A, AA > ,A > A and AB. Write a function that computes A(B − λI) for any λ.
代码:
def cal(lambdaa):
return A*(B-lambdaa*I)
print(
'A = '+str(A),
'B = '+str(B),
'A+A = '+str(A+A),
'A*A.T = '+str(A*A.T),
'A.T*A = '+str(A.T*A),
'A*B = '+str(A*B),
'I = '+str(I),
'A(B − λI) =' + str(cal(5)),
sep='\n\n'
)
输出:
A = [[ 0.04195854 -0.38710441 0.02414068 ... 0.12998791 -1.09805947
-0.49228629]
[ 0.28992152 -0.25311269 1.51203024 ... 0.87529587 0.42869605
-1.01772641]
[-0.01180153 -0.52115087 -1.06669412 ... 0.09279688 -0.18493147
0.47616983]
...
[-1.59083133 2.15117569 0.91646871 ... -0.07697636 0.58731201
0.35030513]
[-0.47906178 0.50324141 -0.14811297 ... -2.20791176 -1.22901378
-0.64334274]
[ 0.91530925 -1.58723489 -0.03392294 ... -0.91057087 -1.15192596
0.89204535]]
B = [[0.82404747 0.34030171 0.04510667 ... 0.97886067 0.96370754 0.22601283]
[0.77107993 0.82404747 0.34030171 ... 0.52657226 0.97886067 0.96370754]
[0.01290139 0.77107993 0.82404747 ... 0.38751267 0.52657226 0.97886067]
...
[0.89640885 0.48097235 0.85224666 ... 0.82404747 0.34030171 0.04510667]
[0.99096391 0.89640885 0.48097235 ... 0.77107993 0.82404747 0.34030171]
[0.42959832 0.99096391 0.89640885 ... 0.01290139 0.77107993 0.82404747]]
A+A = [[ 0.08391708 -0.77420882 0.04828136 ... 0.25997582 -2.19611895
-0.98457258]
[ 0.57984304 -0.50622538 3.02406047 ... 1.75059174 0.85739211
-2.03545282]
[-0.02360306 -1.04230174 -2.13338824 ... 0.18559375 -0.36986295
0.95233966]
...
[-3.18166265 4.30235138 1.83293741 ... -0.15395271 1.17462402
0.70061025]
[-0.95812355 1.00648281 -0.29622593 ... -4.41582353 -2.45802756
-1.28668548]
[ 1.8306185 -3.17446977 -0.06784589 ... -1.82114174 -2.30385192
1.7840907 ]]
A*A.T = [[ 5.27537748e+02 1.34696479e-01 -1.24809925e+01 ... -3.93112193e+01
2.15343154e+01 4.26919125e+01]
[ 1.34696479e-01 5.41951934e+02 1.90018025e+01 ... -1.85704375e+01
3.20290916e+01 -1.35193393e+01]
[-1.24809925e+01 1.90018025e+01 4.60985958e+02 ... 3.89682737e+01
-1.37186793e+01 3.72866375e+01]
...
[-3.93112193e+01 -1.85704375e+01 3.89682737e+01 ... 4.88881188e+02
-1.58172861e+01 3.50264454e+00]
[ 2.15343154e+01 3.20290916e+01 -1.37186793e+01 ... -1.58172861e+01
4.71088145e+02 -1.91876202e+01]
[ 4.26919125e+01 -1.35193393e+01 3.72866375e+01 ... 3.50264454e+00
-1.91876202e+01 5.07528739e+02]]
A.T*A = [[181.59339518 -3.4769892 13.47451816 ... -2.85963471 4.77597727
29.82138842]
[ -3.4769892 205.4581638 -3.30189924 ... 0.36806718 20.29013161
-20.0666659 ]
[ 13.47451816 -3.30189924 181.18403999 ... 11.05256883 -16.15320599
4.64771307]
...
[ -2.85963471 0.36806718 11.05256883 ... 204.89356816 -1.18547402
17.16772537]
[ 4.77597727 20.29013161 -16.15320599 ... -1.18547402 198.06877941
-13.69891785]
[ 29.82138842 -20.0666659 4.64771307 ... 17.16772537 -13.69891785
181.26927499]]
A*B = [[ 4.25658264 6.68511215 14.31291597 ... 15.29729322 22.04987366
7.30773622]
[ -4.38481245 8.37493151 12.30946257 ... 12.13157489 7.13387606
6.31249827]
[ -6.37549539 4.61405134 2.29595035 ... 7.74256197 -8.35348835
-9.42119234]
...
[ 8.90184728 9.78111981 2.04445731 ... 8.66266684 -3.54664847
5.99144587]
[ 0.78702146 -6.53679439 -12.00846036 ... 4.2966042 -5.94174445
9.52998304]
[-20.97034986 -11.18665034 -3.79174184 ... -21.56821008 -5.90990218
-14.0203237 ]]
I = [[1. 0. 0. ... 0. 0. 0.]
[0. 1. 0. ... 0. 0. 0.]
[0. 0. 1. ... 0. 0. 0.]
...
[0. 0. 0. ... 1. 0. 0.]
[0. 0. 0. ... 0. 1. 0.]
[0. 0. 0. ... 0. 0. 1.]]
A(B − λI) =[[ 4.04678994 8.62063419 14.19221256 ... 14.64735367 27.54017104
9.76916766]
[ -5.83442004 9.64049495 4.74931139 ... 7.75509554 4.99039579
11.40113032]
[ -6.31648775 7.21980568 7.62942095 ... 7.27857758 -7.42883099
-11.8020415 ]
...
[ 16.85600391 -0.97475862 -2.53788623 ... 9.04754861 -6.48320852
4.23992024]
[ 3.18233034 -9.05300142 -11.26789553 ... 15.33616302 0.20332445
12.74669674]
[-25.54689611 -3.2504759 -3.62212712 ... -17.01535573 -0.15027239
-18.48055046]]
Exercise 9.2: Solving a linear system
Generate a vector b with m entries and solve Bx = b.
代码:
b = np.random.random(size=(m, 1))
x = np.mat(np.linalg.solve(B, b))
print('x = '+str(x))
输出:
x = [[ 4.38225617e-01]
[-2.62604279e-02]
[ 5.77628412e-02]
...
[-3.41521726e-01]
[ 9.07959804e-02]
[-5.83334792e-01]]
经验证Bx = b
Exercise 9.3: Norms
Compute the Frobenius norm of A: kAk F and the infinity norm of B: kBk ∞ . Also find the largest and
smallest singular values of B.
代码:
AF = np.linalg.norm(A)
Binf = np.linalg.norm(B, np.inf)
s = la.svdvals(B)
print(
'Frobenius norm of A: '+str(AF),
'Infinity norm of B: '+str(Binf),
'The largest singular values of B: '+str(s.max()),
'The smallest singular values of B: ' + str(s.min()),
sep = '\n'
)
输出:
Frobenius norm of A: 315.95558444399796
Infinity norm of B: 259.47901617509274
The largest singular values of B: 253.9171417278695
The smallest singular values of B: 0.04916465790673367
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 time
def power_iteration(A):
n, d = A.shape
count = 0
v = np.ones(d)
ev = 0
begin_time = time.clock()
while True:
count += 1
Av = A.dot(v)
ev_new = np.linalg.norm(Av)
v_new = Av / ev_new
if np.abs(ev - ev_new) < 0.01:
break
v = v_new
ev = ev_new
end_time = time.clock()
print('Time: ' + str(end_time-begin_time))
return ev_new, v_new, count
n = 200
Z = np.random.randn(n, n)
ev, v, count = power_iteration(Z)
print(
'eigenvalue: ' + str(ev),
'eigenvector: ' + str(v),
'number of iteration: ' + str(count),
sep='\n'
)
输出:
Time: 0.0018198326227696195
eigenvalue: 15.81618091989844
eigenvector: [-0.03700825 0.0029305 0.01227961 0.02983836 0.01829982 0.02772856
0.02107871 0.09893484 0.0350257 0.04419609 0.04727178 0.056736
-0.04377146 0.09618859 -0.01862322 0.01417478 -0.02309586 -0.00194714
-0.03120478 -0.02046295 -0.0968875 -0.03675232 0.03318669 0.10242163
0.05797829 0.07278333 0.02818215 -0.01875538 0.12128094 -0.07273369
-0.01526538 -0.17105178 -0.11733502 -0.01788494 -0.0601756 -0.06721134
0.09800221 0.01506051 0.02491852 -0.12336477 0.02428129 -0.02547432
-0.03982852 0.09438698 -0.09507363 -0.10752629 -0.03039604 0.02625876
0.02946467 -0.03132434 -0.00859429 0.02732496 -0.02247782 0.06659113
-0.01402789 0.04534168 0.03565373 0.00755469 -0.0681837 -0.06179226
-0.01214626 -0.03232888 -0.00067058 -0.06898048 0.02596048 -0.04140793
-0.13078973 0.05529428 -0.05294517 0.00409866 0.01130606 0.08446122
0.16050908 0.01138752 0.03961415 -0.1426926 0.05989797 -0.13517472
0.00196125 -0.06839017 0.01876062 0.00799628 0.07576671 0.13489427
-0.08642153 0.01101552 -0.01891506 -0.03754352 -0.11884894 0.05107088
0.00951774 -0.02215805 -0.04574779 -0.08041138 0.08242144 -0.19971271
-0.03669991 -0.0268403 0.07627959 -0.09874723 -0.03435932 -0.00241873
-0.06670095 -0.01555391 0.05338037 -0.21958027 -0.05404186 0.03915557
-0.10399761 0.08424245 0.07007432 0.02019906 -0.00034468 -0.01125532
-0.04812383 -0.07710452 0.09902331 0.00081299 -0.08204022 -0.05218078
0.01894426 -0.0911309 -0.03287158 0.02815841 0.02207783 -0.05711911
0.00441251 -0.21093863 0.0020175 -0.02082109 -0.11669596 -0.06667732
0.07087201 0.03347427 0.02888406 -0.02505937 -0.09300211 0.07533929
0.05001491 -0.02080881 0.09205245 -0.07718376 -0.01268754 0.0559997
0.08342342 0.08266444 -0.05132904 -0.04049663 -0.11326853 -0.14141689
0.1126093 -0.01487698 0.09328524 -0.01179341 0.09013115 0.01899532
0.04128708 0.04779245 0.01369805 -0.05740939 0.04079265 -0.00880687
0.02980492 -0.12020824 0.02103489 -0.04582447 0.05758544 -0.00679323
0.04331108 0.03973211 -0.01544934 0.09382855 -0.06072975 -0.0369994
0.08413063 0.09077527 0.08635241 0.04281156 0.04731971 -0.04836914
0.00405065 0.02900649 0.08715658 0.09654352 0.05667637 0.00851858
0.02996912 -0.12203476 -0.06564825 0.0525448 0.03626988 0.10106428
0.09465614 -0.08401987 -0.20462801 -0.09707087 -0.07242656 -0.0545679
-0.14721702 0.02890988]
number of iteration: 28
所需时间如下图,矩阵大小从100×100到2000×2000,步长为100,对每个大小运行10次取均值。
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 matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
def svd(n, p):
C = np.random.binomial(1, p, n)
for _ in range(n-1):
row = np.random.binomial(1, p, n)
C = np.row_stack((C, row))
return la.svdvals(C)
X = np.arange(100, 1100, 100)
Y = np.arange(0.1, 1.1, 0.1)
X, Y = np.meshgrid(X, Y)
Z = [svd(x,y).max() for x, y in zip(X.flat, Y.flat)]
Z = np.reshape(Z, (10,10))
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.viridis)
plt.xlabel('n')
plt.ylabel('p')
plt.show()
输出:
可知最大奇异值虽n,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.
Hint: Use the built-in functionality of Numpy rather than writing code to find this value manually. In
particular, use brackets and argmin.
代码:
z = np.random.random()
a = np.reshape(A, (1, 100000))
index = np.argmin(np.abs((a-z)))
print(
'z = ' + str(z),
'the element in A that is closest to z: ' + str(a[0,index]),
sep = '\n'
)
输出:
z = 0.7783356000465671
the element in A that is closest to z: 0.7782910522477339