How to compute cholesky to all slice of a tensor?

Question

virginie marchionni 2021-10-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1571563-how-to-compute-cholesky-to-all-slice-of-a-tensor

编辑： Bruno Luong 2021-10-27

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Bruno Luong 2021-10-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1571563-how-to-compute-cholesky-to-all-slice-of-a-tensor#answer_816318

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 个评论
显示无隐藏无

Bruno Luong 2021-10-27

编辑：Bruno Luong 2021-10-27

在 MATLAB Online 中打开

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 2021-10-27

编辑：Bruno Luong 2021-10-27

在 MATLAB Online 中打开

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 2021-10-26

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1571563-how-to-compute-cholesky-to-all-slice-of-a-tensor#answer_816853

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 个评论
显示无隐藏无

virginie marchionni 2021-10-27

在 MATLAB Online 中打开

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 2021-10-27

编辑：Bruno Luong 2021-10-27

在 MATLAB Online 中打开

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