Clear Filters
Clear Filters

How to compute cholesky to all slice of a tensor?

2 views (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
Edited: Bruno Luong on 27 Oct 2021
Using mmx to generate correlated normal random variable as you state in the comment
K=5;
N=3;
M=2;
% Generate dummy CorrTensor for testing
R=randn([N,N,M]);
CorrTensor = pagemtimes(R,"ctranspose",R,"none");
% Factorization R(:,:,p)*R(:,:,p)' == CorrTensor(:,:,p)
R = mmx('chol', CorrTensor, []); % File from here https://www.mathworks.com/matlabcentral/fileexchange/37515-mmx
e = randn(K,N,M);
e = pagemtimes(e, R),
Bruno Luong
Bruno Luong on 27 Oct 2021
Edited: 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
virginie marchionni
virginie marchionni on 27 Oct 2021
Thank you Christine for your comment, I will have a look at pagetimes. In my problem the NxN matrices are correlation matrices and I have M of them for M time istant, I will then have to generate M correlated brownian motions (e)
e = randn(K,N,M);
CorrTensor = Adj_Rho(Rho,dt,Param); % NxNxM
for m=1:M
t_chol = chol(squeeze(CorrTensor(:,:,m))); % NxN
e(:,:,m) = e(:,:,m)*t_chol;
end
e % KxNxM
What I'm looking for is if there is a way to avoid the for-loop.
Bruno Luong
Bruno Luong on 27 Oct 2021
Edited: 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

Find more on Linear Algebra in Help Center and File Exchange

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!