Matlab: How to implement an equation for a dynamical system and time series

Question

I am having difficulty in understanding a technique for clustering and segmentation of biomedical images using the concept of time series. The paper on which the Question is based is : M. Lacomi et. al, Mammographic images segmentation based on chaotic map clustering algorithm download link.

There are a set of N points {r_i} in a D dimensional space. A real variable x_i in [0,1] is asigned to each point and the pair wise interactions J_ij = exp[-(r_i - r_j)^2 / 2a^2] where a is the local length scale. The time evolution of the system is given by

The function f has a close resemblance to the activation function in neural network. It is the logistic map which is a unimodal and univariate discrete in time non-linear dynamical system. I am looking for a faster and efficient (vectorized) way to apply Eq(1) when N = 1 million points that are features of images . t spans from 1 to 10 time evolutions. The way I did but unsure if the code is proper. I am randomnly generating a matrix R of dimension D = 50 and containing 100 data points.

N = 100;
D = 50;
T  =10;
R = rand(N,D);
x = zeros(N,T);
y(1) = rand();
for i = 1:N   %// for loop indicating the number of sample points
    y(i+1) = 1-2*y(i)^2; %/* the iterations of the map f */
    r_1(i) = R(i,:);
    r_2(i) = R(i+1,:);
    sum_j = 0.0;
    for t = 1:T
        x(i,:)= y;
        a = var(r(i));
        J = {exp(-(r_1(i) - r_2(i+1))^2)}/2a;
        sum_j = sum_j+J*(1-2*x(i+1,t)));
        x(i,t) =  (1/c(i))*(sum_j);
    end
end

A small implementation using matrix where each row is the data element and columns are the the dimensions will be very helpful in order to expand the code for multi dimensional images. I am having a tough time to code Eq(1).

Answer 1

I don't think I can give an easy final answer on this one, but I can definitely give you a few helpful tips:

1. Calculate the matrix J once in advance. The information in it is static, so you don't want to recalculate. Idem ditto with Ci.

2. The part sum(Jij*yj) is actually the product of a matrix with a vector. This can be done fastest with linear algebra, i.e. J*y.

3. Vectorize the function f: instead of doing f(x_i) for each element separately, do f(x) for all at once. E.g. f(x) = 1-2*x.^2 with .^ instead of ^ performs the power operator on each element of x.

4. You probably want to put the time-iteration as the outer loop. This is the only sequential calculation you need to do. All the rest, you want to do as much as possible simultaneously (~ parallel) using vectorization and linear algebra.

This should give you a good starting point. If there is other stuff that needs help afterwards, update your question or give a comment. At this point, it's all I can do. Partly, because your code sample isn't very clear/descriptive. Maybe you want to add comments if there are particular reasons for how it is now.

Good luck!

EDIT:

Code sample:

% Jmod is a modification of the matrix J:
%    1. Jmod(i,j) = J(i,j)/C(i) ==> the division by Ci is included
%    2. Jmod(i,i) = 0 ==> the diagonal elements are zero such that the term
%       for i=j is not included in the sum.

% Memory allocation
x = zeros(N,T+1); 
y = zeros(N,T+1);

% Initialization with your choice of x0
x(:,1) = x0;

% Time iterations
for t=1:T 
    y(:,t) = 1 - 2*x(:,t).^2;
    x(:,t+1) = Jmod*y(:,t);         
end

You don't really need both vectors x and y, but I have used them in this sample for clarity.