matlab算法进入死循环

Question

clear all;
%------------------------产生源信号并显示----------------------------------
t=linspace(0,100,100);
s1=sin(t/2);
plot(s1)
s2=1/1000*t.^2.*cos(t/4).^2;
figure;
plot(s2)
s3=exp(-0.2*t).*cos(pi*t);
figure;
plot(s3);
%-------------进行源信号混合---------------------------------------
A=[0.1,0.2,0.7;0.2,0.4,0.4;0.3,0.2,0.5]; %混合矩阵
s=[s1;s2;s3];
Y=A*s;
%-----------去均值和球化-------------------------------------------
meanValue=mean(Y');
Y=Y-diag(meanValue)*ones(size(Y));
covarianceMatrix=cov(Y');
[U,V,W]=svd(covarianceMatrix);
x=diag(V);
x=x.^(-1/2);
V=diag(x);
Z=V*U'*Y; %球化后的矩阵
%------------显示球化后的信号-----------------------------
figure;
plot(Z(1,:));
figure;
plot(Z(2,:));
figure;
plot(Z(3,:));
%---------------------盲源分离算法主程序----------------------------
syms y G %定义y和G为符号变量
G=1/4*y^4;
w1=rand(3,1);
w1=w1./sqrt(sum(w1.^2)); %产生初始的范数为1的随机向量
u=zeros(3,1);
g=diff(G);
for i=1:100
y=w1'*Z(:,i);
u=u+eval(g)*Z(:,i);
end
a=0.1;
u=-u./100-a*w1;
d1=-u; %产生初始的搜索方向向量d
p=zeros(3);
%-----选择初始的步长minb-----------------
for i=1:3
n=0;
c=a.^[0,1,2,3,4];
e=length(c);
while(n<=1000)
minb=c(1);
w2=w1+minb*d1;
h=0;
for k=1:100
y=w2'*Z(:,k);
h=h+eval(G);
end
minF=-h/100-a/2*(norm(w2)^2-1);
for k=2:e
w2=w1+c(k)*d1;
h=0;
for r=1:100
y=w2'*Z(:,r);
h=h+eval(G);
end
h=-h/100-a/2*(norm(w2)^2-1);
if(h<minF)
minb=c(k);
minF=h;
end
end
%-------对迭代求得的向量进行正交归一，进行下一次迭代-------------
w2=w1+minb*d1; %w2为经过迭代求得的向量
w3=zeros(3,1);
for l=1:i-1
w3=w3+(w2'*p(:,l))*p(:,l);
end
w2=w2-w3; %w2进行正交化
w2=w2./sqrt(sum(w2.^2));
u1=zeros(3,1);
for j=1:100
y=w2'*Z(:,j);
u1=u1+eval(g)*Z(:,j);
end
u1=-u1./100-a*w2;
t=-norm(u1)^2/(d1'*u);
b=max(t,0);
d2=-u1+b*d1;
n=n+1;
d1=d2;
w1=w2;
u=u1;
end
p(:,i)=w1; %把求得满足条件的向量存储在p中
q=w1'*Z;
figure;
plot(q); %绘制解混后的信号
end

算法跑不完，我每次终止都在这里
Operation terminated by user during sym/eval (line 11)
In Untitled3 (line 66)
h=h+eval(G);

Answer 1

并不是死循环，只是循环的次数比较多，你没有等到它算完而已。

 

当你强制中断时，也未必每次都在同一个地方，比如这行代码

h=h+eval(G);

就有两处，一处在 for k=1:100 循环内，另一处在 for k=2:e，for r=1:100 的循环内。除了这两行之外，还有一行比较耗时间，是 u1=u1+eval(g)*Z(:,j)。

 

这三行代码都是调用了符号运算，而且循环次数比较多（粗略估算，分别为30万、120万、30万次），差不多99%的运算时间都消耗在这三行代码上了，所以，当你强制中断时，刚好落在这三行上的概率非常大。

 

再耐心等一会儿，应该能够算完的。不过，就代码而言，应该可以进一步优化的，不是必须要等那么久的。

追问

确实跑完了，感觉跑了20多分钟。能给出优化方案吗？跪求。

追答

程序中有些循环可以改成向量化运算，估计效率应该可以获得数倍的提升。

 

不过，我觉得要想提高效率，最主要的还是要在大量循环中尽量避免使用符号运算，可以考虑用匿名函数代替eval，也就是把原代码中有关的几行：

g=diff(G);
...
             h=h+eval(G);
...
                 h=h+eval(G);
...
             u1=u1+eval(g)*Z(:,j);

改成下面这样（增加两行，修改3行）：

g=diff(G);
F=str2func(['@(y)' vectorize(char(G))]);
f=str2func(['@(y)' vectorize(char(g))]); 
...
            h=h+F(y);
...
                h=h+F(y);
...
            u1=u1+f(y)*Z(:,j);

 修改前，在我的电脑上运行一次需要1300多秒，而改后只需要20多秒，提高50倍以上。