-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathiterative.py
201 lines (158 loc) · 6.05 KB
/
iterative.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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
import numpy as np
from lup_funcs import LUP_inplace, solve_from_LUP, solve_LUP
from hh_funcs import householder, householer_qr
def pow_it(A, x0, tol, maxit, store_iterations = False):
"""
For a matrix A, apply the power iteration algorithm with initial
guess x0, until either
||r|| < tol where
r = Ax - lambda*x,
or the number of iterations exceeds maxit.
:param A: an mxm numpy array
:param x0: the starting vector for the power iteration
:param tol: a positive float, the tolerance
:param maxit: integer, max number of iterations
:param store_iterations: if True, then return the entire sequence \
of power iterates, instead of just the final iteration. Default is \
False.
:return x: an m dimensional numpy array containing the final iterate, or \
if store_iterations, an mxmaxit dimensional numpy array containing all \
the iterates.
:return lambda0: the final eigenvalue.
"""
v = x0
if store_iterations:
V = np.zeros((v.size, maxit+1), dtype=A.dtype)
V[:,0] = x0
for i in range(1, maxit+1):
w = A.dot(v)
v = w / np.linalg.norm(w)
if store_iterations:
V[:,i] = v
lambda0 = v.conj().dot(A.dot(v))
if np.linalg.norm(A.dot(v) - lambda0 * v) < tol:
print(f"breaking at {i}!")
break
if store_iterations:
return V, lambda0
else:
return v, lambda0
def inverse_it(A, x0, mu, tol, maxit, store_iterations = False):
"""
For a Hermitian matrix A, apply the inverse iteration algorithm
with initial guess x0, using the same termination criteria as
for pow_it.
:param A: an mxm numpy array
:param mu: a floating point number, the shift parameter
:param x0: the starting vector for the power iteration
:param tol: a positive float, the tolerance
:param maxit: integer, max number of iterations
:param store_iterations: if True, then return the entire sequence \
of inverse iterates, instead of just the final iteration. Default is \
False.
:return x: an m dimensional numpy array containing the final iterate, or \
if store_iterations, an mxmaxit dimensional numpy array containing \
all the iterates.
:return l: a floating point number containing the final eigenvalue \
estimate, or if store_iterations, an m dimensional numpy array containing \
all the iterates.
"""
v = x0
if store_iterations:
V = np.zeros((v.size, maxit+1), dtype=A.dtype)
V[:,0] = v
# Not LU yet but will be when transformed in-place
LU = A - mu * np.eye(*A.shape)
p = LUP_inplace(LU)
for i in range(1, maxit+1):
w = solve_from_LUP(LU, v, p)
v = w / np.linalg.norm(w)
if store_iterations:
V[:,i] = v
lambda0 = v.conj().dot(A.dot(v))
if np.linalg.norm(A.dot(v) - lambda0 * v) < tol:
break
if store_iterations:
return V, lambda0
else:
return v, lambda0
def rq_it(A, x0, tol, maxit, store_iterations = False):
"""
For a Hermitian matrix A, apply the Rayleigh quotient algorithm
with initial guess x0, using the same termination criteria as
for pow_it.
:param A: an mxm numpy array
:param x0: the starting vector for the power iteration
:param tol: a positive float, the tolerance
:param maxit: integer, max number of iterations
:param store_iterations: if True, then return the entire sequence \
of inverse iterates, instead of just the final iteration. Default is \
False.
:return x: an m dimensional numpy array containing the final iterate, or \
if store_iterations, an mxmaxit dimensional numpy array containing \
all the iterates.
:return l: a floating point number containing the final eigenvalue \
estimate, or if store_iterations, an m dimensional numpy array containing \
all the iterates.
"""
v = x0
lambda0 = v.conj().dot(A.dot(v))
identity = np.eye(*A.shape, dtype=A.dtype)
if store_iterations:
V = np.zeros((v.size, maxit+1), dtype=A.dtype)
V[:,0] = v
for i in range(1, maxit+1):
w = solve_LUP(A - lambda0 * identity, v).reshape((v.size,))
v = w / np.linalg.norm(w)
if store_iterations:
V[:,i] = v
lambda0 = v.conj().dot(A.dot(v))
if np.linalg.norm(A.dot(v) - lambda0 * v) < tol:
break
if store_iterations:
return V, lambda0
else:
return v, lambda0
def pure_QR(A, maxit, tol=1e-8, in_place=True, return_iterations=False, termination=None):
"""
For matrix A, apply the QR algorithm and return the result.
:param A: an mxm numpy array
:param maxit: the maximum number of iterations
:param tol: termination tolerance
:param in_place: choice of method
:param return_iterations: boolean to choose whether to return the number of
QR iterations until convergence
:param termination: option to provide custom termination criteria
:return Ak: the result
"""
m = A.shape[0]
if return_iterations:
iterations = 0
# in_place = True avoids creating new NumPy arrays each loop; this however
# does not seem to noticably improve performance
if in_place:
AI = np.c_[A, np.eye(m)]
for i in range(maxit):
householder(AI, kmax=m)
AI[:,:m] = AI[:,:m].dot(AI[:,m:].conj().T)
AI[:,m:] = np.eye(m)
if np.linalg.norm(AI[:,:m][np.tril_indices(m, k=-1)])/m**2 < tol:
break
return AI[:,:m]
# Simplest method, which seems to perform just as well
else:
for i in range(maxit):
Q, R = householder_qr(A)
A = R.dot(Q)
if return_iterations:
iterations += 1
if termination:
check = termination(A)
else:
check = np.linalg.norm(A[np.tril_indices(m, k=-1)]) / m**2 < tol
if check:
break
if return_iterations:
return A, iterations
else:
return A