QR-факторизация с использованием модифицированного Грама Шмидта

Вопрос: для этой задачи вам дан список матриц As, и ваша задача — найти QR-факторизацию для каждой из них.

Реализовать qr_by_gram_schmidt: эта функция принимает в качестве входных данных матрицу A и вычисляет разложение QR, возвращая две переменные, Q и R, где A=QR, с ортогональным Q и нулем R ниже диагонали.

A представляет собой матрицу размера n × m, где n ≥ m (т.е. строк больше, чем столбцов).

Вы должны реализовать эту функцию, используя модифицированную процедуру Грама-Шмидта.

ВХОД:

Как: Список массивов

ВЫХОД:

Qs: Список матриц Q, выведенных qr_by_gram_schmidt, в том же порядке, что и As. Для матрицы A формы n×m Q должна иметь форму n×m.
Rs: Список матриц R, выведенных qr_by_gram_schmidt, в том же порядке, что и As. Для матрицы A формы n×m R должен иметь форму m×m

Я написал код для факторизации QR, который я считаю правильным:

import numpy as np
def qr_by_gram_schmidt(A):
m = np.shape(A)[0]
n = np.shape(A)[1]
Q =  np.zeros((m, m))
R =  np.zeros((n, n)) 
for j in xrange(n):
    v = A[:,j]
    for i in xrange(j):
        R[i,j] = Q[:,i].T * A[:,j]
        v = v.squeeze() - (R[i,j] * Q[:,i])
    R[j,j] =  np.linalg.norm(v)
    Q[:,j] = (v / R[j,j]).squeeze()
return Q, R

Как мне написать цикл для вычисления QR-факторизации каждой из матриц в As и сохранения их в этом порядке?

edit: в коде тоже есть ошибка. Буду признателен, если вы поможете мне в его отладке.

Спасибо


person amateur_programmer    schedule 15.04.2016    source источник


Ответы (1)


Я не проверял ваш код GS, но мне пришлось внести изменения (возможно, неправильные!), чтобы он скомпилировался. Вам просто нужно настроить список ваших матриц, я сделал 2 из них, а затем прокрутите этот список и примените вашу функцию.

импортировать numpy как np

def gs(A):
    m = np.shape(A)[0]
    n = np.shape(A)[1]
    Q =  np.zeros((m, m))
    R =  np.zeros((n, n)) 
    print m,n,Q,R
    for j in xrange(n):
        v = A[:,j]
        for i in xrange(j):
            R[i,j] =  np.dot(Q[:,i].T , A[:,j])   # I made an arbitrary change here!!!
            v = v.squeeze() - (R[i,j] * Q[:,i])
        R[j,j] =  np.linalg.norm(v)
        Q[:,j] = (v / R[j,j]).squeeze()
    return Q, R

As= np.random.rand(2,3,3)  # list of 2 (3x3) matrices
print As

for A in As:
    print gs(A)

Выход:

[[[ 0.9599614   0.02213113  0.43343881]
  [ 0.44202415  0.6816688   0.88321052]
  [ 0.93098107  0.80528361  0.88473308]]

 [[ 0.41794678  0.10762796  0.42110659]
  [ 0.89598082  0.81225543  0.52947205]
  [ 0.0621515   0.59826789  0.14021332]]]
(array([[ 0.68158915, -0.67980134,  0.27075149],
       [ 0.31384477,  0.60583989,  0.73106736],
       [ 0.66101262,  0.41331364, -0.626286  ]]), array([[ 1.40841649,  0.76132516,  1.15743793],
       [ 0.        ,  0.73077208,  0.60610414],
       [ 0.        ,  0.        ,  0.20894464]]))
(array([[ 0.42190511, -0.39510208,  0.81602109],
       [ 0.90446656,  0.121136  , -0.40898205],
       [ 0.06274013,  0.91061541,  0.40846452]]), array([[ 0.99061796,  0.81760207,  0.66535379],
       [ 0.        ,  0.6006613 ,  0.02543844],
       [ 0.        ,  0.        ,  0.18435946]]))
person roadrunner66    schedule 15.04.2016
comment
Спасибо за помощь. Но этот код не дает правильного вывода, и я не могу понять, почему. Форма Q и R на каждой итерации неверна - person amateur_programmer; 15.04.2016