forked from PeterDenton/Eigenvector-Eigenvalue-Identity
-
Notifications
You must be signed in to change notification settings - Fork 0
/
E2.py
159 lines (118 loc) · 4.44 KB
/
E2.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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
import numpy as np
from functools import reduce
# Calculate the norm square of the jth element of the ith eigenvector using the alternate algorithm
# H is a square np.matrix
def vij_sq(H, i, j):
# Get the size of the matrix
n = H.shape[0]
# Get the eigenvalues of the matrix using numpy
eigenvalues = np.linalg.eigvalsh(H)
# Determine the jth minor
Mj = np.delete(np.delete(H, j, 0), j, 1)
# Get the eigenvalues of the minor
minor_eigenvalues = np.linalg.eigvalsh(Mj)
# Calculate the numerator
numerator = 1.
for k in range(n - 1):
numerator *= eigenvalues[i] - minor_eigenvalues[k]
# Calculate the denominator
denominator = 1.
for k in range(n):
if k == i:
continue
denominator *= eigenvalues[i] - eigenvalues[k]
return numerator / denominator
def eig_sq(H):
# Get the size of the matrix
n = H.shape[0]
res = np.zeros((n, n))
for i in range(n):
for j in range(n):
res[i, j] = vij_sq(H, j, i)
return res
# Calculate the norm square with performance improvements
def vij_sq_perf(H, i, j):
# Get the size of the matrix
n = H.shape[0]
# Get the eigenvalues of the matrix using numpy
eigenvalues = np.linalg.eigvalsh(H)
# Determine the jth minor
Mj = np.delete(np.delete(H, j, 0), j, 1)
# Get the eigenvalues of the minor
minor_eigenvalues = np.linalg.eigvalsh(Mj)
# Calculate the numerator
eigenvaluei = eigenvalues[i]
numerator = reduce(lambda x, y: x * y, [eigenvaluei - minor_eigenvalues[k] for k in range(n - 1)])
# Calculate the denominator
denominator = reduce(lambda x, y: x * y, [eigenvaluei - eigenvalues[k] for k in range(n) if k != i])
return numerator / denominator
def eig_sq_perf(H):
# Get the size of the matrix
n = H.shape[0]
res = np.zeros((n, n))
for i in range(n):
for j in range(n):
res[i, j] = vij_sq_perf(H, j, i)
return res
# Calculate the norm square of all the elements of all the eigenvectors using the alternate algorithm
# H is a square np.matrix
def eig_sq_vect(H):
# Get size of the matrix
n = H.shape[0]
# Get eigenvalues of the matrix
eigvals = np.linalg.eigvalsh(H)
# Create a matrix of all the minors minors and a matrix of eigenvalues as required in the algorithm
minors = np.empty((n, n - 1, n - 1), dtype=H.dtype)
eig_matrix = np.empty((n, n - 1))
for i in range(n):
minors[i] = np.delete(np.delete(H, i, 0), i, 1)
eig_matrix[i] = np.delete(eigvals, i)
# Get eigenvalues of all the minors
minor_eigvals = np.linalg.eigvalsh(minors)
# Calculate the numerators and the denominators
numerator = np.empty((n, n))
denominator = np.empty((n,))
for i in range(n):
numerator[i] = np.prod(eigvals[i] - minor_eigvals, axis=1)
denominator[i] = np.prod(eigvals[i] - eig_matrix[i])
return np.divide(numerator.T, denominator.reshape((1, n)))
# Calculate the norm square of all the elements of all the eigenvectors using numpy's built in eigenvector calculator
# H is a square np.matrix
def eig_sq_np(H):
# Get the eigenvalues and eigenvalues of the matrix
eigenvalues, eigenvectors = np.linalg.eigh(H)
return abs(eigenvectors) ** 2
# Prints true for each element that agrees between the E2 method and the method built into numpy
# H is a square np.matrix
# mode determines comparison with normal or improved function
def compare(H, mode=0):
# The tolerance for comparison purposes
tol = 1e-12
# Select function depending upon the value of mode
if mode == 0:
res = eig_sq(H)
elif mode == 1:
res = eig_sq_perf(H)
elif mode == 2:
res = eig_sq_vect(H)
else:
raise IndexError
# Return the comparison matrix
return abs(res - eig_sq_np(H)) < tol
if __name__ == "__main__":
# Generate a random Hermitian matrix
# Dimension of the matrix
n = 8
H = np.random.rand(n, n) + 1j * np.random.rand(n, n)
H = 0.5 * (H + H.T)
# Print the second element of the first eigenvector as calculated in both methods
print("\n--- using the algorithm ---")
basic = eig_sq(H)
print("basic implementation: ", basic[0, 1])
perf = eig_sq_perf(H)
print("performant implementation: ", perf[0, 1])
vect = eig_sq_vect(H)
print("vectorized implementation implementation: ", vect[0, 1])
print("\n--- using numpy ---")
nump = eig_sq_np(H)
print(nump[0, 1])