# How to compute cholesky to all slice of a tensor?

2 views (last 30 days)
virginie marchionni on 25 Oct 2021
Edited: Bruno Luong on 27 Oct 2021
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

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 CommentsShow NoneHide None
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 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.

### More Answers (1)

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 CommentsShow NoneHide None
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 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

### Categories

Find more on Linear Algebra in Help Center and File Exchange

R2020a

### Community Treasure Hunt

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

Start Hunting!