Частично схема Грама-Шмидта на большом наборе векторов

Question

В моделировании я создаю большой набор векторов, к которым применяю схему схемы Грамма-Шмидта. Под частичным я подразумеваю, что сначала я удаляю компоненты векторов, а затем вычисляю их длину (проверьте формулы orth_2 или orth_3 в приведенном ниже коде). Затем длина сохраняется для вычисления желаемого значения (l-значения, см. Код ниже). Только после этого я нормализую вектор.

Выше объясненная процедура показана в следующем коде для трех векторов:

from numpy import random

N = 3                  'Numer of rows and columns. In my real case N is equal to 100'

l1, l2, l3 = 0, 0, 0   'Quantities that i want to compute. Again in my case i have 100s of them l1, l2,...,l100' 


for t in range(10**3):  'Simulate Computation' 

    ''' for i in range(10**5):

        Computation where i compute Matrix A with its vectors A[0], A[1]... '''

    A = random.rand(N, N) 'A is generated randomly to keep the code short'

   'Apply the Gram-Schmidt-Scheme on the vectors A[0], A[1], A[2]. Again in my case i have 100 vectors'

    orth_1 = A[0]                   
    l1    += log(norm(orth_1)) 
    A[0]   = orth_1 / norm(orth_1)        

    orth_2 = A[1] - dot(A[1], A[0]) * A[0] 'Remove any A[0] components from A[1]'
    l2    += log(norm(orth_2))             'Store the values that i want to compute later'
    A[1]   = orth_2 / norm(orth_2)           

    orth_3 = A[2] - (dot(A[2], A[1]) * A[1]) - (dot(A[2], A[0]) * A[0]) 'Remove any A[0] and A[1] components from A[2]'
    l3    += log(norm(orth_3)) 
    A[2]   = orth_3 / norm(orth_3)

Как вы можете видеть формулы для orth_ становиться все длиннее и длиннее по мере увеличения числа векторов. Для N = 3 все еще выполнимо, но для N = 100 становится сложным выписать все формулы для orth_1, orth_2, .., orth_100 для вычисления желаемых значений l1, l2, ..., l100.

У кого-нибудь есть идеи, как я могу распылять этот процесс?

Answer 1

Я скопировал ваш исходный код, чтобы я мог сравнить ответы, чтобы убедиться, что математика сработала. Я считаю, что эта функция возвращает тот же ответ, который вычисляет ваш код.

def removeOrthos(Arr, N, l):
    remove = 0
    for i in range(N):
        remove -= np.dot(Arr[N], Arr[i]) * Arr[i]
    ortho = Arr[N] + remove
    l[N] += np.log(np.linalg.norm(ortho))
    return ortho / np.linalg.norm(ortho)

Вот ваш оригинальный код со сравнениями с моей функцией и linalg.qr ()

def ortho():
    N = 3 #                  'Numer of rows and columns. In my real case N is equal to 100'
    l1, l2, l3 = 0, 0, 0  # 'Quantities that i want to compute. Again in my case i have 100s of them l1, l2,...,l100' 
    L = np.zeros(N)
    for t in range(10**3):  #'Simulate Computation' 

#''' for i in range(10**5):
#Computation where i compute Matrix A with its vectors A[0], A[1]... '''

        A = np.random.rand(N, N) #'A is generated randomly to keep the code short'

#'Apply the Gram-Schmidt-Scheme on the vectors A[0], A[1], A[2]. Again in my case i have 100 vectors'

        B = np.copy(A)
        for b in range(N):
            B[b] = removeOrthos(B, b, L)

        Q,R = np.linalg.qr(np.copy(A).T)    
        orth_1 = A[0]                   
        l1    += np.log(np.linalg.norm(orth_1)) 
        A[0]   = orth_1 / np.linalg.norm(orth_1)        

        orth_2 = A[1] - np.dot(A[1], A[0]) * A[0] #'Remove any A[0] components from A[1]'
        l2    += np.log(np.linalg.norm(orth_2)) #             'Store the values that i want to compute later'
        A[1]   = orth_2 / np.linalg.norm(orth_2)           

        orth_3 = A[2] - (np.dot(A[2], A[1]) * A[1]) - (np.dot(A[2], A[0]) * A[0]) #'Remove any A[0] and A[1] components from A[2]'
        l3    += np.log(np.linalg.norm(orth_3)) 
        A[2]   = orth_3 / np.linalg.norm(orth_3)
    print(l1,l2,l3, L)

    return A, B, Q, R

Ниже обратите внимание на условие False в [2,0], хотя реальные цифры выглядят как совпадение, должна быть некоторая ошибка округления за пределами 7 знаков после запятой.

>>> A, B, Q, R = ortho()
-74.2981307378 -838.247442533 -1628.34240096 [  -74.29813074  -838.24744253 -1628.34240096]
>>> A==B
array([[ True,  True,  True],
       [ True,  True,  True],
       [False,  True,  True]], dtype=bool)
>>> A
array([[ 0.01225148,  0.554242  ,  0.8322654 ],
       [ 0.99376938, -0.09896161,  0.05127395],
       [ 0.1107805 ,  0.82645169, -0.55200116]])
>>> B
array([[ 0.01225148,  0.554242  ,  0.8322654 ],
       [ 0.99376938, -0.09896161,  0.05127395],
       [ 0.1107805 ,  0.82645169, -0.55200116]])
>>> A==Q
array([[False, False, False],
       [False, False, False],
       [False, False, False]], dtype=bool)
>>> Q
array([[-0.01225148,  0.99376938, -0.1107805 ],
       [-0.554242  , -0.09896161, -0.82645169],
       [-0.8322654 ,  0.05127395,  0.55200116]])

Answer 2

Существует довольно простое решение: примените QR-разложение, которое уже присутствует в numpy.linalg (или scipy.linalg)

Q,R = np.linalg.qr(A.T);

Затем в результате

A=R.T*Q.T

столбцы Q или, как вы его используете, строки Q.T содержат ортонормированные приближенные собственные векторы, а диагональ R или R.T=A*Q содержит сокращенную длину, возможно, со знаком. Чтобы получить нормализованное QR-разложение с положительной диагональю в R, используйте

S = np.diag(np.sign(np.diag(R))); Q, R = Q.dot(S), S.dot(R); 

numpy, используя QR-процедуры LAPACK, которые go через более стабильные отражатели Householder. Но результат до знаков / ориентаций одинаков.