Perhitungan efisien dari invers akar kuadrat matriks

Masalah umum dalam statistik adalah menghitung kebalikan akar kuadrat dari matriks pasti positif simetris. Apa cara komputasi yang paling efisien ini?

Saya menemukan beberapa literatur (yang belum saya baca), dan beberapa kode R insidental di sini , yang saya akan mereproduksi di sini untuk kenyamanan

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

Saya tidak sepenuhnya yakin saya mengerti garisnya d = (d+abs(d))/2. Apakah ada cara yang lebih efisien untuk menghitung invers kuadrat matriks kuadrat? Fungsi R eigenmemanggil LAPACK .

Kode yang Anda telah diposting menggunakan dekomposisi eigenvalue dari matriks simetrik untuk menghitung . $A^{-1/2}$

Pernyataan

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

secara efektif mengambil setiap entri negatif dalam d dan menetapkannya ke 0, sambil meninggalkan entri non-negatif sendirian. Artinya, setiap nilai eigen negatif $A$ diperlakukan seolah-olah itu adalah 0. Secara teori, nilai eigen dari A semua harus non-negatif, tetapi dalam praktiknya adalah umum untuk melihat nilai eigen negatif kecil ketika Anda menghitung nilai eigen dari yang seharusnya positif pasti matriks kovarians yang hampir tunggal.

Jika Anda benar-benar membutuhkan kebalikan dari akar kuadrat matriks simetris dari , dan $A$$A$ cukup kecil (tidak lebih besar dari katakanlah 1.000 oleh 1.000), maka ini sama baiknya dengan metode apa pun yang mungkin Anda gunakan.

Dalam banyak kasus, Anda dapat menggunakan faktor Cholesky dari kebalikan dari matriks kovarians (atau secara praktis sama, faktor Cholesky dari matriks kovarians itu sendiri.) Menghitung faktor Cholesky biasanya urutan besarnya lebih cepat daripada menghitung dekomposisi nilai eigen untuk matriks padat dan jauh lebih efisien (baik dalam waktu komputasi dan penyimpanan yang diperlukan) untuk matriks besar dan jarang. Dengan demikian menggunakan faktorisasi Cholesky menjadi sangat diinginkan ketika besar dan jarang. $A$

Dalam pengalaman saya, metode Higham kutub-Newton bekerja jauh lebih cepat (lihat Bab 6 dari Fungsi Matriks oleh N. Higham). Dalam catatan singkat saya ini ada plot yang membandingkan metode ini dengan metode urutan pertama. Juga, kutipan untuk beberapa pendekatan matriks-kuadrat-lain disajikan, meskipun sebagian besar iterasi Newton polar tampaknya bekerja paling baik (dan menghindari melakukan perhitungan vektor eigen).

% 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

Optimalkan kode Anda:

Opsi 1 - Optimalkan kode R Anda:
a. Anda bisa apply()menjalankan fungsinya dbaik dalam satu lingkaran max(d,0)maupun d2[d==0]=0dalam satu lingkaran.
b. Coba operasikan ei$valueslangsung.

Opsi 2 - Gunakan C ++: Tulis
ulang seluruh fungsi dalam C ++ dengan RcppArmadillo. Anda masih dapat memanggilnya dari R.