How to vectorize a for loop but with conditionals inside it ?

How to vectorize a for loop but with conditionals in it ?

    [sizMat1 sizMat2] = size(matrixToMultiply); 

    cumulMatProduct = ones(sizMat1,1); %stores the cumulative Products of chosen Matrices. 
    %gets updated at every iteration

    for ix = 2:length(col1)

        % Depending on if the value is either 0 or 1, pick a matrix;
         if (col1(ix) == 0 )
             cumulProduct = simpleMatrix0 * cumulMatrixProduct;
             matrixToMultiply = matrix1;

         elseif (col1(ix) == 1 )
             matrixToMultiply = matrix2;
         end

         anotherMatrixtoMultiply = diag( exp(constantMatrix) * col2(ix) ); 
         % Another Matrix is created by multiplying a scalar 
         %(picked from the same index ix of a different column col2 having same dimensions as col1)  

         cumulMatrixProduct = matrixToMultiply*anotherMatrixtoMultiply*cumulMatrixProduct; 

    end

[Azzi Abdelmalek]

Is this a complete code?

[nah]

Yes, almost. I could add:

      % matrixToMultiply  is 101 x 101
      % constantMatrix is 101 by 1 
      % anotherMatrixtoMultiply is 101 by 101 
      % cumulMatrixProduct = 101 x 1 (Result )

[Jan]

@nah: Vectorizing is not trivial. Inserting a bug by a typo is as easy and suggesting a method is is much slower than the original version due to the need of huge temporary arrays. Therefore questions about improving the speed require testdata, e.g. created by rand() if this is meaningful, such that we can run and test our ideas before posting.

Should we add this:

matrixToMultiply = rand(101, 101);
constantMatrix = rand(101, 1 );
anotherMatrixtoMultiply = rand(101, 101 );
cumulMatrixProduct = zeros(101, 1);  % Pre-allocation?!

[Jan]

@Readers: Is there a faster was to calculate A * diag(x) * B? I cannot test if A * bsxfun(@times, x, B) is faster than the BLAS librarie's DGEMM with the full matrix.

[nah]

Yes, the above test Data with rand should work. If needed, I can post my data as well (possibly by email or git )

Answer 1

[Jan]

I assume that the bottleneck of your code is: exp(constantMatrix), because exp() is very expensive. So move this calculation out of the loops:

expConstantMatrix = exp(constantMatrix);
for ix = 2:length(col1)
  ...
  anotherMatrixtoMultiply = diag( expConstantMatrix * col2(ix) );

Is this really wanted:

if (col1(ix) == 0 )
  ...
elseif (col1(ix) == 1 )
  ...
end

Or do you mean

if col1(ix) == 0
  ...
else
  ...
end

This might be the same in a standard case, but it is prone to bugs to add an elseif without an else.

[EDITED] Some tests show that this is faster:

anotherVectortoMultiply = expConstantMatrix * col2(ix); 
cumulMatrixProduct = matrixToMultiply * bsxfun(@times, anotherVectortoMultiply, cumulMatrixProduct);

[nah]

Thanks for the suggestions. Yes, the exponential is an expensive op indeed and I've moved it outside. I have an else added already, didn't include it here. Didn't know it could be a bug :-) So, because of the dependency, vectorization is almost impossible I think

[Jan]

If the multiplications of the matrices are the main work, a vectorization of the loop will not be useful. If large temporary arrays must be created, vectorized code can be slower than the loops. So do not try to optimize code, which is not the bottleneck.

[nah]

Actually, the motivation was not to optimize this particular piece, but rather be able to do it in a parallel fashion especially as gpuArray for which converting vectorized code is easier than working with for loops. ( i have a parfor version of the code above, but getting to do such calculations en masse on GPUs is where am stuck )