Wydajne obliczanie odwrotności pierwiastka kwadratowego z macierzy

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 .

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$

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

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.