Cálculo eficiente da raiz quadrada da matriz inversa

15

Um problema comum em estatística é calcular a raiz quadrada inversa de uma matriz definida positiva simétrica. Qual seria a maneira mais eficiente de calcular isso?

Encontrei alguma literatura (que ainda não li) e algum código R incidental aqui , que reproduzirei aqui por conveniência

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

Não tenho muita certeza de entender a linha d = (d+abs(d))/2. Existe uma maneira mais eficiente de calcular a raiz quadrada da matriz inversa? A eigenfunção R chama LAPACK .

linear-algebra numerical-analysis matrix tchakravarty
fonte
Desde que seja real, é igual a . Isso efetivamente elimina quaisquer autovalores negativos que a matriz possa ter. Você precisa de todas as entradas da matriz A ^ {- 1/2} ou é suficiente para poder multiplicar A ^ {- 1/2} por um vetor arbitrário x ? d(d+|d|)/2max(d,0) A - 1 / 2 xA−1/2A−1/2x
precisa
Obrigado por seu comentário. Então, se eu tenho uma matriz PSD, não preciso dessa linha? Meu aplicativo requer o cálculo de formas quadráticas como A−1/2BA−1/2 .
Tchakravarty
Não estou familiarizado com R, mas, dada a linha 7, presumo que ele tenha indexação lógica como o Matlab. Nesse caso, sugiro que você reescreva a linha 5 como d[d<0] = 0, que é mais expressiva.
Federico Poloni
Esse código está correto? Eu o executei em um exemplo simples no matlab e achei a resposta errada. Minha matriz é positiva definida, mas definitivamente não simétrica. Por favor, veja minha resposta abaixo: Transferi o código para o matlab.
Roni 12/05

Respostas:

10

O código que você postou usa a decomposição de autovalor da matriz simétrica para calcular . A−1/2

A declaração

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

efetivamente pega qualquer entrada negativa em d e a define como 0, deixando apenas as entradas não negativas. Ou seja, qualquer autovalor negativo de é tratado como se fosse 0. Na teoria, os autovalores de A devem ser todos não negativos, mas na prática é comum ver autovalores negativos pequenos quando você calcula os autovalores de um valor definido supostamente positivo. matriz de covariância quase singular. A

Se você realmente precisa do inverso da raiz quadrada da matriz simétrica de , e é razoavelmente pequeno (não maior que 1.000 por 1.000), então isso é tão bom quanto qualquer método que você possa usar. AAA

Em muitos casos, você pode usar um fator de Cholesky do inverso da matriz de covariância (ou praticamente o mesmo, o fator de Cholesky da própria matriz de covariância.) A computação do fator de Cholesky é tipicamente uma ordem de magnitude mais rápida que a decomposição de matrizes densas e muito mais eficientes (tanto em tempo computacional quanto no armazenamento necessário) para matrizes grandes e esparsas. Assim, o uso da fatoração de Cholesky se torna muito desejável quando é grande e escasso. A

Brian Borchers
fonte
6
A resposta de Brian dá bons conselhos: use o fator Cholesky (se puder). Há uma outra otimização que você pode fazer em cima disso: não calcular a sua matriz PSD . Frequentemente, você obtém de um cálculo como , com retangular. Nesse caso, é suficiente calcular uma decomposição QR de para obter o fator de Cholesky de , com uma precisão muito melhor. A A = B T B B B R AAAA=BTBBBRA
Federico Poloni
5

Na minha experiência, o método polar-Newton de Higham funciona muito mais rápido (ver Capítulo 6 de Funções de matrizes por N. Higham). Em esta pequena nota da mina há parcelas que comparam esse método para métodos de primeira ordem. Além disso, são apresentadas citações para várias outras abordagens de matriz quadrada-raiz, embora a iteração de Newton polar pareça funcionar melhor (e evita cálculos de vetor próprio).

% 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
fonte
0

Otimize seu código:

Opção 1 - Otimize seu código R:
a. Você pode apply()uma função para dque ambos max(d,0)e d2[d==0]=0em um loop.
b. Tente operar ei$valuesdiretamente.

Opção 2 - Use C ++:
Reescreva toda a função em C ++ com RcppArmadillo. Você ainda poderá chamá-lo de R.

poder
fonte