Performing matrix operation without loop

Hi,

I have two matrices (A and B), each of size, 4000x8000x3 (i,j,k).

The 'k th' dimension of both matrix contains median and 95% CI values.

To be robust by combining the 95%CI of both matrices and calculate the product of these matrices, I distribute each ' ith & jth' value of a matrix across the 'k th' values (95%CI) 1000 times for A and B. and obtain a 1000x1000 matrix for each ith and jth element. From which I calculate the mean and 95% CI for each ith and jth element.

However the below code takes > 4.5 hours to complete, I would request an easier and quicker way to do it. Thank you for the help

tic
for i=1:4000
  for j=1:8000;

      %%%distribute A lognormally
A_err= (log(A(i,j,3)) - log(A(i,j,1)))/1.96;
A_dist= lognrnd(log(A(i,j,1)), A_err,[1,1000]);
A_dist_m=repmat(A_dist,1000,1);


%%%%distribute B

B_err= (log(B(i,j,3)) - log(B(i,j,1)))/1.96;
B_dist= lognrnd(log(B(i,j,1)), B_err,[1,1000]);
% histogram (B_dist)
B_dist=B_dist';
B_dist_m=repmat(B_dist,1,1000);

C=B_dist_m.*A_dist_m;


range = prctile(C(:),[2.5 97.5]);
middle = mean(mean(C,1),2);

pr=[middle,range];
C1 (i,j,:)=C;
  end;

end
toc

[J. Alex Lee]

You can at least generate your A_err and B_err matrices outside the loop (to get 4000x8000 2D matrices)

A_err = (log(A(:,:,3)) - log(A(:,:,1)))/1.96

But i doubt that's the bottleneck.

Repmat is probably taking a while. Doesn't quite get you out of the loop paradigm, but I think this will give you same results instead of transposing and using repmat

A_dist= lognrnd(log(A(i,j,1)), A_err,[1,1000]);
B_dist= lognrnd(log(B(i,j,1)), B_err,[1000,1]); % transpose here so don't need to do it later

C = A_dist .* B_dist; % implicit expansion

[Matt J]

I find it peculiar that you are calculating ensemble estimates of mean and percentiles for data whose statistical distribution you already know? Is there an ulterior motive to this?

[SChow]

Hi Matt, I am trying to combine the uncertainties (95%CI) of A and B> Thought this would be a good way to progress, any suggestions are welcome

Answer 1

[Matt J]

Seems to me you could also just use the known quantile function for a log-normal distribution,

https://en.wikipedia.org/wiki/Log-normal_distribution

That wouldn't require any loops at all and you could compute any confidence interval that you want...

muA=log(A(:,:,1));
muB=log(B(:,:,1));
muC=muA+muB;

sigA=(log(A(:,:,3))-log(A(:,1)) )/1.96;
sigB=( log(B(:,:,3))-log(B(:,1)) )/1.96;
sigC=sigA+sigB;

Quant=@(p) exp( muC + sqrt(2).*erfinv(2*p-1).*sigC ); %quantile function

[SChow]

Thanks Matt!!

That is a great solution!!

Many thanks for such perseverance

Answer 2

[Matt J]

I think this should be faster. Obviously, it will be even faster if you can accept fewer than nSamples=1000 randomized samples for the calculations. Maybe 100 would be enough?

nSamples=1e3;
[M,N,~]=size(A);
e=ones(1,1,nSamples,'single');

res=@(z) reshape(single(z),M,10,[]);

muA=log(A(:,:,1));
muB=log(B(:,:,1));
muC=res(muA+muB);

sigA=(log(A(:,:,3))-log(A(:,1)) )/1.96;
sigB=( log(B(:,:,3))-log(B(:,1)) )/1.96;
sigC=res(sigA+sigB);

[middle,rangeLOW,rangeHIGH]=deal(nan(M,10));


tic
for n=1:size(muC,3)
    
    
    C_dist=lognrnd(muC(:,:,n).*e, sigC(:,:,n).*e);
    
    middle(:,:,n)=mean(C_dist,3);
    
    range = prctile(C_dist,[2.5 97.5],3);
    
    rangeLOW(:,:,n)=range(:,:,1);
    rangeHIGH(:,:,n)=range(:,:,2);
    
end
toc

C1=cat(3,middle,rangeLOW,rangeHIGH);
C1=reshape(C1,M,[],3);

[SChow]

Thanks for the answer Matt. Limiting to ~100 samples would not be great, but its better than nothing

However this does not work, it provides Inf values.

I tried to apply your code with two arbitary matrices

A=cat(3, 1.12.*ones(4000,8000), 1.02.*ones(4000,8000), 1.17.*ones(4000,8000));
B=cat(3, 400.*ones(4000,8000), 200.*ones(4000,8000), 800.*ones(4000,8000));

[Matt J]

Make sure you're using my latest version. I've been tweaking and editing it gradually. With nSamples=1000, I'm getting an expected completion time on my machine of about 17 minutes.

[SChow]

Thanks Matt, I redid with your latest version.

Sorry but the resultant output is of size 4000x10x3. It should ideally be required to be 4000x8000. I tried without your reshape function, but Matlab shows out of memory error

[Matt J]

Make sure the for-loop is running over the full range

for n=1:size(muC,3)