meine Aufgabe ist es eine Power iteration selbst zu implementieren und die Konvergenz zu überprüfen. Dazu habe ich eine Referenz Lösung (Eigenvektor/Eigenwert) mit der Funktion eigh berechnet. Bei jeder Iteration habe ich den relativen Fehler für den Eigenwert (Lama_max - Lama_k) und für den Eigenvektor (norm(eigenvect - z)) berechnet.
Zusätzlich habe ich noch die Geschwindigkeit der Konvergenz über den Quotienten des zweigrößten und größten Eigenwertes ((abs(second_val/lamda_max)**k)) berechnet.
Leider ist das Ergebnis nicht zufriedenstellend, da der Fehler nur ca. die ersten 3 Interationen kleiner wird und dann stagniert (siehe Graph nach dem Kompilieren). Ich finde den Fehler einfach nicht und es lässt mich wirklich verzweifeln...
Würde mich sehr freuen wenn jemand mal über den Code schauen könnte und eventuell den Fehler sieht.
LG,
Num133
Code: Alles auswählen
import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt
def power_iteration(A, z_0):
lamda_err = []
vec_err = []
speed_q = []
k = 0
err = 1e-10
z = z_0
# Compute reference eigenvalues/vectors
vals, vects = LA.eigh(A)
vals = list(vals)
maxcol = list(vals).index(max(vals))
lamda_max = max(vals)
eigenvect = vects[:, maxcol]
eigenvect = eigenvect / LA.norm(eigenvect)
# Speed
vals.remove(max(vals))
second_val = max(vals)
c = abs(second_val/lamda_max)
while True:
# Compute Rayleigh
lamda_k = ((z.T @ A @ z)/(z.T @ z))[0, 0]
# Compute convergence criteria
r = LA.norm((A @ z) - (lamda_k * z))
# Compute errors
lamda_err.append(abs((lamda_max - lamda_k)) / lamda_max)
direction = np.dot(eigenvect, z)
if direction < 0:
eigenvect = -eigenvect
vec_err.append(LA.norm(eigenvect - z) / (LA.norm(eigenvect)))
# Compute speed q
speed_q.append(c**k)
# Compute new vector
z = A @ z
z = z / LA.norm(z)
k = k + 1
if r < err:
break
print("Iterations")
print(k)
print("Eigenvalue error")
print(lamda_err)
print("Eigenvector error")
print(vec_err)
print("Speed of convergence q")
print(speed_q)
print("Approximated Eigenvalue")
print(lamda_k)
print("Approximated Eigenvector")
print(z)
print("Correct max Eigenvalue")
print(lamda_max)
print("Correct Eigenvector")
print(eigenvect)
return z, k, lamda_err, vec_err, speed_q
def plot_graph(k_it, lamda_err, vec_err, speed_q):
x = list(range(0, k_it, 1))
plt.plot(x, vec_err, 'b-', label='Eigenvector Error')
plt.plot(x, lamda_err, 'g-', label='Rayleigh Error')
plt.plot(x, speed_q, 'r-o', label='Convergence Speed')
plt.xlabel("number of iterations")
plt.ylabel("relative error")
plt.yscale('log')
plt.legend()
plt.show()
if __name__ == '__main__':
N = 50
A = np.random.randint(0, 2000, size=(N, N))
A_symm = (A + A.T)/2
z_start = np.random.randint(0, 2000, size=(N, 1))
z_0 = z_start / LA.norm(z_start)
z_vec, k_it, lamda_err, vec_err, speed_q = power_iteration(A, z_0)
plot_graph(k_it, lamda_err, vec_err, speed_q)
