Sparse and binary mask

2018 9 14
3 个回答
回答已采纳
更新时间：2018 9 20
8 次查看（30 天）

0 个投票

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

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Greg
Greg 2018-9-14
In case II, is phase also sparse, or full?
D_coder
D_coder 2018-9-14
phase is full matrix
James Tursa
James Tursa 2018-9-14
What is the class of your image variable? Do you have a supported C compiler installed for a mex routine?
D_coder
D_coder 2018-9-14
class is double,C compiler is present

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Bruno Luong
Bruno Luong 2018-9-15
编辑：Bruno Luong 2018-9-15

1 个投票

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.

17 个评论
显示 15更早的评论 隐藏 15更早的评论

D_coder
D_coder 2018-9-15
编辑：D_coder 2018-9-15
This actually works , can you explain me why we should not use mask but index
D_coder
D_coder 2018-9-15
编辑：D_coder 2018-9-15
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
Walter Roberson 2018-9-15
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
D_coder 2018-9-15
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
Walter Roberson 2018-9-15
编辑：Walter Roberson 2018-9-15
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
D_coder 2018-9-15
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
Bruno Luong 2018-9-16
编辑：Bruno Luong 2018-9-16
"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
Bruno Luong 2018-9-16
编辑：Bruno Luong 2018-9-16
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
Bruno Luong 2018-9-16
编辑：Bruno Luong 2018-9-16
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
D_coder 2018-9-16
i cannot understand what col does
Bruno Luong
Bruno Luong 2018-9-16
col is column number
D_coder
D_coder 2018-9-16
what does col = floor((idx-1)/m))+1; do
Bruno Luong
Bruno Luong 2018-9-16
compute the column numbers
D_coder
D_coder 2018-9-16
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
D_coder 2018-9-16
编辑：D_coder 2018-9-16
essentially my matrix changes everytime when i do matrix = matrix.*phase , I don't think my matrix will update in accumarry
D_coder
D_coder 2018-9-20
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
Bruno Luong 2018-9-20
编辑：Bruno Luong 2018-9-20
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)

请先登录，再进行评论。

更多回答（2 个）

James Tursa
James Tursa 2018-9-15
编辑：James Tursa 2018-9-15

2 个投票

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++;
}
}

8 个评论
显示 6更早的评论 隐藏 6更早的评论

D_coder
D_coder 2018-9-15
i have never used mex before can you show me to run this code along with sample values
James Tursa
James Tursa 2018-9-15
>> mex times_threshold.c % <-- one time only
>> result = times_threshold(matrix,phase,threshold);
D_coder
D_coder 2018-9-15
编辑：D_coder 2018-9-15
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
James Tursa 2018-9-15
编辑：James Tursa 2018-9-15
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
D_coder 2018-9-15
there is no question of threshold in my code as I already do it before entering the loop
Walter Roberson
Walter Roberson 2018-9-15
James's code combines the threshold operation for efficiency.
D_coder
D_coder 2018-9-15
it wont be efficient as in my case i calculate threshold only once
James Tursa
James Tursa 2018-9-16
What code are you using that is "more efficient" than the mex routine?

请先登录，再进行评论。

Walter Roberson
Walter Roberson 2018-9-15

1 个投票

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

9 个评论
显示 7更早的评论 隐藏 7更早的评论

D_coder
D_coder 2018-9-15
some parts of phase may be zero some not. So I should go with the first code you posted right.
D_coder
D_coder 2018-9-15
编辑：D_coder 2018-9-15
how about if I use if M ~= 0 in the first code
Walter Roberson
Walter Roberson 2018-9-15
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
D_coder 2018-9-15
and if M is equivalent to if (M ~= 0) right?
Walter Roberson
Walter Roberson 2018-9-15
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
D_coder 2018-9-15
in your code there should be two ends
D_coder
D_coder 2018-9-15
My matrix is complex hence I am getting the error Complex values cannot be converted to logicals when i use if M
Walter Roberson
Walter Roberson 2018-9-15
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
D_coder 2018-9-15
this is not working it is taking very long time

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by