How to improve matrix indexing efficiency or avoid indexing in Finite difference?

Hi all,

I am trying to optimize the performance of a FDTD (Finite-differnce-time-domain) code implemented in Matlab. As you may know, FDTD usually uses finite difference stencils in a time loop, i.e., the code calculates for example A(2:end,:)-A(1:end-1,:) many places in it. I've found that most of the time Matlab is spending on is not the actual subtracting or adding numbers, but the indexing.

For a very simple example, if I run this:

clear;
clc;
close all;
N1=100;
N2=100;
A=(ones(N1,N2,'single'));
C=(ones(N1,N2,'single'));
tic;
for i=1:1e4
    %C(2:end,:)=C(2:end,:)-A(1:end-1,:);
    C=C-A;
end
toc;

I got Elapsed time is 0.056711 seconds.

Instead if I run the following:

clear;
clc;
close all;
N1=100;
N2=100;
A=(ones(N1,N2,'single'));
C=(ones(N1,N2,'single'));
tic;
for i=1:1e4
    C(2:end,:)=C(2:end,:)-A(1:end-1,:);
    %C=C-A;
end
toc;

I got: Elapsed time is 0.316735 seconds.

That is to say the most of the time Matlab is just doing the matrix indexing, Is there any way to improve this or avoid the indexing but get the same result? This could make my code almost 10 times faster!

Answer 1

[Matt J]

the code calculates for example A(2:end,:)-A(1:end-1,:) many places in it.

For simple cases like this, you should be using the built-in DIFF command,

    D=diff(A,1,1); %equivalent to A(2:end,:)-A(1:end-1,:)

That is to say the most of the time Matlab is just doing the matrix indexing

It is not the indexing itself that is slow. It is that the indexing statements on the right hand side result in the creation of temporary arrays. Indexing on the left hand side should not have this problem, e.g., something like,

for i=1:1e4
    C(2:end,:)=X-Y; 
end

You can work around this this by doing right-hand side indexing outside of the loop, so that temporary arrays will be created only once:

tic;

Ctmp=C(2:end,:); Atmp=A(1:end-1,:);

for i=1:1e4
  Ctmp=Ctmp-Atmp;
end

C(2:end,:)=Ctmp; A(1:end-1,:)=Atmp;

toc

If you have more complicated indexing that varies from iteration to iteration of the loop, you could also try this FEX file, although I haven't used it myself.

[Hao Zhang]

Hi,

Many thanks for your answer. some remarks and other observations: 1. the indexing itself does adding considerable amount of runtime. for example: run this: clear;clc;close all;

N1=100;N2=100;
A=(ones(N1,N2,'single'));
C=(ones(N1,N2,'single'));
X=(ones(N1-1,N2-1,'single'));
Y=(ones(N1-1,N2-1,'single'));
tic;
for i=1:1e4
    %a=diff(A,1,1);
    %C(2:end,1:end-1)=C(2:end,1:end-1)+diff(A(:,1:end-1),1,1);
    %C=C-A;
    %C(2:end,1:end-1)=X-Y;
    a=X-Y;
end
toc;

Elapsed time is 0.052130 seconds.

And if run this:

clear;clc;close all;
N1=100;N2=100;
A=(ones(N1,N2,'single'));
C=(ones(N1,N2,'single'));
X=(ones(N1-1,N2-1,'single'));
Y=(ones(N1-1,N2-1,'single'));
tic;
for i=1:1e4
    %a=diff(A,1,1);
    %C(2:end,1:end-1)=C(2:end,1:end-1)+diff(A(:,1:end-1),1,1);
    %C=C-A;
    C(2:end,1:end-1)=X-Y;
    %a=X-Y;
end
toc;

Elapsed time is 0.169469 seconds.

The diff is a really good idea, it is even faster than no indexing, so if I run this:

clear;clc;close all;
N1=100;N2=100;
A=(ones(N1,N2,'single'));
C=(ones(N1,N2,'single'));
X=(ones(N1-1,N2-1,'single'));
Y=(ones(N1-1,N2-1,'single'));
tic;
for i=1:1e4
    a=diff(A,1,1);
    %C(2:end,1:end-1)=C(2:end,1:end-1)+diff(A(:,1:end-1),1,1);
    %C=C-A;
    %C(2:end,1:end-1)=X-Y;
    %a=X-Y;
end
toc;

Elapsed time is 0.029310 seconds.

But the actual content I have to use many places in the FDTD code is something like

C(2:end,1:end-1)=C(2:end,1:end-1)+A(2:end,1:end-1)-A(1:end-1,1:end-1);

Because usually the boundary needs special treatment such as periodic boundary condition.

If I run this I got Elapsed time is 0.486443 seconds.

By using diff, it is slightly faster, so if I run

C(2:end,1:end-1)=C(2:end,1:end-1)+diff(A(:,1:end-1),1,1);

Then I get Elapsed time is 0.401249 seconds.

I can properly use separated matrices for the inner nodes and the boundaries to avoid the extra indexing, but that will be really complicated and reduce the readability of the code...

I was wondering what Matlab exactly is doing when it sees the indexing? It creates a temporary variable every time it sees an indexing? That seems very stupid to me, couldn't it just get the address of the part of the matrix in the memory and get the numbers out? Or is it a fundamental limit that impossible to optimize anymore?

Any answers or further improvement of the above code are welcome, thanks!!

[Matt J]

See here,

https://www.mathworks.com/matlabcentral/answers/106903-which-matlab-operations-functions-need-speeding-up#comment_558320

[Hao Zhang]

thanks!