# Set of matrix operations for summing around a matrix element

I have a results of simulations, which are discrete complex numbers representing wave function on a NxN grid. I calculate phase of the wavefunction. For every point on that grid I have to perform this set of operations. Example:

|a b c|
|d e f|
|g h i|


Sum of difference between neighbours around element e:

 (a-b)+(b-c)+(c-f)+(f-i)+(i-h)+(h-g)+(g-d)+(d-a)


But every substraction δ have is subjected to this conditions

if δ > π than δ = δ - 2π
if δ < -π than δ = δ + 2π


I wrote a code in Octave, which is extremaly inefficient and slow, but I was hoping there is a way to move from buteforce calculations to faster calculations using some of the functions in image package. I know I can perfom sumation using convolution but I still have problem with other operations.

Code in Octave:

function deltaPhi = phaseDifference(phi1, phi2)
deltaPhi = phi1 - phi2;
if(deltaPhi > pi)
deltaPhi = deltaPhi - 2*pi;
endif

if(deltaPhi < -pi)
deltaPhi = deltaPhi + 2*pi;
endif;
end

function [phase] = checkPhase(M)

phase = zeros(size(M)-2);

for i = 2:size(M,1)-1
for j = 2:size(M,2)-1
phase(i-1,j-1) = phaseDifference(M(i-1,j-1),M(i,j-1)) + phaseDifference(M(i,j-1),M(i+1,j-1)) + phaseDifference(M(i+1,j-1),M(i+1,j)) + phaseDifference(M(i+1,j),M(i+1,j+1)) + phaseDifference(M(i+1,j+1),M(i,j+1)) + phaseDifference(M(i,j+1), M(i-1,j+1)) + phaseDifference(M(i-1,j+1), M(i-1,j)) + phaseDifference(M(i-1,j), M(i-1,j-1));
endfor
endfor

end


The idea is to rewrite this code in OpenCV and use some of that librarys methods.

• You're code is slow in GNU Octave because you write your code that way. It can be much, much faster if you write efficient code in Octave.
– Andy
Jan 30 '18 at 8:39
• @Andy Could you maybe write something that could speed it up, make it more efficient? Jan 31 '18 at 10:53

First some simple things:

• I would use end instead of endfor and endif, to keep the code compatible with MATLAB.

• It is not necessary to end a function with end (though it's not harmful either). The end is only necessary when writing nested functions that have access to the parent function's workspace.

• You should break up very long lines of code using ... (see an example below).

• When looping through a large matrix, always make the first index your inner loop. Matrix elements M(1,j) and M(2,j) are adjacent in memory, M(i,1) and M(i,2) are not. If you switch your two nested loops, your code will be faster for large arrays, as you use the cache better.

Now for the hard part: changing your logic.

|a b c|
|d e f|
|g h i|


Sum of difference between neighbours around element e:

(a-b)+(b-c)+(c-f)+(f-i)+(i-h)+(h-g)+(g-d)+(d-a)


Note that you compute the phase difference between two neighboring cells several times: for example, (a-b) will be computed when determining your value for e, but also in the previous loop iteration, when d was the center element; it will be computed two more times when that pair appears at the bottom of such a neighborhood.

Thus, you should pre-compute these differences and store them:

vdiff = phaseDifference(diff(M,1,1));
hdiff = phaseDifference(diff(M,1,2));


For this to work, phaseDifference must be rewritten, see below.

Now, vdiff has one fewer rows, and hdiff has one fewer columns, than M. We need to take this into account when indexing. Think of these as the elements in between the rows/columns of M.

The sum for one element at (i,j) now becomes:

phase(i-1,j-1) = hdiff(i-1,j-1) + hdiff(i-1,j) ...
- hdiff(i+1,j-1) - hdiff(i+1,j) ...
- vdiff(i-1,j-1) - vdiff(i,j-1) ...
+ vdiff(i-1,j+1) + vdiff(i,j+1);


The full function now becomes:

function phase = checkPhase2(M)
vdiff = phaseDifference(diff(M,1,1));
hdiff = phaseDifference(diff(M,1,2));
phase = zeros(size(M)-2);
for j = 2:size(M,2)-1
for i = 2:size(M,1)-1
phase(i-1,j-1) = hdiff(i-1,j-1) + hdiff(i-1,j) ...
- hdiff(i+1,j-1) - hdiff(i+1,j) ...
- vdiff(i-1,j-1) - vdiff(i,j-1) ...
+ vdiff(i-1,j+1) + vdiff(i,j+1);
end
end


Since Octave is slow with loops, it is likely that this vectorized version will be faster:

function phase = checkPhase3(M)
vdiff = phaseDifference(diff(M,1,1));
hdiff = phaseDifference(diff(M,1,2));
j = 2:size(M,2)-1;
i = 2:size(M,1)-1;
phase = hdiff(i-1,j-1) + hdiff(i-1,j) ...
- hdiff(i+1,j-1) - hdiff(i+1,j) ...
- vdiff(i-1,j-1) - vdiff(i,j-1) ...
+ vdiff(i-1,j+1) + vdiff(i,j+1);


I'm testing on MATLAB R2017a, and the vectorized form is actually slower than the version with loops, because indexing is still quite slow, whereas MATLAB has made huge strides in the last year speeding up code with loops.

Timings for an input of 2000x3000 elements:

• checkPhase(M) (the OP's version) runs in 1.3s,
• and swapping the loop order it runs in 1.0s,
• checkPhase2(M) runs in 0.17s,
• checkPhase3(M) runs in 0.25s.

Your timings on Octave will be vastly different (much slower for the code with loops).

For phaseDifference to work with an array input, it cannot use if the way that the original function does. It is possible to rewrite the same logic using logical indexing to work on a full matrix, but this version is simpler:

function deltaPhi = phaseDifference(deltaPhi)
deltaPhi = mod(deltaPhi + pi, 2*pi) - pi;


Note I replaced the two inputs with a single one, as diff already computes the difference between neighbors.

• Thank you very much this insight, it's like the code achived warp speed. Feb 1 '18 at 14:39
• I have few questions, why do substruct some of the values when you calculate phase? Also do you take into account issue of boundry conditions? Feb 1 '18 at 15:15
• @Glock: This code produces exactly the same results as the code you posted (to numerical accuracy, the order of operations is different, hence you get different rounding errors, to the order of 1e-15 or so), so no, it does not take into account boundary conditions. The output is slightly smaller than the input to avoid those boundary issues. phaseDifference uses mod(x,2*pi), which returns a value in the range [0,2*pi). By subtracting pi first, applying mod, then adding pi, you constrain to the range [-pi,pi). Feb 1 '18 at 16:40
• Thanks. I normaly extend the input matrix by two rows and two columns and than add periodic BC, so it work perfectly with your code. I checked against my data and results are correct. I'm just trying to understand why is it neccesery to substruct as in hdiff(i-1,j-1) + hdiff(i-1,j) ... - hdiff(i+1,j-1) - hdiff(i+1,j) ... Feb 1 '18 at 17:06
• Now I understand, diff substracts from left to right, so you have to substract hdiff(i+1,j-1) in order to keep correct sign. Feb 1 '18 at 17:13