Matlab:如何实现动力系统和时间序列的方程式
Matlab: How to implement an equation for a dynamical system and time series
我很难理解使用时间序列概念对生物医学图像进行聚类和分割的技术。该问题所依据的论文是:M。拉科米等。 al,乳腺 X 光图像分割基于
混沌图聚类算法 download link.
在一个D维space中有一组N个点{r_i}。 [0,1] 中的实变量 x_i 被分配给每个点和成对交互 J_ij = exp[-(r_i - r_j)^2 / 2a^2] 其中 a 是局部长度尺度。系统的时间演化由
给出
函数f
与神经网络中的激活函数非常相似。
logistic map 是一个单峰单变量时间离散非线性动力系统。
我正在寻找一种更快、更有效(矢量化)的方法来应用 Eq(1)
N = 100 万个图像特征点。 t 跨越 1 到 10 次演化。我做的方式但不确定代码是否正确。我随机生成一个维度为 D = 50 并包含 100 个数据点的矩阵 R。
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)。
我认为我无法就此给出一个简单的最终答案,但我绝对可以给您一些有用的提示:
提前计算一次矩阵J
。其中的信息是静态的,所以你不想重新计算。 Ci
.
同上
sum(Jij*yj)
部分实际上是矩阵与向量的乘积。这可以用线性代数最快地完成,即 J*y
.
向量化函数 f
:与其分别对每个元素执行 f(x_i),不如一次对所有元素执行 f(x)。例如。 f(x) = 1-2*x.^2
与 .^
而不是 ^
对 x 的每个元素执行幂运算符。
您可能希望将时间迭代作为外循环。这是您需要执行的唯一顺序计算。所有其余的,你想尽可能同时(〜并行)使用矢量化和线性代数。
这应该会给您一个很好的起点。如果之后还有其他需要帮助的东西,请更新您的问题或发表评论。在这一点上,这是我所能做的。部分原因是您的代码示例不是很 clear/descriptive。如果对现在的情况有特殊原因,也许您想添加评论。
祝你好运!
编辑:
代码示例:
% 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
,但为了清楚起见,我在本示例中使用了它们。
我很难理解使用时间序列概念对生物医学图像进行聚类和分割的技术。该问题所依据的论文是:M。拉科米等。 al,乳腺 X 光图像分割基于 混沌图聚类算法 download link.
在一个D维space中有一组N个点{r_i}。 [0,1] 中的实变量 x_i 被分配给每个点和成对交互 J_ij = exp[-(r_i - r_j)^2 / 2a^2] 其中 a 是局部长度尺度。系统的时间演化由
函数f
与神经网络中的激活函数非常相似。
logistic map 是一个单峰单变量时间离散非线性动力系统。
我正在寻找一种更快、更有效(矢量化)的方法来应用 Eq(1)
N = 100 万个图像特征点。 t 跨越 1 到 10 次演化。我做的方式但不确定代码是否正确。我随机生成一个维度为 D = 50 并包含 100 个数据点的矩阵 R。
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)。
我认为我无法就此给出一个简单的最终答案,但我绝对可以给您一些有用的提示:
提前计算一次矩阵
J
。其中的信息是静态的,所以你不想重新计算。Ci
. 同上
sum(Jij*yj)
部分实际上是矩阵与向量的乘积。这可以用线性代数最快地完成,即J*y
.向量化函数
f
:与其分别对每个元素执行 f(x_i),不如一次对所有元素执行 f(x)。例如。f(x) = 1-2*x.^2
与.^
而不是^
对 x 的每个元素执行幂运算符。您可能希望将时间迭代作为外循环。这是您需要执行的唯一顺序计算。所有其余的,你想尽可能同时(〜并行)使用矢量化和线性代数。
这应该会给您一个很好的起点。如果之后还有其他需要帮助的东西,请更新您的问题或发表评论。在这一点上,这是我所能做的。部分原因是您的代码示例不是很 clear/descriptive。如果对现在的情况有特殊原因,也许您想添加评论。
祝你好运!
编辑:
代码示例:
% 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
,但为了清楚起见,我在本示例中使用了它们。