Szybkie obliczanie / szacowanie niskiego rzędu systemu liniowego

Liniowe układy równań są wszechobecne w statystyce obliczeniowej. Jednym specjalnym systemem, z którym się zetknąłem (np. W analizie czynnikowej) jest system

A x = b

$Ax=b$

gdzie Tutaj jest macierzą diagonalną ze ściśle dodatnią przekątną, jest (z ) symetryczną dodatnią półokreśloną macierzą, a jest arbitralną macierz . Jesteśmy proszeni o rozwiązanie ukośnego układu liniowego (łatwego), który został zakłócony przez matrycę niskiej rangi. Naiwnym sposobem rozwiązania powyższego problemu jest odwrócenie za pomocą wzoru Woodbury'ego

A = D + B Ω B^{T}

$A=D+ B \Omega B^T$

D

$D$

n \times n

$n\times n$

Ω

$\Omega$

m \times m

$m\times m$

m ≪ n

$m\ll n$

B

$B$

n \times m

$n\times m$

A

$A$ . Nie wydaje się to jednak właściwe, ponieważ faktoryzacje Cholesky'ego i QR zwykle mogą znacznie przyspieszyć rozwiązanie układów liniowych (i równań normalnych). Niedawno wpadłem na następujący artykuł , który wydaje się przyjmować podejście Choleskiego, i wspomina o niestabilności numerycznej inwersji Woodbury'ego. Artykuł wydaje się jednak w formie szkicu i nie mogłem znaleźć eksperymentów numerycznych ani badań wspierających. Jaki jest stan techniki w rozwiązaniu opisanego przeze mnie problemu?

— niezadowolony
źródło

Ω^{- 1} + B D^{- 1} B^{T}

$\Omega^{-1}+BD^{-1}B^T$

Ω^{- 1}

$\Omega^{-1}$

m << n

$m<<n$

A

$A$

ϵ

$\epsilon$

\bar{B} = D^{- 1 / 2} B

$\bar{B} = D^{-1/2} B$

(I + \bar{B} Ω {\bar{B}}^{T}) x = \bar{b}

$(I + \bar{B}\Omega \bar{B}^T)x = \bar{b}$

\bar{b} = D^{- 1 / 2} b

$\bar{b} = D^{-1/2} b$

Σ = \bar{B} Ω {\bar{B}}^{T}

$\Sigma = \bar{B} \Omega \bar{B}^T$

Σ

$\Sigma$

m

$m$

m ≪ n

$m \ll n$

m

$m$

x = Q (I + Λ)^{- 1} Q^{T} \bar{b}

$x = Q(I + \Lambda)^{-1} Q^T \bar{b}$

Σ = Q Λ Q^{T}

$\Sigma = Q \Lambda Q^T$

Σ

$\Sigma$

(I + \bar{B} Ω {\bar{B}}^{T}) D^{1 / 2} x = \bar{b}

$(I + \bar{B} \Omega \bar{B}^T) D^{1/2} x = \bar{b}$

x = D^{- 1 / 2} Q (I + Λ)^{- 1} Q^{T} D^{- 1 / 2} b

$x = D^{-1/2} Q (I + \Lambda)^{-1} Q^T D^{-1/2} b$

D^{1 / 2}

$D^{1/2}$

x

$x$ w obu przypadkach.) Zauważ, że wszystkie odwrotności mają matryce diagonalne, a zatem są trywialne.

— kardynał

Ω^{- 1} + B^{T} D^{- 1} B

$\Omega^{-1} + B^T D^{-1} B$

„Obliczenia macierzy” Goluba i van Loana szczegółowo omówiono w rozdziale 12.5.1 na temat aktualizacji faktoryzacji QR i Cholesky'ego po aktualizacjach rangi p.

— Fabian Pedregosa
źródło

Wiem, a odpowiednie funkcje lapack są wymienione zarówno w dokumencie, do którego linkowałem, jak i w książce. Zastanawiam się jednak, jaka jest najlepsza praktyka dla danego problemu, a nie dla ogólnego problemu aktualizacji.

— gappy