How to compute cholesky to all slice of a tensor?

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

Answer 1

[Bruno Luong]

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.

[Bruno Luong]

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]

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.

Answer 2

[Christine Tobler]

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?

[virginie marchionni]

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]

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