独立成分分析（Independent Component Analysis，ICA）原理及代码实现独立成分分析（Independent Component Analysis，ICA）MATLAB代码实现References

62 阅读 0 评论 41 点赞

我是靠谱客的博主任性鞋垫，最近开发中收集的这篇文章主要介绍独立成分分析（Independent Component Analysis，ICA）原理及代码实现独立成分分析（Independent Component Analysis，ICA）MATLAB代码实现References，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

过程监控中会用到很多中方法，如主成分分析（PCA）、慢特征分析（SFA）、概率MVA方法或独立成分分析（ICA）等为主流算法。

其中PCA主要多用于降维及特征提取，且只对正太分布（高斯分布）数据样本有效；SFA被用来学习过程监控的时间相关表示，SFA不仅可以通过监测稳态分布来检测与运行条件的偏差，还可以根据时间分布来识别过程的动态异常，多用于分类分析；概率MVA方法，多以解决动力学、时变、非线性等问题。

今天要介绍的是独立成分分析（ICA），由浅入深，细细道来。

此外文末还附有ICA可实现的代码哟~不要错过

独立成分分析（Independent Component Analysis，ICA）

基本原理

在信号处理中，独立成分分析（ICA）是一种用于将多元信号分离为加性子分量的计算方法。这是通过假设子分量是非高斯信号，并且在统计上彼此独立来完成的。ICA是盲源分离的特例。一个常见的示例应用程序是在嘈杂的房间中聆听一个人的语音的“ 鸡尾酒会问题 ”。

首先，引入一下经典的鸡尾酒宴会问题（Cocktail Party Problem）。

如下图所示，假如在我们的Party中有 $n$ 个人，他(她)们可以同时且互相说话，我们也在房间中一些角落里共放置了n个声音接收器（Microphone）用来记录声音。宴会过后，我们从n个麦克风中得到了一组数据 $mathrm{X}^{(i)}=left(x_{1}^{(i)}, x_{2}^{(i)}, ldots, x_{n}^{(i)}right) ;( i=1, ldots, m)$ ，其中 $i$ 表示采样的时间顺序，也就是说共得到了 $m$ 组采样，每一组采样都是 $n$ 维的，得到的是一个 $n * m$ 的数据矩阵。

我们的目标是：单独从这 $m$ 组采样数据中分辨出每个人说话的信号。

手残灵魂画师原图之鸡尾酒宴会

这里我们还需要另一组信息，那就是声音信号。我们把每个人的声音信号，那么也就是有 $n$ 个信号源 $S=(s_1,s_2,...,s_n)^T，sin mathfrak{R}^n$ ，且每个人发出的声音都是相互独立的。同样我们做了 $m$ 组采样。

假设我们令一个未知的混合系数矩阵（mixing coefficient matrix）为 $A$ ，用来组合叠加信号 $S$ ，
$X = A S$
这里的 $X$ 不是一个向量，是一个矩阵。其中的列向量是 $x^{(i)}=As^{(i)}$ 。

由此我们可知每一个 $x^{(i)}$ 都是由 $s^{(i)}$ 的分量线性表示的。

目前我们需要明确一点，只有 $X$ 是已知的可测量的（由麦克风得到），其中 $A 、 S$ 均为未知条件，我们要想办法根据 $X$ 来推出 $S$ 值。这个过程也称作为盲信号分离。

由于求矩阵的逆在实际运算中会出现一些问题，那么我们令 $H=A^{-1}$ ，则 $s^{(i)}=A^{-1}x^{(i)}=Hx^{(i)}$ 。

其中我们 $H$ 可表示为，
$h_{1}^{T} & \ & ...&\ & h_{n}^{T} & end{bmatrix}$

$h_iin mathfrak{R}^n$ ，将 $H$ 写成向量的形式。由此可得，
$s_{j}^{(i)}=h_{j}^{T}x^{(i)}$

ICA的不确定性（ICA ambiguities）

由于 $h$ 和 $s$ 都不确定，那么在没有先验知识的情况下，无法同时确定这两个相关参数。

比如上面的公式 $s = h x$ 。当 $h$ 扩大两倍时， $s$ 只需要同时扩大两倍即可，等式仍然满足，因此无法得到唯一的 $s$ 。同时如果将人的编号打乱，变成另外一个顺序，那么只需要调换A的列向量顺序即可，因此也无法单独确定 $s$ 。这两种情况称为原信号不确定。

还有一种ICA不适用的情况，那就是信号不能是高斯分布的。假设只有两个人发出的声音信号符合多值正态分布， $s \sim N (0, 1)$ 。

简而言之，不合适，不适用于没先验知识的情况。

ICA 算法

下面直接上ICA算法。

独立成分分析 ICA(Independent Component Correlation Algorithm)是一种函数，X为n维观测信号矢量，S为独立的m（m<=n)维未知源信号矢量，矩阵A被称为混合矩阵。ICA的目的就是寻找解混矩阵W（A的逆矩阵），然后对X进行线性变换，得到输出向量U。

这里使用最大似然估计来解释算法，我们假定每个 $s_i$ 有概率密度 $p_s$ ，那么给定时刻原信号的联合分布就是
$mathrm{p}(mathrm{s})=prod_{i=1}^{n} p_{s}left(s_{i}right)$
此公式代表一个假设前提：每个人发出的声音信号各自独立。

有了 $p (s)$ ，我们可以求得 $p (x)$
$mathrm{p}(mathrm{x})=mathrm{p}_{s}(H x)|mathrm{H}|=|mathrm{H}| prod_{i=1}^{n} p_{s}left(h_{i}{ }^{T} xright)$
左边是每个采样信号 $x$ 的概率，右边是每个原信号概率的乘积的 $∣ H ∣$ 倍。

若没有先验知识，我们无法求得 $H$ 和 $s$ 。

因此我们需要知道 $p_s(s_i)$ ，我们打算选取一个概率密度函数赋给 $s$ ，但是我们不能选取高斯分布的密度函数。在概率论里我们知道密度函数p(x)由累计分布函数（cdf）F(x)求导得到。F(x)要满足两个性质是：单调递增和在[0,1]。我们发现sigmoid函数很适合，定义域负无穷到正无穷，值域0到1，缓慢递增。我们假定 $s$ 的累积分布函数符合sigmoid函数
$g(s)=frac{1}{1+e^{-s}}$
求导可得，
$p_{s}(s)=g^{prime}(s)=frac{e^{s}}{left(1+e^{s}right)^{2}}$
这就是 $s$ 的密度函数。此时的 $s$ 是实数。

要是我们预先知道 $s$ 的分布函数，那就不用假设了，但在未知的情况下，sigmoid函数能够在大多数问题上取得不错的效果。

由于上式中 $p_s(s)$ 是个对称函数，因此E[s]=0（s的均值为0），那么E[x]=E[As]=0，x的均值也是0。

现在我们知道了 $p_s(s)$ ，下面开始求 $H$ 。

采样后的训练样本为 $mathrm{X}^{(i)}=left(x_{1}^{(i)}, x_{2}^{(i)}, ldots, x_{n}^{(i)}right) ;( i=1, ldots, m)$ ，使用前面得到的 $x$ 的概率密度函数，得其样本对数似然估计：
$ell(H)=sum_{i=1}^{m}left(sum_{j=1}^{n} log g^{prime}left(h_{j}^{T} x^{(i)}right)+log |H|right)$
其中，括号里的一大堆为 $p(x^{(i)})$ ，然后再对 $H$ 进行求导操作。在上式中包含有行列式，对行列式|W|进行求导的方法可参考这里。

最终得到的求导结果公式（~~很复杂很繁琐~~–心情）：
$gleft(h_{1}^{T} x^{(i)}right) \ 1-2 gleft(h_{2}^{T} x^{(i)}right) \ vdots \ 1-2 gleft(h_{n}^{T} x^{(i)}right) end{array}right] x^{(i)^{T}}+left(H^{T}right)^{-1}right)$
其中 $α$ 表示的是梯度上升速率，可自定义。

当通过多次迭代后，可求出 $H$ ，便可得到 $s^{(i)}=Hx^{(i)}$ 来还原出原始信号。

举个Paper的栗子

下面是 $s = 2$ 的原始信号：
在这里插入图片描述
下面为我们观测到的信号：

然后，再通过ICA还原后的信号为：

MATLAB代码实现

MATLAB代码：
Fast ICA

% Input：X 行变量维数，列采样个数；需要对原始矩阵转置
% Output：Sources重构的原信号, Q白化矩阵, P白化信号解混矩阵
function [Sources, Q, P] = FastICA(X, P)
% 白化处理
[dim, numSample] = size(X);
Xcov = cov(X');
[U, lambda] = eig(Xcov);
Q = lambda^(-1/2)*U';
Z = Q*X;
% FastICA
maxiteration = 10000; %最大迭代次数
error = 1e-5; % 收敛误差
% P = randn(dim,dim); % 随机初始化P，并按照列更新

for k = 1:dim
    Pk = P(:,k);
    Pk = Pk./norm(Pk); % 向量归一化
    lastPk = zeros(dim,1); % 0不需要再归一化
    count = 0;
    while abs(Pk - lastPk)&abs(Pk + lastPk) > error
        count = count + 1;
        lastPk = Pk;
        g = tanh(lastPk'*Z); % g(y)函数
        dg = 1 - g.^2; % g(y)的一阶导函数
%-------------------------------核心公式------------------------------------        
        Pk = mean(Z.*repmat(g,dim,1), 2) - repmat(mean(dg),dim,1).*lastPk;
        Pk = Pk - sum(repmat(Pk'*P(:,1:k-1),dim,1).*P(:,1:k-1),2);
        Pk = Pk./norm(Pk);
%--------------------------------------------------------------------------       
        if count == maxiteration
            fprintf('第%d个分量在%d次迭代内不收敛！n',k,maxiteration);
            break;
        end
    end
    P(:,k) = Pk;
end
    Sources = P'*Z;
% end

此外还有基于故障诊断的ICA算法代码实现->在这里

下面给出部分代码：

% 基于ICA的故障诊断
clear;clc;close all;
load('MPD2000.mat');

Xnormal = MPD0';
% 数据归一化
[dim, numSample] = size(Xnormal);
XnormalMean = mean(Xnormal, 2);
XnormalStd = std(Xnormal, 0, 2);
XnormalNorm = normalization(Xnormal, XnormalMean, XnormalStd);

% 正常数据计算 解混矩阵W
% P = rand(dim,dim)*100;
load('P.mat');
[S, Q, P] = FastICA(XnormalNorm, P);
W = P'*Q;
% ------------------------利用2范数大小对W的行重新排列---------------------
Wnorm = zeros(dim,1);
for k = 1:dim
    Wnorm(k) = norm(W(k,:));
end
[Wnorm, indices] = sort(Wnorm, 'descend');

% -------------------------确定主导成分Sd与参与成分Se----------------------
threshold = 0.80;
percentage = cumsum(Wnorm)./sum(Wnorm);
for k = 1:dim
    if percentage(k) > threshold
        break;
    end
end

效果如下：

References

https://www.cnblogs.com/jerrylead/archive/2011/04/19/2021071.html
https://baike.baidu.com/item/ICA/4972405?fr=aladdin
https://spaces.ac.cn/archives/2383
A. Hyva¨rinen, E. Oja*. Independent component analysis: algorithms and applications.