(Twitter)
0
Zur Lösung des linearen Ausgleichsproblems ∥Ax − ⃗b∥2 mit ⃗b∈R^m, A∈R^m×n, m≥n und rang (A)=n, wollen wir nun die QR-Zerlegung verwenden. Dazu implementieren wir den Schneller Least-Squares mit Householder-Algorithmus. Gehen Sie wie folgt vor:
(a) Implementieren Sie dazu folgende Methoden:
- Householder_Vektor(x) liefert einen Householder-Vektor zum Vektor x zurück.
- Householder_Mult(A,v) multipliziert die Matrix A mit der Householder-Matrix zum Householder Vektor v und liefert das Produkt zurück
- Householder_QR(A) berechnet eine QR-Zerlegung der Householder-Matrix A bei der der essentielle Teil der Householder-Matrizen unter der Diagonalen von A gespeichert wird und liefert die resultierende Matrix zurück.
- Householder_Least_Squares(A,b) berechnet die kleinste Quadrate Lösung x von Ax=b und liefert x sowie das Residuum ρLS zurück.
Beachte, dass in diesen Algorithmen die Dimensionen der Vektoren und Matrizen variieren. Tipp: Dafür dürfen Kopien von Teilvektoren oder Teilmatrizen geeigneter Dimensionen angelegt werden.
(b) Testen Sie Ihr Verfahren, in dem Sie mithilfe der QR-Zerlegung ein Ausgleichsproblem mit
den folgenden Daten lösen und vergleichen Sie Ihr Ergebnis mit der Lösung, die Sie über
die Normalengleichung berechnen (mit ihrer Implementierung vom letzten Übungsblatt):
A=(1 1 1
δ 0 0
0 δ 0
0 0 δ)
,b =(3
δ
δ
δ)
δ = 10^-7.
Beim Vergleich beider Verfahren soll der relative Fehler zur exakten Lösung x = (1, 1, 1)^T
in der ℓ2-Norm bestimmt werden.
meine lösung lautet/
import numpy as np
def Householder_Vektor(x):
v = np.copy(x)
v[0] = x[0] + np.sign(x[0]) * np.linalg.norm(x)
v = v / np.linalg.norm(v)
return v
def Householder_Mult(A, v):
return A - 2.0 * np.outer(v, np.dot(v, A))
def Householder_QR(A):
m, n = A.shape
R = np.copy(A)
Q = np.eye(m)
for i in range(min(m-1, n)):
v = Householder_Vektor(R[i:, i])
H = np.eye(m)
H[i:, i:] -= 2.0 * np.outer(v, v)
R = np.dot(H, R)
Q = np.dot(H, Q)
return Q, R
def Householder_Least_Squares(A, b):
Q, R = Householder_QR(A)
# Löse das System R * x = Q^T * b für x
y = np.dot(Q.T, b)[:A.shape[1]]
x_ls = np.linalg.solve(R[:A.shape[1]], y)
# Berechne das Residuum
residuum = np.linalg.norm(np.dot(A, x_ls) - b)
return x_ls, residuum
# Test
delta = 1e-7
A = np.array([[1, 1, 1], [delta, 0, 0], [0, delta, 0], [0, 0, delta]])
b = np.array([3, delta, delta, delta])
# Householder Least Squares
x_ls, residuum_ls = Householder_Least_Squares(A, b)
# Normalengleichung
A_normal_eq = np.dot(A.T, A)
b_normal_eq = np.dot(A.T, b)
x_normal_eq = np.linalg.solve(A_normal_eq, b_normal_eq)
residuum_normal_eq = np.linalg.norm(np.dot(A, x_normal_eq) - b)
# Berechne den relativen Fehler
exact_solution = np.array([1, 1, 1])
relative_error_ls = np.linalg.norm(x_ls - exact_solution) / np.linalg.norm(exact_solution)
relative_error_normal_eq = np.linalg.norm(x_normal_eq - exact_solution) / np.linalg.norm(exact_solution)
print("Householder Least Squares Lösung:", x_ls)
print("Householder Least Squares Residuum:", residuum_ls)
print("Normalengleichung Lösung:", x_normal_eq)
print("Normalengleichung Residuum:", residuum_normal_eq)
print("Relative Fehler (Householder Least Squares):", relative_error_ls)
print("Relative Fehler (Normalengleichung):", relative_error_normal_eq)
kann jemand den Code nachsehen
Danke
(a) Implementieren Sie dazu folgende Methoden:
- Householder_Vektor(x) liefert einen Householder-Vektor zum Vektor x zurück.
- Householder_Mult(A,v) multipliziert die Matrix A mit der Householder-Matrix zum Householder Vektor v und liefert das Produkt zurück
- Householder_QR(A) berechnet eine QR-Zerlegung der Householder-Matrix A bei der der essentielle Teil der Householder-Matrizen unter der Diagonalen von A gespeichert wird und liefert die resultierende Matrix zurück.
- Householder_Least_Squares(A,b) berechnet die kleinste Quadrate Lösung x von Ax=b und liefert x sowie das Residuum ρLS zurück.
Beachte, dass in diesen Algorithmen die Dimensionen der Vektoren und Matrizen variieren. Tipp: Dafür dürfen Kopien von Teilvektoren oder Teilmatrizen geeigneter Dimensionen angelegt werden.
(b) Testen Sie Ihr Verfahren, in dem Sie mithilfe der QR-Zerlegung ein Ausgleichsproblem mit
den folgenden Daten lösen und vergleichen Sie Ihr Ergebnis mit der Lösung, die Sie über
die Normalengleichung berechnen (mit ihrer Implementierung vom letzten Übungsblatt):
A=(1 1 1
δ 0 0
0 δ 0
0 0 δ)
,b =(3
δ
δ
δ)
δ = 10^-7.
Beim Vergleich beider Verfahren soll der relative Fehler zur exakten Lösung x = (1, 1, 1)^T
in der ℓ2-Norm bestimmt werden.
meine lösung lautet/
import numpy as np
def Householder_Vektor(x):
v = np.copy(x)
v[0] = x[0] + np.sign(x[0]) * np.linalg.norm(x)
v = v / np.linalg.norm(v)
return v
def Householder_Mult(A, v):
return A - 2.0 * np.outer(v, np.dot(v, A))
def Householder_QR(A):
m, n = A.shape
R = np.copy(A)
Q = np.eye(m)
for i in range(min(m-1, n)):
v = Householder_Vektor(R[i:, i])
H = np.eye(m)
H[i:, i:] -= 2.0 * np.outer(v, v)
R = np.dot(H, R)
Q = np.dot(H, Q)
return Q, R
def Householder_Least_Squares(A, b):
Q, R = Householder_QR(A)
# Löse das System R * x = Q^T * b für x
y = np.dot(Q.T, b)[:A.shape[1]]
x_ls = np.linalg.solve(R[:A.shape[1]], y)
# Berechne das Residuum
residuum = np.linalg.norm(np.dot(A, x_ls) - b)
return x_ls, residuum
# Test
delta = 1e-7
A = np.array([[1, 1, 1], [delta, 0, 0], [0, delta, 0], [0, 0, delta]])
b = np.array([3, delta, delta, delta])
# Householder Least Squares
x_ls, residuum_ls = Householder_Least_Squares(A, b)
# Normalengleichung
A_normal_eq = np.dot(A.T, A)
b_normal_eq = np.dot(A.T, b)
x_normal_eq = np.linalg.solve(A_normal_eq, b_normal_eq)
residuum_normal_eq = np.linalg.norm(np.dot(A, x_normal_eq) - b)
# Berechne den relativen Fehler
exact_solution = np.array([1, 1, 1])
relative_error_ls = np.linalg.norm(x_ls - exact_solution) / np.linalg.norm(exact_solution)
relative_error_normal_eq = np.linalg.norm(x_normal_eq - exact_solution) / np.linalg.norm(exact_solution)
print("Householder Least Squares Lösung:", x_ls)
print("Householder Least Squares Residuum:", residuum_ls)
print("Normalengleichung Lösung:", x_normal_eq)
print("Normalengleichung Residuum:", residuum_normal_eq)
print("Relative Fehler (Householder Least Squares):", relative_error_ls)
print("Relative Fehler (Normalengleichung):", relative_error_normal_eq)
kann jemand den Code nachsehen
Danke
Diese Frage melden
gefragt
abdull
Punkte: 48
Punkte: 48