3d Matrix - Extract Subarray and Multiply by Conjugate Transpose without Forloops

I have a 3d matrix A(i,j,k) of size [1:100,1:10000,1:989]. On the kth index I want to extract the 989 elements into a vector u and form the product u*ctranspose(u), for each of the indices.

Using a double for loop, which is not something one should do in Matlab,

% this is an evil double for loop - I want to avoid doing this.
for ii=1:100
    for jj=1:10000
        u=squeeze(A(ii,jj,:));
        % somehow compute and store uu in a vectorized way, but how?
        uu=u*ctranspose(u); % note that uu is a 989x989 matrix, not a vector.
    end %ii
end %jj

This would be really slow. Is there a vectorized way to the above, so I am not doing a double for loop?

[Matt J]

Yes, but the result will consume 7 TB in double floats

989^2*10000*100*8/2^40

ans = 7.1168

Outer-products are a bad idea for many other reasons, too. What are you going to do with this even if you could store it?

[Science Machine]

Let's suppose we slice off the i (do 1 element at a time instead of 100). Then the variable is 70 GB.

[Matt J]

70 GB still seems like a lot. Do you have that much RAM?

[Science Machine]

yes, i have locally 196 gb ram, and on the cluster I have about 1 tb ram.

[Matt J]

Well, outer products are still a very inefficient thing to do in most circumstances. It would be advisable for you to elaborate on why you think you need this, and how you would use it in subsequent steps.

[Bruno Luong]

What a waste of memory (and cpu) to explicitly expand rank-1 matrices.

[Science Machine]

@Bruno Luong I was confused to whom your comment was addressed.

[Bruno Luong]

Of course yours

[Science Machine]

It would be cheaper to use correlation (xcorr) but I would prefer to do things exactly to the letter until I have a full working procedure. This is only a step of the procedure and I don't want to introduce any uncertainty, however small.

[Matt J]

@Science Machine The connection between xcorr and what you're attempting is not as apparent as you seem to think. A rank-1 matrix doesn't represent a correlation operation.

[Science Machine]

It was not clear what @Bruno Luong meant - does he think it's okay to modify a single step of a procedure? And what is his alternative. I also learned in my matrix methods class that outer product is rank-1, but what he wants to do with that information.

@Matt J thanks for your help. that seems to work. Nevermind about correlations.

[Matt J]

It was not clear what @Bruno Luong meant - does he think it's okay to modify a single step of a procedure? And what is his alternative.

That depends on what you plan to do with uu. Note for example that if you have column vectors u and x, then this,

y=u*(u'*x);

produces the same result as this,

y=(u*u')*x;

but the former is much more efficient, both in terms of memory consumption and CPU time.

[Science Machine]

I am finding a spectral density matrix, I am not solving a linear equation!

[Matt J]

My comment wasn't related to equation-solving. It is relevant to whether you will be doing matrix multiplications with uu.

[Science Machine]

Interesting. I see now. So if I have $u(t,r) u^*(t',r)r$ this should be a way to avoid the outer product still? by doing $u(t,r)(u^*(t',r)r)$

[Matt J]

Probably, but you haven't explained what t and r are.

[Science Machine]

t and ras parameters are the ii and jj indices from above u of size 1:100 and 1:10000, and the vector $$ has parameters jj also.

Answer 1

[Matt J]

[m,n,p]=size(A);

Ar=reshape(A,[],p).';

UU = reshape(Ar,p,1,[]).*conj(reshape(Ar,1,p,[]));

UU=reshape(UU,p,p,m,n); %obtain the (i,j)-th outer product as UU(:,:,i,j)

[Science Machine]

what is out?

UU=reshape(out,p,p,m,n);
Unrecognized function or variable 'out'. 

do you mean,

out=reshape(UU,p,p,m,n);  

?

[Matt J]

A typo. I fixed it.

[Science Machine]

How did you learn to do such things? I knew reshape function but it did not occur to me to use it here.