ICA algorithm
2014-02-25 11:04
357 查看
ICA 即为 独立成分分析 类比于主成分分析PCA。
与PCA不同的是,他适用于非高斯性信号。典型应用背景为cocktail party effect。
![](http://simg.sinajs.cn/blog7style/images/common/sg_trans.gif)
algorithm" TITLE="ICA algorithm" />
相比于PCA,效果十分明显。
假设观测信号为Y(NxN),待提取信号为X,维度相同。Y=W*X;算法目标就是找到W逆,或是近似。一般实现基于两种思想,一是互信息最小,或是非高斯性最大。
互信息最小
互信息 I(X;Y)=H(X)-H(X|Y);而H(X|Y)=H(X,Y)-H(Y)
算法流程:
1>
随机初始化W(0)
2>
W(t+1)=W(t) + n(t)(I-f(Y)YT)W(t)
3> 不收敛,重复2
其中n为一般函数,确定W更新步长,f根据情况分,超高斯分布,f(Y)=Tanh(Y),亚高斯分布,f(Y)=Y-Tanh(Y)
2.非高斯性最大
既然假设服从非高斯分布,则需要指标衡量其非高斯性程度--负熵
N(X)=H(Xgaussian) - H(X)
,由于计算比较困难,通常只能逼近。
代码:
%%%%%%%%%%%%%%%%%%%%%%%%%% 初始化
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;clear all;close all;
%%%%%%%%%%%%%% 读入原始图像,混合,并输出混合图像
%%%%%%%%%%%%%%%%%%
% 读入混合前的原始图片并显示
t=0:1/100:9;
I1=sin(t);
I2=randn(1,901);
I3=square(4*t);
subplot(4,3,1),plot(I1),title('输入信号1'),
subplot(4,3,2),plot(I2),title('输入信号2'),
subplot(4,3,3),plot(I3),title('输入信号3'),
% 将其组成矩阵
S=[I1;I2;I3];
% 图片个数即为变量数,图片的像素数即为采样数
% 因此S_all是一个变量个数*采样个数的矩阵
Sweight=randn(size(S,1));
%
取一随机矩阵,作为信号混合的权矩阵
MixedS=Sweight*S;
%
得到三个混合信号矩阵
% 将混合矩阵重新排列并输出
subplot(4,3,4),plot(MixedS(1,:)),title('混合信号1'),
subplot(4,3,5),plot(MixedS(2,:)),title('混合信号2'),
subplot(4,3,6),plot(MixedS(3,:)),title('混合信号3'),
MixedS_bak=MixedS;
%
将混合后的数据备份,以便在恢复时直接调用
%%%%%%%%%%%%%%%%%%%%%%%%%% 标准化
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MixedS_mean=zeros(3,1);
for i=1:3
MixedS_mean(i)=mean(MixedS(i,:));
end
% 计算MixedS的均值
for i=1:3
for
j=1:size(MixedS,2)
MixedS(i,j)=MixedS(i,j)-MixedS_mean(i);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%% 白化
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MixedS_cov=cov(MixedS');
% cov为求协方差的函数
[E,D]=eig(MixedS_cov);
%
对图片矩阵的协方差函数进行特征值分解
Q=inv(sqrt(D))*(E)';
% Q为白化矩阵
MixedS_white=Q*MixedS;
%
MixedS_white为白化后的图片矩阵
IsI=cov(MixedS_white');
% IsI应为单位阵
%%%%%%%%%%%%%%%%%%%%%%%% FASTICA算法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X=MixedS_white;
%
以下算法将对X进行操作
[VariableNum,SampleNum]=size(X);
numofIC=VariableNum;
%
在此应用中,独立元个数等于变量个数
B=zeros(numofIC,VariableNum);
%
初始化列向量w的寄存矩阵,B=[b1 b2 ...
bd]
for r=1:numofIC
%
迭代求取每一个独立元
i=1;maxIterationsNum=100;
%
设置最大迭代次数(即对于每个独立分量而言迭代均不超过此次数)
IterationsNum=0;
b=rand(numofIC,1)-.5;
%
随机设置b初值
b=b/norm(b);
% 对b标准化
while
i<=maxIterationsNum+1
if i == maxIterationsNum
% 循环结束处理
fprintf('\n第%d分量在%d次迭代内并不收敛。', r,maxIterationsNum);
break;
end
bOld=b;
a2=1;
u=1;
t=X'*b;
g=(exp(2.*t)-1)./(exp(2.*t)+1);
dg=4*exp(2.*t)./(exp(2.*t)+1).^2;
% g的一阶导数
b=((1-u)*t'*g*b+u*X*g)/SampleNum-mean(dg)*b;
% 核心公式,参见理论部分公式2.52
b=b-B*B'*b;
% 对b正交化
b=b/norm(b);
if abs(abs(b'*bOld)-1)<1e-9
% 如果收敛,则保存b
B(:,r)=b;
break;
end
i=i+1;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%% 数据复原并构图
%%%%%%%%%%%%%%%%%%%%%%%%
ICAedS=B'*Q*MixedS_bak;
% 参见理论部分公式2.55
% 将混合矩阵重新排列并输出
subplot(4,3,7),plot(ICAedS(1,:)),title('ICA输出信号1'),
subplot(4,3,8),plot(ICAedS(2,:)),title('ICA输出信号2'),
subplot(4,3,9),plot(ICAedS(3,:)),title('ICA输出信号3'),
%%%%%%%%%%%%%%%%%%%%%%%%% PCA计算并构图
%%%%%%%%%%%%%%%%%%%%%%%%
[V,D]=eig(MixedS_cov);
Vtmp=zeros(size(V,1),1);
for j=1:2
for i=1:2
if D(i,i)
tmp=D(i,i);Vtmp=V(:,i);
D(i,i)=D(i+1,i+1);V(:,i)=V(:,i+1);
D(i+1,i+1)=tmp;V(:,i+1)=Vtmp;
end
end
end
t1=(MixedS'*V(:,1))';
t2=(MixedS'*V(:,2))';
t3=(MixedS'*V(:,3))';
subplot(4,3,10),plot(t1),title('PCA输出信号1'),
subplot(4,3,11),plot(t2),title('PCA输出信号2'),
subplot(4,3,12),plot(t3),title('PCA输出信号3'),
![](http://simg.sinajs.cn/blog7style/images/common/sg_trans.gif)
algorithm" TITLE="ICA algorithm" />
因为一般数据需要去相关性,通常要白化处理,白化推导过程如下:
![](http://simg.sinajs.cn/blog7style/images/common/sg_trans.gif)
algorithm" TITLE="ICA algorithm" />
与PCA不同的是,他适用于非高斯性信号。典型应用背景为cocktail party effect。
![](http://simg.sinajs.cn/blog7style/images/common/sg_trans.gif)
algorithm" TITLE="ICA algorithm" />
相比于PCA,效果十分明显。
假设观测信号为Y(NxN),待提取信号为X,维度相同。Y=W*X;算法目标就是找到W逆,或是近似。一般实现基于两种思想,一是互信息最小,或是非高斯性最大。
互信息最小
互信息 I(X;Y)=H(X)-H(X|Y);而H(X|Y)=H(X,Y)-H(Y)
算法流程:
1>
随机初始化W(0)
2>
W(t+1)=W(t) + n(t)(I-f(Y)YT)W(t)
3> 不收敛,重复2
其中n为一般函数,确定W更新步长,f根据情况分,超高斯分布,f(Y)=Tanh(Y),亚高斯分布,f(Y)=Y-Tanh(Y)
2.非高斯性最大
既然假设服从非高斯分布,则需要指标衡量其非高斯性程度--负熵
N(X)=H(Xgaussian) - H(X)
,由于计算比较困难,通常只能逼近。
代码:
%%%%%%%%%%%%%%%%%%%%%%%%%% 初始化
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;clear all;close all;
%%%%%%%%%%%%%% 读入原始图像,混合,并输出混合图像
%%%%%%%%%%%%%%%%%%
% 读入混合前的原始图片并显示
t=0:1/100:9;
I1=sin(t);
I2=randn(1,901);
I3=square(4*t);
subplot(4,3,1),plot(I1),title('输入信号1'),
subplot(4,3,2),plot(I2),title('输入信号2'),
subplot(4,3,3),plot(I3),title('输入信号3'),
% 将其组成矩阵
S=[I1;I2;I3];
% 图片个数即为变量数,图片的像素数即为采样数
% 因此S_all是一个变量个数*采样个数的矩阵
Sweight=randn(size(S,1));
%
取一随机矩阵,作为信号混合的权矩阵
MixedS=Sweight*S;
%
得到三个混合信号矩阵
% 将混合矩阵重新排列并输出
subplot(4,3,4),plot(MixedS(1,:)),title('混合信号1'),
subplot(4,3,5),plot(MixedS(2,:)),title('混合信号2'),
subplot(4,3,6),plot(MixedS(3,:)),title('混合信号3'),
MixedS_bak=MixedS;
%
将混合后的数据备份,以便在恢复时直接调用
%%%%%%%%%%%%%%%%%%%%%%%%%% 标准化
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MixedS_mean=zeros(3,1);
for i=1:3
MixedS_mean(i)=mean(MixedS(i,:));
end
% 计算MixedS的均值
for i=1:3
for
j=1:size(MixedS,2)
MixedS(i,j)=MixedS(i,j)-MixedS_mean(i);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%% 白化
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MixedS_cov=cov(MixedS');
% cov为求协方差的函数
[E,D]=eig(MixedS_cov);
%
对图片矩阵的协方差函数进行特征值分解
Q=inv(sqrt(D))*(E)';
% Q为白化矩阵
MixedS_white=Q*MixedS;
%
MixedS_white为白化后的图片矩阵
IsI=cov(MixedS_white');
% IsI应为单位阵
%%%%%%%%%%%%%%%%%%%%%%%% FASTICA算法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X=MixedS_white;
%
以下算法将对X进行操作
[VariableNum,SampleNum]=size(X);
numofIC=VariableNum;
%
在此应用中,独立元个数等于变量个数
B=zeros(numofIC,VariableNum);
%
初始化列向量w的寄存矩阵,B=[b1 b2 ...
bd]
for r=1:numofIC
%
迭代求取每一个独立元
i=1;maxIterationsNum=100;
%
设置最大迭代次数(即对于每个独立分量而言迭代均不超过此次数)
IterationsNum=0;
b=rand(numofIC,1)-.5;
%
随机设置b初值
b=b/norm(b);
% 对b标准化
while
i<=maxIterationsNum+1
if i == maxIterationsNum
% 循环结束处理
fprintf('\n第%d分量在%d次迭代内并不收敛。', r,maxIterationsNum);
break;
end
bOld=b;
a2=1;
u=1;
t=X'*b;
g=(exp(2.*t)-1)./(exp(2.*t)+1);
dg=4*exp(2.*t)./(exp(2.*t)+1).^2;
% g的一阶导数
b=((1-u)*t'*g*b+u*X*g)/SampleNum-mean(dg)*b;
% 核心公式,参见理论部分公式2.52
b=b-B*B'*b;
% 对b正交化
b=b/norm(b);
if abs(abs(b'*bOld)-1)<1e-9
% 如果收敛,则保存b
B(:,r)=b;
break;
end
i=i+1;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%% 数据复原并构图
%%%%%%%%%%%%%%%%%%%%%%%%
ICAedS=B'*Q*MixedS_bak;
% 参见理论部分公式2.55
% 将混合矩阵重新排列并输出
subplot(4,3,7),plot(ICAedS(1,:)),title('ICA输出信号1'),
subplot(4,3,8),plot(ICAedS(2,:)),title('ICA输出信号2'),
subplot(4,3,9),plot(ICAedS(3,:)),title('ICA输出信号3'),
%%%%%%%%%%%%%%%%%%%%%%%%% PCA计算并构图
%%%%%%%%%%%%%%%%%%%%%%%%
[V,D]=eig(MixedS_cov);
Vtmp=zeros(size(V,1),1);
for j=1:2
for i=1:2
if D(i,i)
tmp=D(i,i);Vtmp=V(:,i);
D(i,i)=D(i+1,i+1);V(:,i)=V(:,i+1);
D(i+1,i+1)=tmp;V(:,i+1)=Vtmp;
end
end
end
t1=(MixedS'*V(:,1))';
t2=(MixedS'*V(:,2))';
t3=(MixedS'*V(:,3))';
subplot(4,3,10),plot(t1),title('PCA输出信号1'),
subplot(4,3,11),plot(t2),title('PCA输出信号2'),
subplot(4,3,12),plot(t3),title('PCA输出信号3'),
![](http://simg.sinajs.cn/blog7style/images/common/sg_trans.gif)
algorithm" TITLE="ICA algorithm" />
因为一般数据需要去相关性,通常要白化处理,白化推导过程如下:
![](http://simg.sinajs.cn/blog7style/images/common/sg_trans.gif)
algorithm" TITLE="ICA algorithm" />
相关文章推荐
- ICA IBM Content A…
- GVF&nbsp;-&nbsp;a&nbsp;new&nbsp;snake&nbsp;algorithm&nbsp;for&nbsp;…
- think_java第二章&nbsp;一切皆对象(1)
- vc&nbsp;bmp对象与Opencv&nbsp;Iplimage对象…
- IDES Access Key破解 (转帖)
- Java IO流读写文件
- (转)Monkey test的使用
- SQL&nbsp;Block&nbsp;输出选项
- C语言&nbsp;结构体的内存对齐问题与位域
- OO ALV常用功能完整简例(热键单击…
- S3C2410&nbsp;Linux&nbsp;IIS音频设备驱动分…
- DB LUW 与 SAP LUW(一)
- C#&nbsp;集合类&nbsp;Array&nbsp;Arraylist&nbsp;List&nbsp;H…
- Android--获取电池信息(Battery in…
- linux&nbsp;下&nbsp;opencv2.0&nbsp;的编译与安装&nbsp;…
- Native SQL 整理
- JavaScript&nbsp;RegExp&nbsp;对象
- Deep learning:五(regulariz…
- 修复Ubuntu 10.10中的Alt+Printcre…
- Video&nbsp;for&nbsp;Linux&nbsp;Two&nbsp;API&nbsp;Specific…