Multidimensional arrays of 3-D objects: how to vectorise inner products

Question

In this question, I discussed two custom functions to multiply arrays of 3x3 matrices and 3x1 vectors, preserving the structure of the 3-dimensional (matrix) inner product and making the whole process as computationally efficient and fast as possible.

I have now generalised those functions to multidimensional arrays (NxN) of 3x4 matrices and 3x1 vectors. Here are the functions I have written, which make use of for loops.

BlockScalar

This function should multiply the ij element (a scalar) of the (NxN matrix) nv by the ij element (3x3 matrix) of A (NxNx3x3 matrix). So it's essentially a multidimensional version of the product of a matrix by a scalar.

function [B] = BlockScalar(nv,A)

    N=size(nv,1);
    B=zeros(N,N,3,3);

    for i=1:N
        for j=1:N 
            B(i,j,:,:)= nv(i,j).*A(i,j,:,:);
        end
    end

end

--------BlockScalar Example

Inputs:

N=2;
A = shiftdim( repmat( eye(3,3), 1, 1, N, N ), 2 );
nv=[1 2; 3 4];

Output:

BlockScalar(nv,A)

ans(:,:,1,1) =

 1     2
 3     4

ans(:,:,2,1) =

 0     0
 0     0

ans(:,:,3,1) =

 0     0
 0     0

ans(:,:,1,2) =

 0     0
 0     0

ans(:,:,2,2) =

 1     2
 3     4

ans(:,:,3,2) =

 0     0
 0     0

ans(:,:,1,3) =

 0     0
 0     0

ans(:,:,2,3) =

 0     0
 0     0

ans(:,:,3,3) =

 1     2
 3     4

BlockMatrix

This second function does not work at the moment because I am struggling to implement the matrix product A*u between the ij-th element (which is a 3x3 matrix) of A and a column vector that contains the 3 components of the ijth element of u. As you may easily see, I would like this to be a multidimensional generalisation of the matrix*vector product in 3-D.

 function [B] = BlockMatrix(A,u)

    N = size(u,2);
    B = zeros(N,N,3);

  for i=1:N
     for j=1:N
             B(i,j,:)= reshape(reshape(A(i,j,:,:),[3,3])*reshape(u(i,j,:),[1 3]),size(u));
     end

 end

-------BlockMatrix Example

If the input is a generalised identity matrix (NxN elements each of which is a 3x3 identity matrix), and an NxN matrix made of 3x1 vectors:

 N=2;   

 A = 4.*shiftdim( repmat( eye(3,3), 1, 1, N, N ), 2 );

 c = ones(2,2);
 V(1,1,:)=[1 2 3];
 u = c.*V;

The desired output is clearly an object which has the structure of V (NxN matrix made of 3x1 vectors) where each of the elements is the matrix product of reshape(A(i,j,:,:),[3 3]) and reshape(V(i,j,:),[1 3]). That is:

i=1;j=1;
reshape(B(i,j,:),[3,1])


ans =

 4
 8
 12

for any i and j.

Full output, for completeness:

B(:,:,1) =

 4     4
 4     4

B(:,:,2) =

 8     8
 8     8

B(:,:,3) =

 12     12
 12     12

Questions

I struggle to (0) make BlockMatrix work; (1) find a way to properly vectorise this, and (2) I am not even particularly sure the vectorised version would be faster.

Any help in answering the points above will be much appreciated.

Answer 1

For the first function:

B = bsxfun(@times, A, nv);

For the second:

B = sum(bsxfun(@times, A, reshape(u, [size(u,1) size(u,2) 1 size(u,3)])), 4);