Thread Subject: vectorise a reshape loop

Hi, I am having trouble vectorising the below for loop. Can
anyone help me please? I am expecting M to consist of 3 rows
of 2x2 matrices each with 16 deep. Thanks in advance.

for p=1:3
q=p;
M=reshape([a(:,q) c(:,q) b(:,q) d(:,q)].', 2, 2, []);
end

where a,b,c,d are 3x16 matrice such as:
(imaginary parts can be ignored)
   1.0000 + 0.0000i 0.9999 + 0.0000i 0.9999 + 0.0001i
   1.0000 + 0.0000i 0.9999 + 0.0000i 0.9999 + 0.0001i
   1.0000 + 0.0000i 0.9999 + 0.0000i 0.9999 + 0.0001i
   1.0000 + 0.0000i 0.9999 + 0.0000i 0.9999 + 0.0001i
   1.0004 + 0.0000i 1.0003 + 0.0000i 1.0002 + 0.0001i
   1.0004 + 0.0000i 1.0003 + 0.0000i 1.0002 + 0.0001i
   1.0003 + 0.0000i 1.0003 + 0.0000i 1.0002 + 0.0001i
   1.0003 + 0.0000i 1.0003 + 0.0000i 1.0002 + 0.0001i
   1.0051 + 0.0000i 1.0050 + 0.0000i 1.0050 + 0.0001i
   1.0068 + 0.0000i 1.0068 + 0.0000i 1.0067 + 0.0001i
   1.0111 + 0.0000i 1.0110 + 0.0000i 1.0110 + 0.0001i
   1.0178 + 0.0000i 1.0178 + 0.0000i 1.0177 + 0.0001i
   1.0500 + 0.0000i 1.0499 + 0.0000i 1.0498 + 0.0000i
   1.1666 + 0.0000i 1.1666 + 0.0000i 1.1665 + 0.0000i
   1.2666 + 0.0000i 1.2666 + 0.0000i 1.2665 + 0.0000i
   3.2499 + 0.0000i 3.2497 + 0.0000i 3.2493 + 0.0000i

"Dave Brackett" <davebrackett@hotmail.com> wrote in message
<g5krk2$77v$1@fred.mathworks.com>...
> Hi, I am having trouble vectorising the below for loop. Can
> anyone help me please? I am expecting M to consist of 3 rows
> of 2x2 matrices each with 16 deep. Thanks in advance.

Try that again? What is the shape of M?
Is it 3x2x2x16? This is inconsistent with
what you have shown below.


> for p=1:3
> q=p;
> M=reshape([a(:,q) c(:,q) b(:,q) d(:,q)].', 2, 2, []);
> end
>
> where a,b,c,d are 3x16 matrice such as:
> (imaginary parts can be ignored)
> 1.0000 + 0.0000i 0.9999 + 0.0000i 0.9999 + 0.0001i

Um, that is a 16x3 matrix.

Anyway, what loop?

Since you overwrite M on each pass
through, that non-loop does nothing
but define M at the last value of p.

We can't help you unless you will think
clearly enough about what you want to
explain it.

You might start by preallocating M to be
the proper size. Use zeros.

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in
message <g5kthu$svb$1@fred.mathworks.com>...
> "Dave Brackett" <davebrackett@hotmail.com> wrote in message
> <g5krk2$77v$1@fred.mathworks.com>...
> > Hi, I am having trouble vectorising the below for loop. Can
> > anyone help me please? I am expecting M to consist of 3 rows
> > of 2x2 matrices each with 16 deep. Thanks in advance.
>
> Try that again? What is the shape of M?
> Is it 3x2x2x16? This is inconsistent with
> what you have shown below.
>
>
> > for p=1:3
> > q=p;
> > M=reshape([a(:,q) c(:,q) b(:,q) d(:,q)].', 2, 2, []);
> > end
> >
> > where a,b,c,d are 3x16 matrice such as:
> > (imaginary parts can be ignored)
> > 1.0000 + 0.0000i 0.9999 + 0.0000i 0.9999 + 0.0001i
>
> Um, that is a 16x3 matrix.
>
> Anyway, what loop?
>
> Since you overwrite M on each pass
> through, that non-loop does nothing
> but define M at the last value of p.
>
> We can't help you unless you will think
> clearly enough about what you want to
> explain it.
>
> You might start by preallocating M to be
> the proper size. Use zeros.
>
> John
>
>

Ok, sorry for any confusion, I will try again.

I want a 3D matrix, consisting of in this case, 3 2x2
matrices (as shown below) plus 15 additional 2x2 matrices
stacked behind them, giving 3 rows of 16 2x2 matrices all
contained within one matrix.

[a1 b1 [a2 b2 [a3 b3
 c1 d1] c2 d2] c3 d3]

The amended 'for' loop shown below intends to do this. I
have changed it so M should store each 2x2 matrix, but it is
isn't right as I'm not sure how to do it right.

for p=1:3
q=p;
r=1:16;
M(:,q,r)=reshape([a(:,q) c(:,q) b(:,q) d(:,q)].', 2, 2, []);
end

As you rightly point out a,b,c, and d are indeed 16x3
matrices. Preallocation of M can come later.

Hopefully that explanation is a bit clearer. Thanks.

"John D'Errico" <woodchips@rochester.rr.com> wrote in
message <g5kthu$svb$1@fred.mathworks.com>...
> "Dave Brackett" <davebrackett@hotmail.com> wrote in message
> <g5krk2$77v$1@fred.mathworks.com>...
> > Hi, I am having trouble vectorising the below for loop. Can
> > anyone help me please? I am expecting M to consist of 3 rows
> > of 2x2 matrices each with 16 deep. Thanks in advance.
>
> Try that again? What is the shape of M?
> Is it 3x2x2x16? This is inconsistent with
> what you have shown below.
>
>
> > for p=1:3
> > q=p;
> > M=reshape([a(:,q) c(:,q) b(:,q) d(:,q)].', 2, 2, []);
> > end
> >
> > where a,b,c,d are 3x16 matrice such as:
> > (imaginary parts can be ignored)
> > 1.0000 + 0.0000i 0.9999 + 0.0000i 0.9999 + 0.0001i
>
> Um, that is a 16x3 matrix.
>
> Anyway, what loop?
>
> Since you overwrite M on each pass
> through, that non-loop does nothing
> but define M at the last value of p.
>
> We can't help you unless you will think
> clearly enough about what you want to
> explain it.
>
> You might start by preallocating M to be
> the proper size. Use zeros.
>
> John
>
>

Ok, sorry for any confusion, I will try again.

I want a 3D matrix, consisting of in this case, 3 2x2
matrices (as shown below) plus 15 additional 2x2 matrices
stacked behind them, giving 3 rows of 16 2x2 matrices all
contained within one matrix.

[a1 b1 [a2 b2 [a3 b3
 c1 d1] c2 d2] c3 d3]

The amended 'for' loop shown below intends to do this. I
have changed it so M should store each 2x2 matrix, but it is
isn't right as I'm not sure how to do it right.

for p=1:3
q=p;
r=1:16;
M(:,q,r)=reshape([a(:,q) c(:,q) b(:,q) d(:,q)].', 2, 2, []);
end

As you rightly point out a,b,c, and d are indeed 16x3
matrices. Preallocation of M can come later.

Hopefully that explanation is a bit clearer. Thanks.

cut out quotes from previous posts as they are making the
post unreadable, also sorry for double post of previous
reply - don't know what happened there!

ok, i have spent some more time and the below for loop now
shows what I am trying to achieve. However, the cells are
causing me some trouble as i don't get 2x2 matrices for each
prodH entry. I have read that cells don't allow
vectorisation so is there a better way that someone can
suggest? Thanks.

M=cell(1,3,16);
prodH=cell(1,3);

for p=1:3
    q=p;
    for r=1:16;
        M(1,q,r)={reshape([ai_col(r,q) ci_col(r,q)
bi_col(r,q) di_col(r,q)].', 2, 2, [])}
    end
    prodH(1,q)={[ndfun('mprod', M(1,q,:))]}
end

"Dave Brackett" <davebrackett@hotmail.com> wrote in message
<g5n4kd$jtp$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in
> message <g5kthu$svb$1@fred.mathworks.com>...
> > "Dave Brackett" <davebrackett@hotmail.com> wrote in message
> > <g5krk2$77v$1@fred.mathworks.com>...
> > > Hi, I am having trouble vectorising the below for loop. Can
> > > anyone help me please? I am expecting M to consist of 3 rows
> > > of 2x2 matrices each with 16 deep. Thanks in advance.
> >
> > Try that again? What is the shape of M?
> > Is it 3x2x2x16? This is inconsistent with
> > what you have shown below.
> >
> >
> > > for p=1:3
> > > q=p;
> > > M=reshape([a(:,q) c(:,q) b(:,q) d(:,q)].', 2, 2, []);
> > > end
> > >
> > > where a,b,c,d are 3x16 matrice such as:
> > > (imaginary parts can be ignored)
> > > 1.0000 + 0.0000i 0.9999 + 0.0000i 0.9999 + 0.0001i
> >
> > Um, that is a 16x3 matrix.
> >
> > Anyway, what loop?
> >
> > Since you overwrite M on each pass
> > through, that non-loop does nothing
> > but define M at the last value of p.
> >
> > We can't help you unless you will think
> > clearly enough about what you want to
> > explain it.
> >
> > You might start by preallocating M to be
> > the proper size. Use zeros.
> >
> > John
> >
> >
>
> Ok, sorry for any confusion, I will try again.
>
> I want a 3D matrix, consisting of in this case, 3 2x2
> matrices (as shown below) plus 15 additional 2x2 matrices
> stacked behind them, giving 3 rows of 16 2x2 matrices all
> contained within one matrix.

Do you want a 4-d matrix, that is 2x2x3x16?

Or does this totally unclear explanation
request a 2x2x48 matrix? What do you
mean when you say 3 rows 2x2 matrices?
Does this ask for a 3x2x2 array? Or for a
2x2x3 array?

For some reason, you do not accept what
it is you are trying to do, nor have you
taken the time or effort to explain it clearly.
Think clearly. If you want a single matrix
that contains 3 sets of 2x2 matrices, 16
times, then you have a 3x2x2x16 matrix.
A 4-d matrix.

Is this what you want?


> [a1 b1 [a2 b2 [a3 b3
> c1 d1] c2 d2] c3 d3]
>
> The amended 'for' loop shown below intends to do this. I
> have changed it so M should store each 2x2 matrix, but it is
> isn't right as I'm not sure how to do it right.
>
> for p=1:3
> q=p;
> r=1:16;
> M(:,q,r)=reshape([a(:,q) c(:,q) b(:,q) d(:,q)].', 2, 2, []);
> end

Again, this loop merely (tries to) overwrite
M on each pass through. Is this what you
want? I don't think so from your description.

 
> As you rightly point out a,b,c, and d are indeed 16x3
> matrices. Preallocation of M can come later.

No. Preallocation should come FIRST, for
several good reasons. Not the least of
which is your acceptance of what you are
trying to do.

 
> Hopefully that explanation is a bit clearer. Thanks.

Clear now. As mud.

What is your goal? If it is as I have suggested
it might be, then you can do everything in
one line. I'll do it in two lines to make it
more understandable.

M = reshape([a(:),b(:),c(:),d(:)],16,3,2,2);
M = permute(M,[3 4 2 1]);

This will create a single 2x2x3x16 array.

So far, I still have no idea if this is what
you really want. So slow down. Think
clearly about what you are doing. Then
explain it clearly. Communication is a
valuable skill. Learn it.

John

I think the communication issue stems from the fact i didn't
understand the number of array dimensions could be above 3.
To me, a matrix with rows, columns, and pages was always a
3D matrix (as it forms a 3D 'block'). Evidently this is not
the case.

The reshape and permute functions you suggest do achieve
what I wanted thanks:
M=reshape([ai_col(:),bi_col(:),ci_col(:),di_col(:)],16,3,2,2);
M=permute(M,[3 4 2 1]);

For each of the 3 sets of 16 2x2 arrays contained within M I
want to find the product in the dimension of the 16 2x2
matrices. I have been using your ndfun mex file as follows:

prodH=zeros(2,2,3);
for p=1:3
    q=p;
    prodH(1,q)=ndfun('mprod',M(:,:,p,:))
end

However, i get the error that 'assignment has more
non-singleton rhs dimensions than non-singleton subscripts'.
Could you show me how to do this correctly pls? Thanks for
your help.

"Dave Brackett" <davebrackett@hotmail.com> wrote in message
<g5ngf2$fdp$1@fred.mathworks.com>...
> I think the communication issue stems from the fact i didn't
> understand the number of array dimensions could be above 3.
> To me, a matrix with rows, columns, and pages was always a
> 3D matrix (as it forms a 3D 'block'). Evidently this is not
> the case.

Ah. This was the problem. Yes, I don't think
there is any real limit on the number of
dimensions. Yes, there may be a limit, but
it is probably some large number.

 
> The reshape and permute functions you suggest do achieve
> what I wanted thanks:
> M=reshape([ai_col(:),bi_col(:),ci_col(:),di_col(:)],16,3,2,2);
> M=permute(M,[3 4 2 1]);
>
> For each of the 3 sets of 16 2x2 arrays contained within M I
> want to find the product in the dimension of the 16 2x2
> matrices. I have been using your ndfun mex file as follows:

Its not my ndfun.


> prodH=zeros(2,2,3);
> for p=1:3
> q=p;
> prodH(1,q)=ndfun('mprod',M(:,:,p,:))
> end
>
> However, i get the error that 'assignment has more
> non-singleton rhs dimensions than non-singleton subscripts'.
> Could you show me how to do this correctly pls? Thanks for
> your help.

What would be the result of this product?
Perhaps just a 2x2 array? How do you
expect a 2x2 array to fit into a single
array element of prodH? Why would you
index the result as prodH(1,q) ?

If that result is a 2x2 array, and you have
preallocated prodH as 2x2x3, then should
you not do it as

   prodH(:,:,q)=ndfun('mprod',M(:,:,p,:));

John

> Its not my ndfun.

sorry, not sure why i suggested it was yours. it is from
http://www.mit.edu/~pwb/matlab/ndfun/ndfun.m

 
> What would be the result of this product?
> Perhaps just a 2x2 array? How do you
> expect a 2x2 array to fit into a single
> array element of prodH? Why would you
> index the result as prodH(1,q) ?
>
> If that result is a 2x2 array, and you have
> preallocated prodH as 2x2x3, then should
> you not do it as
>
> prodH(:,:,q)=ndfun('mprod',M(:,:,p,:));

after i had posted that i realised it was wrong and my line
in the for loop now looks like:
prodH(:,:,q)=ndfun('mprod',M(:,:,q,:))

however, i still get the same error. i am expecting 3 2x2
arrays, 1 for each product of the 16 arrays.

"Dave Brackett" <davebrackett@hotmail.com> writes:

> I think the communication issue stems from the fact i didn't
> understand the number of array dimensions could be above 3.
> To me, a matrix with rows, columns, and pages was always a
> 3D matrix (as it forms a 3D 'block'). Evidently this is not
> the case.
>
> The reshape and permute functions you suggest do achieve
> what I wanted thanks:
> M=reshape([ai_col(:),bi_col(:),ci_col(:),di_col(:)],16,3,2,2);
> M=permute(M,[3 4 2 1]);
>
> For each of the 3 sets of 16 2x2 arrays contained within M I
> want to find the product in the dimension of the 16 2x2
> matrices. I have been using your ndfun mex file as follows:
>
> prodH=zeros(2,2,3);
> for p=1:3
> q=p;
> prodH(1,q)=ndfun('mprod',M(:,:,p,:))
> end
>
> However, i get the error that 'assignment has more
> non-singleton rhs dimensions than non-singleton subscripts'.
> Could you show me how to do this correctly pls? Thanks for
> your help.

prodH is 2x2x3. You index it with only two indices. What do you want
the third index to be?

For each thing you are trying to assign (in your case, only one!),
figure out the dimensions of the right-hand side. Try it first
mentally, and if you're not sure, assign it to a variable and ask MATLAB
to tell you what the size is.

Then, look at where you want to assign it. Your assignment indices
must provide a place in the output which exactly matches the dimensions
of the right side.

For instance, I can't say:

A(1, 2:4) = B(1:4, 2);

Why?
left side is 1x3, right side is 4x1

I can't say:
A(1,1,2:4) = B(1,2:4)

Why?
left side is 1x1x3, right side is 1x3

-Peter

In article <g5nj77$j6j$1@fred.mathworks.com>,
John D'Errico <woodchips@rochester.rr.com> wrote:
>"Dave Brackett" <davebrackett@hotmail.com> wrote in message
><g5ngf2$fdp$1@fred.mathworks.com>...
>> I think the communication issue stems from the fact i didn't
>> understand the number of array dimensions could be above 3.
>> To me, a matrix with rows, columns, and pages was always a
>> 3D matrix (as it forms a 3D 'block'). Evidently this is not
>> the case.

>Ah. This was the problem. Yes, I don't think
>there is any real limit on the number of
>dimensions. Yes, there may be a limit, but
>it is probably some large number.

The number of dimensions is -at least- 2^16+1

foo = [5;7];
bar = permute(foo,[2:65537,1]);

I have it working on 2^20+1 but it is taking a fair bit of time.
--
  "The slogans of an inadequate criticism peddle ideas to fashion"
                                              -- Walter Benjamin

> figure out the dimensions of the right-hand side. Try it
first
> mentally, and if you're not sure, assign it to a variable
and ask MATLAB
> to tell you what the size is.
>
> Then, look at where you want to assign it. Your
assignment indices
> must provide a place in the output which exactly matches
the dimensions
> of the right side.
>
> For instance, I can't say:
>
> A(1, 2:4) = B(1:4, 2);
>
> Why?
> left side is 1x3, right side is 4x1
>
> I can't say:
> A(1,1,2:4) = B(1,2:4)
>
> Why?
> left side is 1x1x3, right side is 1x3
>
> -Peter


ok, the size of my right side is 2x2x90x16 and the required
size for my left is 2x2x90x1 so I have written this as:

prodH=zeros(2,2,length(freq),1);
for p=1:length(freq)
    q=p;
    prodH(:,:,q,1)=ndfun('mprod',M(:,:,q,:))
end

I still get:

??? Assignment has more non-singleton rhs dimensions than
non-singleton
subscripts

Error in ==> new_for_gatool_2 at 110
        prodH(:,:,q,1)=ndfun('mprod',M(:,:,q,:))

What am i still doing wrong? Thanks for your help.

In article <g61k6f$1gj$1@fred.mathworks.com>,
Dave Brackett <davebrackett@hotmail.com> wrote:

> prodH(:,:,q,1)=ndfun('mprod',M(:,:,q,:))

>I still get:

>??? Assignment has more non-singleton rhs dimensions than
>non-singleton
>subscripts

try squeeze(M(:,:,q,:))

M(:,:,q,:) is a 4D array whose third dimension happens to be 1.
Matlab only automatically collapses trailing singleton array
dimensions.
--
  "Nothing recedes like success." -- Walter Winchell

"Dave Brackett" <davebrackett@hotmail.com> writes:

>> figure out the dimensions of the right-hand side. Try it first
>> mentally, and if you're not sure, assign it to a variable and ask
>> MATLAB to tell you what the size is.
>>
>> Then, look at where you want to assign it. Your assignment indices
>> must provide a place in the output which exactly matches the
>> dimensions of the right side.
>>
>> For instance, I can't say:
>>
>> A(1, 2:4) = B(1:4, 2);
>>
>> Why?
>> left side is 1x3, right side is 4x1
>>
>> I can't say:
>> A(1,1,2:4) = B(1,2:4)
>>
>> Why?
>> left side is 1x1x3, right side is 1x3

Pay VERY careful to the examples above, again. In the second example, I
am not asking for the size of B, but for the size of entire right-hand
side. B is at least 1x4, but the indexing makes it 1x3.

> ok, the size of my right side is 2x2x90x16 and the required
> size for my left is 2x2x90x1 so I have written this as:
>
> prodH=zeros(2,2,length(freq),1);
> for p=1:length(freq)
> q=p;
> prodH(:,:,q,1)=ndfun('mprod',M(:,:,q,:))
> end

The size of the left side IS NOT 2x2x90x1. q is a scalar, so the left
side is 2x2x1x1. The right side, I assume, is 2x2x1x16. But don't guess
or believe me. If you can't figure out your dimension issues, I meant
what I said about "if you're not sure, assign it to a variable and ask
MATLAB to tell you what the size is".

tmp = ndfun('mprod', M(:,:,q,:));
size(tmp)

I could take a guess at fixing it, but only you know what you want
where. Actually, I think ndfun('mprod') is not doing what you think
it's is doing. Do what I suggest above, and I think you'll see what I
mean. Then go back and read the help for ndfun, carefully.


-Peter

> The size of the left side IS NOT 2x2x90x1. q is a scalar,
> so the left side is 2x2x1x1. The right side, I assume, is
> 2x2x1x16. But don't guess or believe me. If you can't
> figure out your dimension issues, I meant what I said
> about "if you're not sure, assign it to a variable and
> ask MATLAB to tell you what the size is".
> tmp = ndfun('mprod', M(:,:,q,:));
> size(tmp)

I am getting quite confused with what I thought would be a
simple thing. prodH and M are contained within a loop of 90
iterations, so the size of the resulting prodH matrix after
running the loop should be 2x2x90x1.

M is of the size 2x2x90x16 and i am iterating through this
using a loop to find the product of the 4th dimension
(length 16) for each of the indices in the 3rd dimension
(length 90).

This would then result in prodH of the size 2x2x90x1 with
each 2x2 matrix being the product of the 16 in M.

I hope that makes sense. That is what I am trying to achieve
and thought that ndfun('mprod') would be enable that. If my
loop is confusing and seems to suggest something different
then apologies.

Please could you show me how to write this correctly as I am
starting to go mad with it? Thanks.

"Dave Brackett" <davebrackett@hotmail.com> writes:

>> The size of the left side IS NOT 2x2x90x1. q is a scalar,
>> so the left side is 2x2x1x1. The right side, I assume, is
>> 2x2x1x16. But don't guess or believe me. If you can't
>> figure out your dimension issues, I meant what I said
>> about "if you're not sure, assign it to a variable and
>> ask MATLAB to tell you what the size is".
>> tmp = ndfun('mprod', M(:,:,q,:));
>> size(tmp)
>
> I am getting quite confused with what I thought would be a
> simple thing. prodH and M are contained within a loop of 90
> iterations, so the size of the resulting prodH matrix after
> running the loop should be 2x2x90x1.
>
> M is of the size 2x2x90x16 and i am iterating through this
> using a loop to find the product of the 4th dimension
> (length 16) for each of the indices in the 3rd dimension
> (length 90).
>
> This would then result in prodH of the size 2x2x90x1 with
> each 2x2 matrix being the product of the 16 in M.
>
> I hope that makes sense. That is what I am trying to achieve
> and thought that ndfun('mprod') would be enable that. If my
> loop is confusing and seems to suggest something different
> then apologies.
>
> Please could you show me how to write this correctly as I am
> starting to go mad with it? Thanks.

Until now, I have been unable to understand exactly what you want where.
With this explanation, I think I can help.

I admit that I am still frustrated that you refuse to follow my
debugging advice by looking at the size of the matrix generated by
ndfun, and then reading the help for ndfun. So instead of just writing
the answer, I'll go through the exercise here, hopefully providing an
illustration of the debugging process I was asking you to try.

M = rand(2,2,90,16);
tmp = ndfun('mprod', M(:,:,1,:));
size(tmp)

ans =
   2 2 16


Hmm, why is it 2x2x16?

"help ndfun"
    'mprod' behaves differently. It cumulatively multiplies a set
    of matrices and produces a single 2D output. The equivalent
    code is:
        C = A(:,:,1);
        for i=2:N
            C = C * A(:,:,i);
        end
    2D inputs return themselves. Inputs with more than 3
    dimensions collapse the 3rd dimension only. So with an A of
    size [2 2 7 3 4],
        C = ndfun('mprod', A);
    is equivalent to
        for i=1:3
            for j=1:4
                C(:,:,i,j)=ndfun('mprod',M(:,:,:,i,j));
            end
        end
    and C will have size [2 2 3 4].


Reading that carefully, especially the second example, suggests that
inputting a 2x2x1x16 matrix will produce a 2x2x16 output, essentially
doing nothing. Since ndfun('mprod') works down the third dimension, you
want your x16 size to be on that third dimension. "Squeeze" will do
that in this case, but it's just a special case of permute, so whichever
you want.

So this should work:

prodH=zeros(2,2,length(freq),1);
for p=1:length(freq)
    q=p;
    prodH(:,:,q,1)=ndfun('mprod', squeeze(M(:,:,q,:)))
end

Let's check:

Right side:

M(:,:,q,:) is 2x2x1x16
squeeze(M(:,:,q,:)) is 2x2x16
   (you can always ignore trailing 1 dimensions)
output of ndfun is 2x2

Left side:
prodH(:,:,q,1) is 2x2

Good!

OK, one more thing. Again looking carefully at the help for ndfun, the
example shows a case where the mprod is repeated for higher-dimension
indices. So it appears that if you put your x90 dimension at the end
(in the fourth dimension), ndfun will automatically do what you want,
which is to collapse the 3rd dimension (16) by multiplying the 2x2
matrices, then do this for all of the 4d "pages". So you want to switch
3rd and 4th dimensions of M:

This should do everything, sans loop:

prodH = ndfun('mprod', permute(M, [1 2 4 3]));


-Peter

Thank you for your clear explanation and help with this. I
appreciate your time and effort.