我想知道如何在矩阵数组上执行快速Cholesky因式分解。假设我有一个3x3x1000个pd矩阵的数组A,并且希望找到这1000个矩阵的cholesky分解。我知道我可以很容易地写出
for n=1:1000
cl(:,:,n) = chol(A(:,:,n);
end
但是,由于for循环,这是一个相当慢的过程。有没有任何方法来加速这一点,并完全避免'for‘循环。任何想法都会受到赞赏。
我的一个想法是将数组转换为块对角线矩阵,并接受它的chol分解,因此将有一个cholesky分解的块矩阵,但是,如果不需要创建一个完整的3000x3000矩阵,我还不知道如何执行这个矩阵(这样的矩阵大小也会使过程变慢)。我不熟悉在matlab上使用稀疏矩阵,这可能是解决方案之一。
发布于 2014-09-11 18:59:55
正如您所描述的那样,通过构造稀疏矩阵确实有可能获得一些计算时间(可能是因为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)所有使用数组行号、冒号、值(所有长度相同)的非零元素,使C( rowinds (k),colinds (K))=值(K),即
values = A(:);
C=sparse(rowinds,colinds,values,D*N,D*N);
问题是如何构造数组行和colinds,在您的示例中,这些数组将看起来像rowinds =1 1 1 2 2 2 3 3 3 4 4 4 5 5 6 6 6 7 7 8 8 8 9.糖尿病肾病和结肠症=1 2 3 1 2 3 1 3 3 4 4 5 6 6 4 5 6 6 6 7 8 8DN-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.
我认为,应该有一种更快的方法来构造指数数组。
https://stackoverflow.com/questions/25787713
复制相似问题