Z czego korzysta LAPACK

9

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$

H = I - τ v v^{T},

$H=I-\tau v v^T,$

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

H = I - 2 v v^{T},

$H = I-2vv^T,$

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)

linear-algebra matrix lapack

— geza
źródło

7

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.

— rchilton1980
źródło

Dziękuję za odpowiedź! Nie rozumiem tego

τ

$\tau$ jest po prostu

1 \times 1

$1 \times 1$ rozmiar

T

$\mathbf T$ . W cytowanej pracy, w algorytmie 5,

Y

$\mathbf Y$ jest

v

$v$ , i

T

$\mathbf T$ wynosi -2. Tak więc kończy się ona wersją podręcznika, a nie wersją LAPACK. Czy coś mi umknęło?

— geza

2

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.

— Federico Poloni
źródło

1

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.

— VorKir
źródło