Full Matrix Multiplication Vectorization

Hi,

I'm interested in vectorizing one of my for loops to improve the efficiency. I don't want to perform a element-wise matrix multiplication for this, I want to perform a full matrix multiplication. Can anyone help? A sample of my code is similar the the one shown below

    for x=i:-1:1                            
        if x==j
            V(:,:,x)=dj;                   
        else
            V(:,:,x);                       
        end
        if x==i                                                         
            U=V(:,:,x);
        else
            U=V(:,:,x)*U;                   
        end
    end

Ok to expand on my initial question and provide a little more background. i and j both range from 1 to 6, V is a 4 x 4 x 6 matrix, dj is the derivative of the matrix V(:,:,j) (I realize this can be coded better). I then want to multiply all the layers of the matrix V together (which is the for loop), from i to 1 but replace the layer of the matrix equal to j with it's derivative, dj (which is computed by the x==j if statement). Hopefully that makes a bit more sense - Thanks

if i < j
    U=zeros(4,4);                           
else
    dj= [0 -1 0 0;...                       
        1 0 0 0;...
        0 0 0 0;...
        0 0 0 0] * V(:,:,j);
      for x=i:-1:1                           
          if x==j
              V(:,:,x)=dj;                 
          else
              V(:,:,x);                      
          end
          if x==i
              U=V(:,:,x);
          else
              U=V(:,:,x)*U;                   
          end
      end
 end

[Geoff Hayes]

Will - you will need to provide more code than the above in order to better explain the problem that you wish to improve upon. In your for loop, you are iterating from -1 to 1 using x as the indexing variable. Then your code uses x to index the third dimension of V which will fail when x is -1 or 0. Please describe what you are trying to do here.

As well, the first else statement doesn't perform any assignment. What should be happening to V(:,:,x) in this case?

[Will]

In answer to your specific question Geoff, the index of x will never fail as it will run from a maximum of 6 to a minimum of 1. The aim of this for loop is to multiply the layers of the matrix at index x together replacing the layer x when it is equal to j with dj.

In the final 'else' statement that you refer to, you are correct in saying that there is no assignment. I believe that including an else statement within an 'if' statement is classed as good coding practice. V(:,:,x) should not change if the 'if' statement is false.

Answer 1

[Geoff Hayes]

Will - if you have an else statement that does nothing, then there is no need to include it. With that in mind, you could reduce your code to something simpler.

 if i < j
    U=zeros(4,4);                           
 else
    dj= [0 -1 0 0;...                       
         1 0 0 0;...
         0 0 0 0;...
         0 0 0 0] * V(:,:,j);

    % do stuff
 end

Now to elaborate on the do stuff part of the code: look at what we know. We enter the else if j is less than or equal to i. As both i and j run from 1 through to six, then we can replace the

 if x==j
     V(:,:,x)=dj;                 
 else
     V(:,:,x);                      
 end

with

 V(:,:,j) = dj;

There is no need to put this code in the for loop since the condition of x==j will only occur once so we can avoid the if condition on each iteration of the for loop.

Now that leaves the following code in the else

 if x==i
     U=V(:,:,x);
 else
     U=V(:,:,x)*U;                   
 end

Again, you don't need to check the condition x==i since the loop starts iterating at i, so you could do the initialization of U outside of the for loop. So the do stuff block now becomes

 V(:,:,j) = dj;
 U        = V(:,:,i);
 for x=i-1:-1:1
     U = V(:,:,x)*U;
 end

So the original 12 lines of code now becomes 5, and it may become a little more clear on how to vectorize this if you feel that it is still necessary at this point.

As an aside, you may want to consider renaming your i and j variables as these are also used to represent the imaginary number.

[Will]

Thanks Geoff - I really appreciate the time you've spent on this, a great explanation and answer