统计学中的一个常见问题是计算对称正定矩阵的平方根逆。计算这个最有效的方法是什么？
我在这里遇到了一些文献（我还没有读过）和一些附带的 R 代码，为了方便我将在这里复制
# 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))
}
我不完全确定我理解这条线d = (d+abs(d))/2。有没有更有效的方法来计算矩阵平方根逆？Reigen函数调用 LAPACK。
您发布的代码使用对称矩阵的特征值分解来计算。
该声明
d=(d+abs(d))/2
有效地接受 d 中的任何负条目并将其设置为 0，同时不理会非负条目。也就是说，的任何负特征值都被视为 0。理论上，A 的特征值都应该是非负的，但在实践中，当您计算假定为正定的特征值时，通常会看到小的负特征值几乎是奇异的协方差矩阵。
的对称矩阵平方根的逆，并且相当小（不大于 1,000 x 1,000），那么这与您可能使用的任何方法一样好。
在许多情况下，您可以改用协方差矩阵的逆矩阵的 Cholesky 因子（或实际上相同的协方差矩阵本身的 Cholesky 因子。）计算 Cholesky 因子通常比计算特征值分解快一个数量级密集矩阵，并且对于大型和稀疏矩阵更有效（在计算时间和所需存储方面）。大且稀疏 时，使用 Cholesky 分解变得非常可取。
根据我的经验，Higham 的极牛顿法工作得更快（参见N. Higham的矩阵函数的第 6 章）。在我的这篇简短笔记中，有一些图表将此方法与一阶方法进行了比较。此外，还引用了其他几种矩阵平方根方法，尽管大多数极坐标牛顿迭代似乎效果最好（并且避免进行特征向量计算）。
% 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
优化您的代码：
选项 1 - 优化您的 R 代码
：您可以在一个循环中
同时使用apply()一个函数。湾。尝试直接操作。dmax(d,0)d2[d==0]=0ei$values
选项 2 - 使用 C++：
用RcppArmadillo. 您仍然可以从 R 调用它。