Z czego korzysta LAPACK

Procedura QR LAPACK przechowuje Q jako reflektory Householdera. Skaluje wektor odbicia pomocą , więc pierwszy element wyniku staje się , więc nie trzeba go przechowywać. I przechowuje osobny wektor , który zawiera potrzebne współczynniki skali. Zatem macierz reflektorów jest taka:$v$$1/v_1$$1$$\tau$

gdzie nie jest znormalizowane. Natomiast w podręcznikach matryca reflektorów jest$v$

gdzie jest znormalizowane.$v$

Dlaczego skala LAPACK z , zamiast normalizowania go?$v$$1/v_1$

Potrzebne miejsce jest takie samo (zamiast , należy zapisać ), a następnie zastosowanie można wykonać szybciej, ponieważ nie ma potrzeby mnożenia przez (mnożenie przez w wersji podręcznika można zoptymalizować, jeśli zamiast prostej normalizacji, jest skalowane przez ).$\tau$$v_1$$H$$\tau$$2$$v$$\sqrt 2/\|v\|$

(Powodem mojego pytania jest to, że piszę procedurę QR i SVD i chciałbym poznać przyczynę tej decyzji, czy muszę ją przestrzegać, czy nie)

To zablokowany wariant Householder-QR napędza ten projekt. Jeśli spojrzysz na książkę Goluba i Van Loana (rozdział 5.2 lub więcej), mówią o tym, jak k-iteracje algorytmu mogą być blokowane razem poprzez gromadzenie poszczególnych reflektorów w reflektorze rangi w postaci , gdzie zarówno i są macierzami „o wysokim chuście” o rozmiarze . Ten algorytm działa więcej, ale w praktyce jest szybszy, ponieważ jest bogaty w wywołania gemm (). Niestety, marnowanie pamięci jest marnotrawione z powodu konieczności niezależnego reprezentowania i$\mathbf I + \mathbf W \mathbf Y^{\mathrm T}$$\mathbf W$$\mathbf Y$$n \times k$$\mathbf W$$\mathbf Y$

W późniejszym artykule (cytowanym poniżej) Van Loan opisuje bardziej wydajną „symetryczną” strukturę danych, odbłyśnik blokowy w postaci . Tutaj wciąż jest , ale wymóg flop / storage dla formowania został wyeliminowany poprzez wprowadzenie , małej górnej trójkątnej macierzy. Chociaż potrzeba pomnożenia przez wprowadza niewielką ilość dodatkowej pracy, zazwyczaj jest to zysk netto, ponieważ .$\mathbf I + \mathbf Y \mathbf T \mathbf Y^{\mathrm T}$$\mathbf Y$$n \times k$$\mathbf W$$\mathbf T$$k \times k$$\mathbf T$$k << n$

W LAPACK, niezablokowany algorytm jest tak naprawdę ograniczającym przypadkiem algorytmu blokowego, aż do wyboru symboli (co prowadzi nas do , małej wersji Trójkąt ).$k \rightarrow 1$$\tau$$1\times1$$\mathbf T$

Cytowanie: Schreiber, Robert i Charles Van Loan. „Wydajna dla przechowywania reprezentacja WY dla produktów transformacji Householder”. SIAM Journal on Scientific and Statistics Computing 10.1 (1989): 53-57.

Nie musisz przechowywać $\tau$, możesz ponownie obliczyć go z pozostałej części wektora. (Możesz ponownie obliczyć$v_1$ z innych pozycji również w znormalizowanej wersji, ale jest to wyraźnie niestabilne obliczenie z powodu tych odejmowań).

W rzeczywistości możesz ponownie użyć dolnej trójkątnej części $R$ przechować $v_2,...v_n$, aby faktoryzacja została obliczona w pełni na miejscu. Lapack bardzo dba o te lokalne wersje algorytmów.

Moja sugestia jest oparta na dokumentacji dla Intel MKL https://software.intel.com/en-us/mkl-developer-reference-c-geqrf . Wygląda jak wartości na przekątnej magazynu wyjściowego R i powyżej, więc dla Q pozostaje tylko niższy trójkąt. Naturalnym wydaje się wykorzystanie dodatkowej pamięci dla czynników skalujących.