快速胆甾醇分解矩阵阵列

Question

我想知道如何对矩阵阵列执行快速Cholesky分解。假设我有一个3x3x1000阵列A的pd矩阵，并希望找到这1000个矩阵的cholesky因式分解。我知道我可以很容易地写出类似于

for n=1:1000 
    cl(:,:,n) = chol(A(:,:,n); 
end 

但是，由于for循环，这是一个相当缓慢的过程。有什么方法可以加速这一点，并完全避免“for”循环。任何想法都表示赞赏。

我的一个想法是将数组转换为块对角矩阵并对其进行chol因式分解，因此会产生一个块分解矩阵，但我不知道如何执行此操作而无需创建一个完整的3000x3000矩阵（和这样的矩阵大小可以使程序也慢）。我不熟悉在matlab上使用稀疏矩阵，也许这可能是解决方案之一。

Answer 1

看来确实是有由像你描述构建稀疏矩阵获得一些计算时间的潜力（可能的，因为MATLAB理解这种矩阵的很好的块结构）：

N=10000; 
D=3; 
A=zeros(D,D,N); 
for n=1:N, a=rand(D); A(:,:,n)=a'*a; end 

tic 
cl=zeros(D,D,N); 
for n=1:N, cl(:,:,n)=chol(A(:,:,n)); end 
toc 
% Elapsed time is 0.0919061 seconds. 

B=zeros(D*N); 
for n=1:N, B(D*(n-1)+1:D*(n-1)+D,D*(n-1)+1:D*(n-1)+D)=A(:,:,n); end 
C=sparse(B); 
tic, cl2 = chol(C); toc 
% Elapsed time is 0.006457 seconds. 

现在的问题是如何构造稀疏矩阵C.上面描述的方法非常缓慢，如果D和N较大（我猜，在你的情况下，D和N比你的例子大），需要大量的内存。一般来说，要在matlab中构造一个稀疏矩阵，你必须指定（1）它的维数，即在你的情况下DN x DN（2）所有非零元素使用数组rowinds，colinds，values（所有长度相同），使得C （rowinds（K），colinds（K））=值（K），即

values = A(:); 
C=sparse(rowinds,colinds,values,D*N,D*N); 

的问题是你如何构建阵列rowinds和colinds，而你的情况看起来像rowinds = [1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8 9 9 9 ... DN DN DN]和colinds = [1 2 3 1 2 3 1 2 3 4 5 6 4 5 6 4 5 6 7 8 9 7 8 9 ... DN-2 DN-1 DN]。这是一个可能的解决方案

tic, 
temp = repmat(1:D*N,D,1); rowinds = temp(:); 
cols = zeros(D,D*N); 
for j=1:D 
    tempj = repmat(j:D:D*N-D+j,D,1); cols(j,:)=tempj(:)'; 
end 
colinds=cols(:); 
C=sparse(rowinds,colinds,A(:),D*N,D*N); 
toc 
% Elapsed time is 0.017159 seconds. 
tic, 
cl2 = chol(C); 
toc 
% Elapsed time is 0.006926 seconds. 

我相信，应该有一个更快的方式来构建索引数组。