How to compute cholesky to all slice of a tensor?

1 view (last 30 days)
Hi All,
I was looking for a fast and efficient way to compute the cholesky decomposition to all faces of a tensor and collect them into another tensor.
The easiest approach would be the following one, but I would like to know if it exist a more "elegant" way to do it:
A=rand(N,N,M) % suppose each NxN matrix to be positive semi definite
B=zeros(N,N,M)
for k=1:M
B(:,:,k)=chol(A(:,:,k));
end
Thanks in advance

Accepted Answer

Bruno Luong
Bruno Luong on 25 Oct 2021
For real semi definite matrices, you can take a look at https://www.mathworks.com/matlabcentral/fileexchange/37515-mmx
Unfortunately this package doesn't support complex matrices.
  2 Comments
Bruno Luong
Bruno Luong on 27 Oct 2021
I bench mark pagesvd and mmx('chol',...), on my computer mmx is much faster (8 times)
N=5;
M=100000;
% Generate dummy CorrTensor for testing
R=randn([N,N,M]);
CorrTensor =pagemtimes(R,"none",R,"ctranspose");
% Factorization T(:,:,p)'*T(:,:,p) == CorrTensor(:,:,p)
% using pagesvd
tic
[U,S,V]=pagesvd(CorrTensor,'vector');
T = sqrt(S).*pagectranspose(U); % T(:,:,p)'*T(:,:,p) == CorrTensor(:,:,p)
toc % Elapsed time is 0.156520 seconds.
% using mmx
tic
R = mmx('chol', CorrTensor, []);
toc % Elapsed time is 0.019043 seconds.

Sign in to comment.

More Answers (1)

Christine Tobler
Christine Tobler on 26 Oct 2021
There isn't another way to do this right now. We have functions that do other linear algebra operations to each page of an ND-array: pagemtimes, pagesvd.
Could you say some more about where the matrix A is coming from in your application, and what your next steps are? Perhaps pagemtimes could be useful for working with the result of those Cholesky factorizations?
  2 Comments
Bruno Luong
Bruno Luong on 27 Oct 2021
You can use pagesvd to factorize your Correlation matrix (I suspect this is slower than calling chol in a for-loop):
K=5;
N=3;
M=2;
% Generate dummy CorrTensor for testing
R=randn([N,N,M]);
CorrTensor =pagemtimes(R,"ctranspose",R,"none");
% Factorization T(:,:,p)'*T(:,:,p) == CorrTensor(:,:,p)
[U,S,~]=pagesvd(CorrTensor,'vector');
T = sqrt(S).*pagectranspose(U); % T(:,:,p)'*T(:,:,p) == CorrTensor(:,:,p)
% Generate random correlated normal vectors
e = randn(K,N,M);
e = pagemtimes(e, T),
e =
e(:,:,1) = 1.3038 1.5687 -0.3056 0.7351 0.6160 -1.3979 -1.3023 -2.4529 0.4638 -0.7750 -0.7354 1.1630 0.5316 0.4257 -0.1150 e(:,:,2) = -1.4826 0.3785 3.7542 -0.5790 -0.3740 0.3464 -0.0701 -0.4828 -1.1557 0.0878 1.5143 2.1562 -0.9627 -0.4035 1.0495

Sign in to comment.

Categories

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!