Matlab:如何实现动力系统和时间序列的方程

Question

我很难理解一种利用时间序列的概念对生物医学图像进行聚类和分割的技术。这个问题所依据的论文是: M. Lacomi et。基于混沌映射聚类算法下载链接的乳腺图像分割。

在D维空间中有一组N点{r_i} .将0,1中的一个实变量x_i分配到每个点，并对成对交互作用J_ij = exp-(r_i - r_j)^2 / 2a^2，其中a是局部长度尺度。给出了系统的时间演化。

​

f函数与神经网络中的激活函数非常相似。物流地图是时间上的单峰单变量离散非线性动力系统.我正在寻找一种更快、更有效(矢量化)的方法来应用Eq(1)，当N=100万点时，这是图像的特征。T从1到10次演化。我这样做了，但不确定代码是否正确。我随机地生成维数D= 50的矩阵R，包含100个数据点。

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

一个使用矩阵的小型实现，其中每一行都是数据元素，列是维度，这对于扩展多维图像的代码非常有帮助。我很难编码Eq(1)。

Answer 1

我不认为我能给出一个简单的最后答案，但我绝对可以给你一些有用的建议：

1. 预先计算一次矩阵J。它中的信息是静态的，所以您不想重新计算。Ci的同上。
2. 部分sum(Jij*yj)实际上是带有向量的矩阵的乘积。用线性代数，即J*y，可以最快地做到这一点。
3. 对函数f进行向量化:不对每个元素分别执行f(x_i)，而是同时对所有元素执行f(x)。例如，f(x) = 1-2*x.^2与.^而不是^一起对x的每个元素执行幂运算符。
4. 您可能希望将时间迭代作为外部循环。这是您唯一需要进行的顺序计算。其余的，你想要尽可能多的同时(~并行)使用矢量化和线性代数。

这应该给你一个好的起点。如果有其他事情需要帮助，然后，更新你的问题或给出评论。在这一点上，这是我所能做的。部分原因是您的代码示例不是很清楚/描述性很强。也许你想添加评论，如果有特殊的原因，它现在是这样的。

祝好运!

编辑：

代码示例：

% 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

您实际上并不需要向量x和y，但为了清晰起见，我在这个示例中使用了它们。