# How to reduce computational time in DDE simulation on Matlab

I need to simulate a network of nodes. The weights of the edges are being given in a matrix . Due to the non-zero distances between the nodes, we consider time-delays, which are computed given the distances and an arbitrary velocity. The delay values have already been calculated. Given that $$A$$ is symmetric, for a $$n$$ by $$n$$ matrix, we have $$m=\frac{n(n-1)}{2}$$ different time-delays.

load('A.mat'); % A is a symmetric matrix of dimensions n by n
load('tau.mat'); % tau is the delay array vector with n(n-1)/2 components
tf = 10;
sol = dde23(@(t,y,z)(dde(t,y,z,A)),tau,@(t)(history(t,A)),[0 tf]);
t = linspace(0,tf,200);
y = deval(sol,t);
figure
plot(t,y)


The network needs to be expressed in its analytic form as a dynamical system:

$$\dot{x}=A_{0}x+A_{1}x(t-\tau_{1})+...+A_{m}x(t-\tau_{m})$$

, where $$A_{0}$$ is the diagonal matrix of the initial $$A$$ matrix, where it has been considered that since the distance between a node and itself is zero, there is no delay. In the absence of a way to initialize $$m$$ different matrices of $$n$$ by $$n$$, i created a $$P$$ matrix as below. Each symmetric pair of $$A$$ matrix would be put in its respective place of $$P$$ and multiplied with the respective delay.

function dydt = dde(t,y,yd,A)
n = length(A);
A0 = diag(diag(A));
m = n*(n-1)/2;
P = zeros(n,n*m);
k = 0;
for i = 1:n
for j = 1:n
if j>i
P(i,j+n*k) = A(i,j);
P(j,i+n*k) = A(i,j);
k = k+1;
end
end
end
dydt = A0*y;
for c = 0:m-1
dydt = dydt + P(:,n*c+1:(c+1)*n)*yd(:,c+1);
end
end
function y = history(t,A)
y = rand(length(A),1);
end


If the code is technically correct, it is not efficient, since the time it does to produce results increase exponentially even for small number of nodes. For instance, for a network of 7 nodes it takes 2.1 seconds, for 8 nodes 7.6 seconds, for 9 nodes 30 seconds and for 10 nodes almost 121 seconds. This makes its use prohibitive for large networks.

Does this exponentially increasing computational cost make sense?

Is there anything to do to decrease the time of computation?

• The construction of the large, almost zero matrix $P$ is a significant overhead. We already had that discussion. It is also possible that with increasing complexity of the model the system becomes more stiff, requiring smaller step sizes. So eliminating $P$ will reduce the base number of function evaluation time, but not the dynamic of the run times. Apr 3 at 5:10
• So, if I get you correctly, I need to find an alternative way of creating the additional matrices. However, i am not certain if Matlab gives such an opportunity. Does that mean that I should quit trying this with Matlab and try another programming language that might allow initialization of an arbitrary number of matrices? Apr 3 at 15:30

To make my comments more explicit:

Try with the following code for the DDE function that eliminates P

function dydt = dde(t,y,yd,A)
n = length(A);
dydt = diag(A).*y;
k = 0;
for i = 1:n
for j = (i+1):n
dydt(i) += A(i,j)*yd(j,k);
dydt(j) += A(i,j)*yd(i,k);
k = k+1;
end
end
end


This eliminates the $$O(n^3)$$ matrix-vector multiplications of P with yd, the effort of filling the non-zero entries of P is about the same as directly updating the derivatives vector, the most costly is the explicit double loop.

The suspected result of this change is a general reduction by a factor 2 to 10(?) in the run-time and a somewhat flatter progression of the run-time dependence on n.

And yes, for experiments with larger configurations you should use compiled code, in matlab, or using JitCDDE in python (which does a code transformation to C or Fortran and uses the compiled routines) or using the respective julia-lang library.