Sparse and binary mask

I am working on saving computations and speeding up my Matlab code. Here's what I do Perform thresholding on an image and make most of its components as zero run it inside a loop with element by element multiplication with another matrix(full matrix) of the same size

Here is a detailed version of what I am doing

CASE I:When I dont avoid multiplication with zeros

 matrix(matrix<threshold) = 0 ;
 while the condition is true
 matrix = matrix.*phase;
 other lines of code ;
 end

Result: From MATLAB profiler I see that matrix = matrix.*phase; takes 80 seconds

CASE II:When I avoid multiplication with zeros using SPARSE

matrix(matrix<threshold) = 0 ;
matrix = sparse(matrix);
while the condition is true
matrix = matrix.*phase;
other lines of code ;
end

Result: From MATLAB profiler I see that matrix = matrix.*phase takes 110 seconds

CASE III:When I avoid multiplication with zeros using a mask

matrix(matrix<threshold) = 0 ;
mask = matrix~=0;
while the condition is true
matrix(mask) = matrix(mask).*phase(mask);
other lines of code ;
end

Result: From MATLAB profiler I see that matrix = matrix.*phase takes 130 seconds

Why is my code slow if I try to avoid multiplications with zeros? Since the matrix after thresholding has 75% zeros my speed should be up by 75% as I am avoiding multiplication with zeros. But my time is increased Is there any other way to achieve avoiding multiplication with zeros and speeding up of the code?

MY PRIORITY IS TO SPEED UP THE CODE EXECUTION AND SAVE MUTLIPLICATIONS

[Greg]

In case II, is phase also sparse, or full?

[D_coder]

phase is full matrix

[James Tursa]

What is the class of your image variable? Do you have a supported C compiler installed for a mex routine?

[D_coder]

class is double,C compiler is present

Answer 1

[Bruno Luong]

Instead of matrix(mask) = matrix(mask).*phase(mask)

you might work on index

idx = find(mask & matrix~=0);
while ...
   matrix(idx) = matrix(idx).*phase(idx);
end

You might want to compute phase(idx) outside the loop if they do not change, or update idx if the matrix elements become zero. As you don't give any detail, we have no idea if those are applicable for you.

If matrix does not share any copy, you might do element inplace replacement using MEX.

[D_coder]

This actually works , can you explain me why we should not use mask but index

[D_coder]

after doing the multiplication operation , inside the while loop I am accumulating the sum of matrix in a different element.

while 
matrix(idx) = matrix(idx).*phase(idx);
new_Matrix(a,:) = integer*sum(matrix,1) ;%want to sum matrix ignoring zeros in each row

Is there any way to ignore zeros in each row and sum elements of each column without making the computation slower (by using same thing as you used above)

[Walter Roberson]

If you are referring to individual zeros, then it is unlikely that you could do any better than just adding them anyhow. Detecting them would require a comparison and conditional logic, which at the MATLAB level would be slower than just doing the additions. At the HDL level, maybe you could get a marginal improvement by OR'ing all of the bits together with gates to test for zero and using that as a control signal, but it would not be very much.

But for floating point on an x64 architecture, the tests would be slower than not doing the tests. See https://stackoverflow.com/questions/23203710/floating-point-operations-per-cycle-intel

[D_coder]

I was thinking of same indices idx in the above code already calculated and use it for calculating sum if i do

integer*sum(matrix(idx),1)

it is superfast almost 9 times but problem is matrix(idx) is a vector and i want a matrix

[Walter Roberson]

sum(matrix(idx),1) would probably calculate a scalar. matrix(idx) would typically be a column vector when idx is determined by find() or logical indexing, and sum() of column vector along the first dimension is a scalar. This is not the same thing as integer*sum(matrix,1) but skipping the zeros: for 2D matrix, sum(matrix,1) is a row vector with as many entries as the number of columns in matrix.

[D_coder]

i know this, thats what i was saying matrix(idx) is a vector not a matrix so will result in a scalar. I guess the only way to tackle this issue is to store indices of each row of the matrix 'matrix'

[Bruno Luong]

"Is there any way to ignore zeros in each row and sum elements of each column without making the computation slower?"

If you store the matrix as sparse it will do it that way. But you might pay the price somewhere else.

[Bruno Luong]

You might split the calculation of sum(matrix,1) as two quantities:

matrix_base = matrix;
matrix_base(mask) = 0;
sumbase = sum(matrix_base,1);
matrix_change = matrix;
matrix_change(~mask) = 0;
idx = find(mask)
while ...
   matrix_change(idx) = matrix_change(idx).*phase(idx);
   sumchange = sum(matrix_change,1);
   summatrix = sumbase + sumchange ;
   ...
end

Again you might store matrix_change as sparse. The efficiency of such thing depends on the structure of the mask. You have to experience with your case to know what strategy is the best for you.

[Bruno Luong]

Another way is to avoid building matrix all together:

...
[m,n] = size(matrix);
col = floor((idx-1)/m))+1;
idx = idx(:);
phaseidx = phase(idx);
while ...
     sumchange = accumarray(col, matrix_change(idx).*phaseidx, [n,1]);
     summatrix = sumbase + sumchange.';
     ...
end

[D_coder]

i cannot understand what col does

[Bruno Luong]

col is column number

[D_coder]

what does col = floor((idx-1)/m))+1; do

[Bruno Luong]

compute the column numbers

[D_coder]

i also agree that accum array is an excellent way of doing it but the resultant matrix I am getting is way smaller than what sum gives me

new_Matrix(a,:) = integer*sum(matrix,1) %gives me new_Matrix size as equal to matrix

Why am I getting a different size?

[D_coder]

essentially my matrix changes everytime when i do matrix = matrix.*phase , I don't think my matrix will update in accumarry

[D_coder]

how does the accumarray function that you wrote changes when I want to do sum of individual columns instead of individual rows.The idx find creates a problem here as it finds index column wise

a = [ 0 0 3 0; 2 5 0 1;0 0 0 3];
b  = [ 4 3 2 1; 5 3 6 2; 7 1 -2 3];
[r c] = size(a);
prod = a.*b; 
row_sum = sum(prod,2);%must have same result as sum_accum
idx = find(a);
ared = a(idx);
bred = b(idx);
cred = ared.*bred;
row = ceil(idx/c);
sum_accum = accumarray(row,cred,[r 1]);%must have same result as row_sum
both sum_accum and row_sum give different result.

[Bruno Luong]

a = [ 0 0 3 0; 2 5 0 1; 0 0 0 3];
b  = [ 4 3 2 1; 5 3 6 2; 7 1 -2 3];
idx = find(a)

[m,n] = size(a);
col = floor((idx-1)/m)+1;
sumcol = accumarray(col,a(idx).*b(idx)).'
% check
prod = a.*b;
sum(prod,1)

row = mod((idx-1),m)+1;
sumrow = accumarray(row,a(idx).*b(idx))
% check
sum(prod,2)

Answer 2

[James Tursa]

Try this simple C mex routine. It avoids the overhead you are creating with your sparse alternate stuff. Could maybe be improved further by multi-threading for large variables but it is currently not coded that way. (Caution, no input checking so it will crash MATLAB if not used properly)

/*
 * File: times_threshold.c
 *
 * C = times_threshold(A,B,threshold)  does
 *
 * C = A .* B but only where A > threshold (C = 0 in the other spots)
 *
 * A = full real double variable
 * B = full real double variable (same size as A)
 * C = full real double variable (same size as A)
 * threshold = numeric scalar
 *
 * CAUTION: Bare bones with no argument checking! Will crash MATLAB if the
 *          inputs are not exactly as expected.
 */

#include "mex.h"
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *Apr, *Bpr, *Cpr;
    double threshold;
    mwSize n;

      n = mxGetNumberOfElements(prhs[0]);
      plhs[0] = mxCreateNumericArray(mxGetNumberOfDimensions(prhs[0]),
                                     mxGetDimensions(prhs[0]),
                                     mxDOUBLE_CLASS, mxREAL);
      Apr = mxGetPr(prhs[0]);
      Bpr = mxGetPr(prhs[1]);
      Cpr = mxGetPr(plhs[0]);
      threshold = mxGetScalar(prhs[2]);
      while( n-- ) {
          if( *Apr > threshold ) {
              *Cpr = *Apr * *Bpr;
          }
          Apr++; Bpr++; Cpr++;
      }
  }

[D_coder]

i have never used mex before can you show me to run this code along with sample values

[James Tursa]

>> mex times_threshold.c   % <-- one time only
>> result = times_threshold(matrix,phase,threshold);

[D_coder]

how will the above code change if I want it to run as

while a>b
matrix = matrix.*phaseshift
some other lines of code
end

Since my variables are large numel(matrix) = 1000000 how will your code change when multithreading is introduced

[James Tursa]

Simply this:

>> matrix = times_threshold(matrix,phase,threshold);

Unless you are asking how to modify matrix inplace?

For multi-threading, the source code would need to be rewritten.

[D_coder]

there is no question of threshold in my code as I already do it before entering the loop

[Walter Roberson]

James's code combines the threshold operation for efficiency.

[D_coder]

it wont be efficient as in my case i calculate threshold only once

[James Tursa]

What code are you using that is "more efficient" than the mex routine?

Answer 3

[Walter Roberson]

matrix<threshold and matrix~=0 both generate values the same size as matrix. matrix(mask) and phase(mask) and matrix(mask).*phase(mask) generate values with the same number of entries as the number of non-zero values in mask. Those variables get stored in memory and have to be processed.

At the MATLAB level, matrix = matrix.*phase gets dispatched to the high performance libraries that operate on multiple cores. One temporary value the same size as matrix would be involved to hold the results.

The high performance libraries would probably be scheduling multiple multiplies simultaneously, perhaps with SIMD type instructions, possibly able to graduate one result per cycle, or whatever the maximum graduation rate is.

When you select part of the array it has to be conditionally moved to other memory locations, and then that array processed, and then selectively moved back. That involves memory bandwidth and comparisons and loops. Those cost, and the tradeoffs are potentially just not worth the effort compared to SIMD or vector instructions.

You would face different tradeoffs in your FPGA.

To reduce multiplications, it might be worth creating an explicit loop, along the lines of

new_matrix = zeros(size(matrix));
for idx = 1 : numel(matrix)
  M = matrix(idx); P = phase(idx);
  if M && P; new_matrix(idx) = M .* P;
end

If phase is certain to be non-zero then

new_matrix = zeros(size(matrix));
for idx = 1 : numel(matrix)
  M = matrix(idx);
  if M; new_matrix(idx) = M .* phase(idx);
end

[D_coder]

some parts of phase may be zero some not. So I should go with the first code you posted right.

[D_coder]

how about if I use if M ~= 0 in the first code

[Walter Roberson]

if M && P

is equivalent to

if (M ~= 0) && (P ~= 0)

except in the case where M might be Nan, in which case the explicit comparison will work but the one without will fail. Though if you are expecting nans you might want to test for them, as nan times anything is always nan. On the other hand, if nan are a possibility, that drives up the die area required for FPGA: processing NaN is relatively expensive.

[D_coder]

and if M is equivalent to if (M ~= 0) right?

[Walter Roberson]

Correct,

if M

is equivalent to

if M ~=0

except in the case where M might have NaN. If NaN are present then

if M

will complain that NaN cannot be converted to logicals, but

if M ~= 0

will identify NaN as being non-zero.

[D_coder]

in your code there should be two ends

[D_coder]

My matrix is complex hence I am getting the error Complex values cannot be converted to logicals when i use if M

[Walter Roberson]

new_matrix = zeros(size(matrix));
for idx = 1 : numel(matrix)
  M = matrix(idx);
  if M ~= 0; new_matrix(idx) = M .* phase(idx); end
end

[D_coder]

this is not working it is taking very long time