Calcul efficace de l'inverse de la racine carrée de la matrice

Un problème courant en statistique est le calcul de l'inverse de la racine carrée d'une matrice définie positive symétrique. Quelle serait la façon la plus efficace de calculer cela?

Je suis tombé sur de la littérature (que je n'ai pas encore lue) et un code R accessoire ici , que je reproduirai ici pour plus de commodité

# 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))
}

Je ne suis pas tout à fait sûr de comprendre la ligne d = (d+abs(d))/2. Existe-t-il un moyen plus efficace de calculer l'inverse de la racine carrée de la matrice? La eigenfonction R appelle LAPACK .

Le code que vous avez publié utilise la décomposition de valeurs propres de la matrice symétrique pour calculer . $A^{-1/2}$

La déclaration

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

prend effectivement toute entrée négative dans d et la met à 0, tout en laissant les entrées non négatives seules. Autrement dit, toute valeur propre négative de est traitée comme si elle était de 0. En théorie, les valeurs propres de A doivent toutes être non négatives, mais en pratique, il est courant de voir de petites valeurs propres négatives lorsque vous calculez les valeurs propres d'un certain supposé positif matrice de covariance presque singulière. $A$

Si vous avez vraiment besoin de l'inverse de la racine carrée de la matrice symétrique de , et que est raisonnablement petit (pas plus grand que disons 1000 par 1000), alors c'est à peu près aussi bon que n'importe quelle méthode que vous pourriez utiliser. $A$$A$

Dans de nombreux cas, vous pouvez utiliser à la place un facteur Cholesky de l'inverse de la matrice de covariance (ou pratiquement le même, le facteur Cholesky de la matrice de covariance elle-même.) Le calcul du facteur Cholesky est généralement un ordre de grandeur plus rapide que le calcul de la décomposition des valeurs propres pour matrices denses et beaucoup plus efficaces (à la fois en temps de calcul et en stockage requis) pour les matrices grandes et clairsemées. Ainsi, l'utilisation de la factorisation de Cholesky devient très souhaitable lorsque est grand et clairsemé. $A$

D'après mon expérience, la méthode polar-newton de Higham fonctionne beaucoup plus rapidement (voir le chapitre 6 des fonctions des matrices de N. Higham). Dans ma courte note, il y a des graphiques qui comparent cette méthode aux méthodes de premier ordre. En outre, des citations de plusieurs autres approches de matrice-racine carrée sont présentées, bien que principalement l'itération polaire de Newton semble fonctionner le mieux (et évite de faire des calculs de vecteurs propres).

% 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

Optimisez votre code:

Option 1 - Optimisez votre code R:
a. Vous pouvez apply()une fonction à la dfois max(d,0)et d2[d==0]=0en une seule boucle.
b. Essayez d'opérer ei$valuesdirectement.

Option 2 - Utilisez C ++:
Réécrivez toute la fonction en C ++ avec RcppArmadillo. Vous pourrez toujours l'appeler depuis R.