-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathrandom_matrices.py
110 lines (101 loc) · 3.86 KB
/
random_matrices.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg as la
def random_mat(m):
"""
Return a random mxm matrix whose entries are drawn from the real normal distribution with mean 0, devation 1/sqrt(m)
"""
return np.random.randn(m,m) / m**.5
def random_upper_triangular(m):
"""
Random mxm upper triangular matrix, with normal entries
"""
M = random_mat(m)
return np.triu(M)
def plot_svals(svals, title=""):
svals = np.array(svals)
m = -np.log(np.min(svals))/np.log(2)
bins = [2**(i) for i in xrange(-int(m)-1, 1)]
# print m, bins, svals
plt.hist(svals, bins, normed = 1, cumulative=1)
plt.xscale('log')
plt.title(title)
plt.show()
#explore eigenvalues of random mxm matrices
def square_mats():
ms = []
avg_rads = []
avg_norms = []
svals_list = []
for exp in xrange(3, 7):
m = 2**exp
ms.append(m)
rs = []
two_norms = []
smallest_singular = []
for i in xrange(100):
A = random_mat(m)
two_norms.append(la.norm(A, ord=2))
smallest_singular.append(la.norm(A, ord=-2))
evals = la.eigvals(A)
rs.append(np.max(np.abs(evals)))
plt.scatter(np.real(evals), np.imag(evals))
plt.title("Eigenvalue distribution for random {}x{} matrices".format(m,m))
plt.xlabel("Real part")
plt.ylabel("Complex part")
plt.axes().set_aspect('equal', 'datalim')
plt.show()
svals_list.append(smallest_singular)
avg_rads.append(sum(rs)/len(rs))
avg_norms.append(sum(two_norms)/len(rs))
print "(a) It looks like the eigen-values fall in the complex unit ball! The spectral radius is around 1."
plt.plot(ms, avg_rads, label="Average spectral radii")
plt.plot(ms, avg_norms, label="Average 2-norms")
plt.legend()
plt.ylim(1,3)
plt.xlabel("Dimension")
plt.show()
print "(b) It appears that the 2-norm is bounded in the limit, but larger than the spectral radius."
for m, svals in zip(ms, svals_list):
plot_svals(svals, title="Cumulative distribution of sigma_min for m={}".format(m))
print "(c) It appears that as m scales, the smallest singular values become uniformly smaller."
def upper_triangular_mats():
ms = []
avg_rads = []
avg_norms = []
svals_list = []
for exp in xrange(3, 7):
m = 2**exp
ms.append(m)
rs = []
two_norms = []
smallest_singular = []
for i in xrange(100):
A = random_upper_triangular(m)
two_norms.append(la.norm(A, ord=2))
smallest_singular.append(la.norm(A, ord=-2))
evals = la.eigvals(A)
rs.append(np.max(np.abs(evals)))
plt.scatter(np.real(evals), np.imag(evals))
plt.title("Eigenvalues for random upper triangular {}x{} matrices".format(m,m))
plt.xlabel("Real part")
plt.ylabel("Complex part")
plt.axes().set_aspect('equal', 'datalim')
plt.show()
svals_list.append(smallest_singular)
avg_rads.append(sum(rs)/len(rs))
avg_norms.append(sum(two_norms)/len(rs))
print "(d)(a) The eigen-values are always real, and the spectral radius tends to zero."
plt.plot(ms, avg_rads, label="Average spectral radii")
plt.plot(ms, avg_norms, label="Average 2-norms")
plt.legend()
plt.ylim(0,3)
plt.xlabel("Dimension")
plt.show()
print "(d)(b) The spectral radius goes to zero, while the two-norm seems to be boundedly increasing."
for m, svals in zip(ms, svals_list):
plot_svals(svals, title="Cumulative distribution of sigma_min for m={}, upper triangular mats".format(m))
print "(d)(c) It appears that as m scales, the smallest singular values become uniformly smaller. Also, they go to zero very quickly."
if __name__ == "__main__":
square_mats()
upper_triangular_mats()