Algorithmus für das Gram-Schidtsche Orthogonalisierungsverfahren

Wir können das vorangegangene Beispiel jetzt als Ausgangspunkt für die algorithmische Umsetzung nehmen:

def GramSchmidt(A):
    M = A.shape[0]
    N = A.shape[1]
    assert(M >= N)
    Q = numpy.zeros((M,N))
    for k in range(0,N):
        Q[:,k] = A[:,k]
        for j in range(0,k):
            Q[:,k] = Q[:,k] - numpy.dot(Q[:,j],A[:,k])*Q[:,j]
        Q[:,k] = (1.0 / linalg.norm(Q[:,k])) * Q[:,k]
    return Q

Wir testen es anhand des vorherigen Beispiels, wobei wir jetzt A als die Matrix der Spalten-Vektoren \(a_k\) speichern:

A = numpy.zeros((3,3))
A[:,0] = [0,2,0]
A[:,1] = [8,4,0]
A[:,2] = [10,6,12]
print "Gram Schmidt von"
print A

Q = GramSchmidt(A)
print "liefert"
print Q
print "Q^T*Q"
print numpy.dot(Q,Q)

Varianten

Auch hier können wir wieder die Konstruktion inplace durchführen:

def GramSchmidt_inplace(A):
    M = A.shape[0]
    N = A.shape[1]
    assert(M >= N)
    for k in range(0,N):
        ak = A[:,k]
        for j in range(0,k):
            A[:,k] = A[:,k] - numpy.dot(A[:,j],ak)*A[:,j]
        A[:,k] = (1.0 / linalg.norm(A[:,k])) * A[:,k]
    return A

Q = GramSchmidt_inplace(numpy.copy(A))
print "inplace Variante liefert"
print Q

und mit Hilfe von numpy läßt sich die Summe auch durch entsprechende Matrix-Vektor-Produkte ersetzen:

def GramSchmidt_inplace_compact(A):
    M = A.shape[0]
    N = A.shape[1]
    assert(M >= N)
    for k in range(0,N):
        A[:,k] = A[:,k] - numpy.dot(A[:,:k],numpy.dot(A[:,k],A[:,:k]))
        A[:,k] = (1.0 / linalg.norm(A[:,k])) * A[:,k]
    return A

Q = GramSchmidt_inplace_compact(numpy.copy(A))
print "inplace compact Variante liefert"
print Q

Auch hier ist wieder (wie bei der Cholesky-Zerlegung) das Problem, dass wir spaltenweise durch die Matrix laufen, diese aber meist zeilenweise gespeichert ist. Dies führt zu ungünstigen Speicherzugriffen.