Wydajne obliczanie odwrotności pierwiastka kwadratowego z macierzy

15

Częstym problemem w statystyce jest obliczanie pierwiastka kwadratowego odwrotnego symetrycznej dodatniej macierzy określonej. Jaki byłby najbardziej efektywny sposób obliczenia tego?

Natknąłem pewnym literaturze (które nie zostały jeszcze przeczytane) oraz jakiegoś przypadkowego kodu R tutaj , które będę tutaj dla wygody odtworzenia

# function to compute the inverse square root of a matrix
fnMatSqrtInverse = function(mA) {
  ei = eigen(mA)
  d = ei$values
      d = (d+abs(d))/2
      d2 = 1/sqrt(d)
      d2[d == 0] = 0
      return(ei$vectors %*% diag(d2) %*% t(ei$vectors))
}

Nie jestem do końca pewien, czy rozumiem tę linię d = (d+abs(d))/2. Czy istnieje bardziej skuteczny sposób obliczania odwrotności pierwiastka kwadratowego macierzy? Funkcja R eigenwywołuje LAPACK .

linear-algebra numerical-analysis matrix

— tchakravarty
źródło

d

$d$

(d + | d |) / 2

$(d+|d|)/2$

max (d, 0)

$\max(d,0)$

A^{- 1 / 2}

$A^{-1/2}$

A^{- 1 / 2}

$A^{-1/2}$

x

$x$

@DanielShapero Dziękujemy za komentarz. Więc jeśli mam matrycę PSD, nie potrzebuję tej linii? Moja aplikacja wymaga obliczania kwadratowych form, takich jak .

A^{- 1 / 2} B A^{- 1 / 2}

$A^{-1/2}BA^{-1/2}$

— tchakravarty

Nie znam R, ale w wierszu 7 zakładam, że ma on logiczne indeksowanie jak Matlab. Jeśli tak, proponuję przepisać wiersz 5 jako d[d<0] = 0, co jest bardziej wyraziste.

— Federico Poloni

Czy ten kod jest poprawny? Uruchomiłem go na prostym przykładzie w Matlabie i stwierdziłem, że odpowiedź jest błędna. Moja matryca jest pozytywna, ale zdecydowanie nie symetryczna. Zobacz moją odpowiedź poniżej: przesłałem kod do Matlaba.

— roni

10

W opublikowanym kodzie zastosowano rozkład wartości własnej macierzy symetrycznej do obliczenia . $A^{-1/2}$

Wyrok

d = (d + abs (d)) / 2

efektywnie przyjmuje wszelkie ujemne wpisy id i ustawia je na 0, pozostawiając nieujemne wpisy w spokoju. Oznacza to, że każda ujemna wartość własna jest traktowana tak, jakby wynosiła 0. Teoretycznie wartości własne A powinny być nieujemne, ale w praktyce często zdarza się widzieć małe ujemne wartości własne podczas obliczania wartości własnych z rzekomo pozytywnego określonego macierz kowariancji, która jest prawie pojedyncza. $A$

Jeśli naprawdę potrzebujesz odwrotności pierwiastka kwadratowego z macierzy symetrycznej , a jest dość małe (nie większe niż powiedzmy 1000 na 1000), to jest to tak dobra jak każda metoda, której możesz użyć. $A$ $A$

W wielu przypadkach można zamiast tego użyć współczynnika Cholesky'ego odwrotności macierzy kowariancji (lub praktycznie takiego samego, współczynnika Cholesky'ego samej macierzy kowariancji). Obliczanie współczynnika Cholesky'ego jest zwykle o rząd wielkości szybsze niż obliczanie rozkładu wartości własnej dla gęste matryce i znacznie bardziej wydajne (zarówno w czasie obliczeń, jak i wymaganym przechowywaniu) dla dużych i rzadkich matryc. Zatem zastosowanie faktoryzacji Choleskiego staje się bardzo pożądane, gdy jest duże i rzadkie. $A$

— Brian Borchers
źródło

6

Odpowiedź Briana daje dobrą radę: zamiast tego użyj współczynnika Cholesky'ego (jeśli możesz). Nie ma innego optymalizacji można zrobić na początku, że: nie obliczyć swoją PSD macierzy . Często otrzymujesz z obliczeń takich jak , z prostokątemW tym przypadku jest to wystarczające do obliczenia rozkładu QR w celu uzyskania współczynnika Choleskiego z , o znacznie większej dokładności.

A

$A$

A

$A$

A = B^{T} B

$A=B^TB$

B

$B$

B

$B$

R

$R$

A

$A$

— Federico Poloni

5

Z mojego doświadczenia wynika, że metoda Higham-biegunowa-Newtona działa znacznie szybciej (patrz Rozdział 6 Funkcji Matryc autorstwa N. Highama). W tej krótkiej notatce znajdują się wykresy porównujące tę metodę do metod pierwszego rzędu. Przedstawiono również cytowania kilku innych podejść do macierzy-pierwiastka kwadratowego, chociaż najlepiej wydaje się, że najlepiej działa biegunowa iteracja Newtona (i unika się wykonywania obliczeń wektorów własnych).

% compute the matrix square root; modify to compute inverse root.
function X = PolarIter(M,maxit,scal)
  fprintf('Running Polar Newton Iteration\n');
  skip = floor(maxit/10);
  I = eye(size(M));
  n=size(M,1);
  if scal
    tm = trace(M);
    M  = M / tm;
  else
    tm = 1;
  end
  nm = norm(M,'fro');

  % to compute inv(sqrt(M)) make change here
  R=chol(M+5*eps*I);

  % computes the polar decomposition of R
  U=R; k=0;
  while (k < maxit)
    k=k+1;
    % err(k) = norm((R'*U)^2-M,'fro')/nm;
    %if (mod(k,skip)==0)
    %  fprintf('%d: %E\n', k, out.err(k));
    %end

    iU=U\I;
    mu=sqrt(sqrt(norm(iU,1)/norm(U,1)*norm(iU,inf)/norm(U,inf)));
    U=0.5*(mu*U+iU'/mu);

   if (err(k) < 1e-12), break; end
  end
  X=sqrt(tm)*R'*U;
  X = 0.5*(X+X');
end

— suvrit
źródło

0

Zoptymalizuj swój kod:

Opcja 1 - Zoptymalizuj swój kod R:
a. Można apply()funkcja celu d, który będzie zarówno max(d,0)i d2[d==0]=0na jednej pętli.
b. Spróbuj operować ei$valuesbezpośrednio.

Opcja 2 - użyj C ++:
Przepisz całą funkcję w C ++ za pomocą RcppArmadillo. Nadal będziesz mógł zadzwonić do niego z R.

— moc
źródło