Klassisches und modifiziertes Gram-Schmidt-Verfahren

Hier werden die Plots aus den Folien zur Wiederholung der Linearen Algebra (LinAWdh.pdf) erstellt. Gleichzeitig sind hier die Musterlösungen zur Aufgabe 41. Diese Lektion wurde nicht in der Vorlesung gezeigt.

In [1]:
import numpy as np;
from scipy.linalg import norm;
from matplotlib import pyplot as plt;

floatFormatter = lambda x: "{0: .2f}".format(x);
np.set_printoptions(formatter = { 'float_kind': floatFormatter });

M = 5;
N = M - 1;
A = np.random.rand(M, N);

Klassisches Gram-Schmidt-Verfahren

Version 1

In [2]:
def GramSchmidt_V1(A):
    Q = np.zeros_like(A) # Float-Nullen der gleichen Größe, wie A
    R = np.zeros((A.shape[1], A.shape[1])) # Platz für die obere Dreiecksmatrix
    for k in range(A.shape[1]):
        Q[:, k] = A[:, k]  
        if k != 0: # Verfahren aus der Vorlesung
            for j in range(k):
                R[j, k] = Q[:, j].T@Q[:, k]
            Q[:, k] = Q[:, k] - Q[:, :k]@R[:k, k]
        R[k, k] = np.sqrt(Q[:, k]@(Q[:, k].T))
        Q[:, k] = Q[:, k]/R[k, k]
    return Q, R

Version 2 (vektorisiert)

In [3]:
def GramSchmidt_V2(A):
    m, n = A.shape;
    Q = 0. * A;        # Float-Nullen der gleichen Größe, wie A
    R = 0. * A[:n, :]; # Platz für die obere Dreiecksmatrix
    for kk in range(n):
        a_k = A[:, kk];
        # Orthogonalisierungsfaktoren
        # der letzte ist Null, da Q in Spalte j noch Nullen enthält
        r_k = a_k @ Q[:, :kk+1];
        # Orthogonalisierung durchführen (alles in einem Schritt)
        q_k = a_k - Q[:, :kk] @ r_k[:-1];
        # Normieren
        r_k[-1] = norm(q_k);
        # Eintragen
        Q[:, kk] = q_k / r_k[-1];
        R[:kk+1, kk] = r_k;
    return Q, R;

Vergleich

In [4]:
Q_, R_ = GramSchmidt_V1(A);
Q, R = GramSchmidt_V2(A);
print("Q =\n{0}\n\nR =\n{1}\n".format(Q, R));
print("Version 1 und Version 2 stimmen ueberein",np.allclose(Q_,Q) and np.allclose(R_,R))
print("max(| A - Q*R |) = {0:g}, max(| Q^T*Q - Id |) = {1:g}".format(
        abs(A - Q@R).max(), abs(Q.T @ Q - np.eye(N)).max()));
Q =
[[ 0.51  0.79 -0.17  0.27]
 [ 0.54 -0.17 -0.33 -0.75]
 [ 0.52 -0.59 -0.07  0.56]
 [ 0.16  0.06  0.83 -0.22]
 [ 0.40 -0.03  0.41  0.04]]

R =
[[ 1.68  1.24  0.62  1.47]
 [ 0.00  0.43  0.02  0.08]
 [ 0.00  0.00  0.64  0.39]
 [ 0.00  0.00  0.00  0.47]]

Version 1 und Version 2 stimmen ueberein True
max(| A - Q*R |) = 1.11022e-16, max(| Q^T*Q - Id |) = 5.56477e-16

Modifiziertes Gram-Schmidt-Verfahren

In [5]:
M = 5;
N = M - 1;
A = np.random.rand(M, N);

Version 1

In [6]:
def GramSchmidtMod_V1(A):
    Q = np.zeros_like(A)
    R = np.zeros((A.shape[1],A.shape[1]))
    for k in range(A.shape[1]):
        Q[:, k] = A[:, k]
        for j in range(k):
            R[j, k] = Q[:, j].T@Q[:, k]
            Q[:, k] -= Q[:,j]*R[j, k]
        R[k, k] = np.sqrt(Q[:, k]@(Q[:, k].T))
        Q[:, k] = Q[:,k]/R[k, k]
    return Q,R

Version 2

In [7]:
def GramSchmidtMod_V2(A):
    m, n = A.shape;
    Q = 0 * A;        # Nullen des gleichen Datentyps, wie A
    R = 0 * A[:n, :];
    for kk in range(n):
        q_k = A[:, kk].copy();
        for jj in range(kk):
            R[jj, kk] = q_k @ Q[:, jj];
            q_k -= R[jj, kk] * Q[:, jj];
        R[kk, kk] = norm(q_k);
        # Eintragen
        Q[:, kk] = q_k / R[kk, kk];
    return Q, R

Vergleich

In [8]:
Q, R = GramSchmidtMod_V1(A);
Q_, R_ = GramSchmidtMod_V2(A);
print("Q =\n{0}\n\nR =\n{1}\n".format(Q, R));
print("Version 1 und Version 2 stimmen ueberein",np.allclose(Q_,Q) and np.allclose(R_,R))
print("max(| A - Q*R |) = {0:g}, max(| Q^T*Q - Id |) = {1:g}".format(
        abs(A - Q@R).max(), abs(Q.T @ Q - np.eye(N)).max()));
Q =
[[ 0.21 -0.23 -0.10  0.92]
 [ 0.59  0.72 -0.10  0.10]
 [ 0.50 -0.01 -0.03 -0.30]
 [ 0.58 -0.65  0.01 -0.20]
 [ 0.10  0.06  0.99  0.10]]

R =
[[ 1.03  1.03  0.58  0.97]
 [ 0.00  0.43 -0.07 -0.50]
 [ 0.00  0.00  0.77  0.77]
 [ 0.00  0.00  0.00  0.66]]

Version 1 und Version 2 stimmen ueberein True
max(| A - Q*R |) = 2.77556e-17, max(| Q^T*Q - Id |) = 3.69321e-16

Vergleich von klassischem und modifiziertem QR

In [9]:
M = 400;
N = M;
A = np.random.rand(M, N);  # Zufallsmatrix
A1 = A.astype(np.float32); # Einfache Genauigkeit
abs(A1 - A).max()

Q1, R1 = GramSchmidt_V1(A1);
print("max(| A - Q1*R1 |) = {0:g}, max(| Q1^T * Q1 - Id |) = {1:g}".format(
        abs(A1 - Q1@R1).max(), (Q1.T @ Q1 - np.eye(N)).max()));

Q2, R2 = GramSchmidtMod_V1(A1);
print("max(| A - Q2*R2 |) = {0:g}, max(| Q2^T * Q2 - Id |) = {1:g}".format(
        abs(A1 - Q2@R2).max(), (Q2.T @ Q2 - np.eye(N)).max()));

        
res1 = np.triu(Q1.transpose() @ Q1 - np.eye(N)).max(1);
res2 = np.triu(Q2.transpose() @ Q2 - np.eye(N)).max(1);


plt.figure(1)
plt.semilogy(res1, ".", label = "klassisch");
plt.semilogy(res2, "x", label = "modifiziert")
plt.plot([ 0, N ], [ 1e-8, 1e-8 ]); # Maschinengenauigkeit (zur Referenz)
plt.ylim([ 8e-9, 1 ]);
plt.legend(numpoints = 1, fancybox = True);
plt.xlabel("k");
plt.savefig("gramschmidt.pdf");

#%%

M = 40;
N = M;
A = np.random.rand(M, N);
A1 = A.astype(np.float32);
abs(A1 - A).max()

Q1, R1 = GramSchmidt_V1(A1);
Q2, R2 = GramSchmidtMod_V2(A1);

orth1 = abs(Q1.T @ Q1);
orth1[orth1 == 0] = 1e-8;
orthLog1 = np.log10(orth1);
orth2 = abs(Q2.T @ Q2);
orth2[orth2 == 0] = 1e-8;
orthLog2 = np.log10(orth2);

plt.figure(2)
plt.subplot(1, 2, 1);
h = plt.imshow(orthLog1, cmap = "gray", interpolation = "nearest");
h.set_clim(vmin = -8, vmax = 0);
plt.title("Klassisch")
plt.subplot(1, 2, 2);
h = plt.imshow(orthLog2, cmap = "gray", interpolation = "nearest");
h.set_clim(vmin = -8, vmax = 0);
plt.title("Modifiziert");
plt.savefig("gramschmidt-orth.pdf");
max(| A - Q1*R1 |) = 6.81743e-08, max(| Q1^T * Q1 - Id |) = 0.00303227
max(| A - Q2*R2 |) = 8.34548e-07, max(| Q2^T * Q2 - Id |) = 5.36689e-05