Анализ независимых компонент

Данные представлены наблюдаемым случайным вектором ${\displaystyle {\boldsymbol {x}}=(x_{1},\ldots ,x_{m})^{T}}$, а скрытые компоненты случайным вектором ${\displaystyle {\boldsymbol {s}}=(s_{1},\ldots ,s_{n})^{T}}$. Задачей построения алгоритма является преобразование наблюдаемых данных ${\displaystyle {\boldsymbol {x}}}$ с помощью статического преобразования ${\displaystyle {\boldsymbol {W}}}$ в наблюдаемый вектор максимально независимых компонент ${\displaystyle {\boldsymbol {s}}={\boldsymbol {W}}{\boldsymbol {x}}}$, измеренных некоторой функцией независимости ${\displaystyle F(s_{1},\ldots ,s_{n})}$.

Компоненты ${\displaystyle x_{i}}$ наблюдаемого случайного вектора ${\displaystyle {\boldsymbol {x}}=(x_{1},\ldots ,x_{m})^{T}}$ генерируются как сумма независимых компонент ${\displaystyle s_{k}}$, ${\displaystyle k=1,\ldots ,n}$:

${\displaystyle x_{i}=a_{i,1}s_{1}+\cdots +a_{i,k}s_{k}+\cdots +a_{i,n}s_{n}}$

взвешенных весами ${\displaystyle a_{i,k}}$.

Та же генерирующая модель может быть записана в векторном виде как ${\displaystyle {\boldsymbol {x}}=\sum _{k=1}^{n}{\boldsymbol {s}}_{k}{\boldsymbol {a}}_{k}}$, где наблюдаемый случайный вектор ${\displaystyle x}$ представлен базисными векторами ${\displaystyle {\boldsymbol {a}}_{k}=({\boldsymbol {a}}_{1,k},\ldots ,{\boldsymbol {a}}_{m,k})^{T}}$. Базисные вектора ${\displaystyle {\boldsymbol {a}}_{k}}$ образуют столбцы матрицы смешивания ${\displaystyle {\boldsymbol {A}}=({\boldsymbol {a}}_{1},\ldots ,{\boldsymbol {a}}_{n})}$ и генерирующая формула может быть записана как ${\displaystyle {\boldsymbol {x}}={\boldsymbol {A}}{\boldsymbol {s}}}$, где ${\displaystyle {\boldsymbol {s}}=(s_{1},\ldots ,s_{n})^{T}}$.

Если дана модель и реализации ${\displaystyle x_{1},\ldots ,x_{N}}$ случайного вектора ${\displaystyle {\boldsymbol {x}}}$, задачей является оценка как матрицы смешивания ${\displaystyle {\boldsymbol {A}}}$, так и источников ${\displaystyle {\boldsymbol {s}}}$. Это делается путём адаптивного вычисления векторов ${\displaystyle {\boldsymbol {w}}}$ и установления функции цены, которая либо максимизирует негауссовость вычисленного ${\displaystyle s_{k}={\boldsymbol {w}}^{T}{\boldsymbol {x}}}$ или минимизирует взаимную информацию. В некоторых случаях априорное знание распределения вероятности источников может быть использовано в функции цены.

Исходные источники ${\displaystyle {\boldsymbol {s}}}$ могут быть извлечены путём умножения наблюдаемых сигналов ${\displaystyle {\boldsymbol {x}}}$ на обратную к матрице смешивания ${\displaystyle {\boldsymbol {W}}={\boldsymbol {A}}^{-1}}$, которая известна также как не смешивающая матрица. Здесь предполагается, что матрица смешивания квадратная (${\displaystyle n=m}$). Если число базисных векторов больше размерности наблюдаемых векторов ${\displaystyle n>m}$, задача является переопределённой, но остаётся разрешимой с помощью псевдообратной матрицы.

С добавочным предположением о нулевом среднем и не коррелирующим гауссовым шумом ${\displaystyle n\sim N(0,\operatorname {diag} (\Sigma ))}$, модель АНК принимает форму ${\displaystyle {\boldsymbol {x}}={\boldsymbol {A}}{\boldsymbol {s}}+n}$.

Смесь источников не обязательно должна быть линейной. Используя нелинейную функцию смешивания ${\displaystyle f(\cdot |\theta )}$ с параметрами ${\displaystyle \theta }$ нелинейной моделью АНК будет ${\displaystyle x=f(s|\theta )+n}$.

Независимые компоненты различимы с точностью до перестановки и масштабирования источников. Эта различимость требует, чтобы:

• Максимум один из источников ${\displaystyle s_{k}}$ был гауссовым,
• Число наблюдаемых смесей ${\displaystyle m}$ должно быть не меньше числа компонент ${\displaystyle n}$: ${\displaystyle m\geqslant n}$. Это эквивалентно высказыванию, что матрица смеси ${\displaystyle {\boldsymbol {A}}}$ должна иметь полный ранг, чтобы существовала обратная ей смесь.

Смеси сигналов имеют тенденцию к получению гауссовой плотности вероятности, а сигналы источников имеют тенденцию к негауссовой плотности вероятности. Каждый источник сигнала может быть выделен из набора смесей сигналов путём вычисления скалярного произведения вектора весов и той смеси сигналов, на которой это скалярное произведение даёт ортогональную проекцию смеси сигналов. Следующая задача заключается в нахождении вектора весов. Один из методов — поиск наилучшей проекции[2][7].

Поиск наилучшей проекции ищет одну проекцию за шаг, при условии, что выделенный сигнал будет настолько негауссовым, насколько это возможно. Это контрастирует с АНК, который обычно выделяет M сигналов одновременно из M смесей сигналов, что требует оценки ${\displaystyle M\times M}$ несмешивающей матрицы. Одним из практических преимуществ поиска наилучшей проекции над АНК является то, что может выделяться менее M сигналов, если требуется, где каждый источник сигнала выделяется из смеси M сигналов используя M-элементный вектор весов.

Мы можем использовать коэффициент эксцесса для извлечения сигнала с несколькими источниками путём нахождения правильных векторов весов с использованием поиска наилучшей проекции.

Коэффициент эксцесса плотности вероятности сигнала, для конечной выборки вычисляется как

${\displaystyle K={\frac {\operatorname {E} [(\mathbf {y} -\mathbf {\overline {y}} )^{4}]}{(\operatorname {E} [(\mathbf {y} -\mathbf {\overline {y}} )^{2}])^{2}}}-3}$

где ${\displaystyle \mathbf {\overline {y}} }$ является выборочным средним ${\displaystyle \mathbf {y} }$ выделенных сигналов. Константа 3 обеспечивает, чтобы гауссовы сигналы имели нулевой коэффициент эксцесса, супергауссовы сигналы имели положительный коэффициент эксцесса, а субгауссовы сигналы имели отрицательный коэффициент эксцесса. Знаменатель равен дисперсии ${\displaystyle \mathbf {y} }$ и он обеспечивает, чтобы измеренный коэффициент эксцесса получал дисперсию сигнала. Целью поиска наилучшей проекции является максимизация коэффициента эксцесса и сделать выделенный сигнал настолько ненормальным, насколько возможно.

Используя коэффициент эксцесса как меру ненормальности мы можем теперь проверить насколько коэффициент эксцесса сигнала ${\displaystyle \mathbf {y} =\mathbf {w} ^{T}\mathbf {x} }$, извлечённого из набора M смесей ${\displaystyle \mathbf {x} =(x_{1},x_{2},\ldots ,x_{M})^{T}}$, изменяется по мере того, как вектор весов ${\displaystyle \mathbf {w} }$ вращается вокруг начала координат. Если задано, что каждый источник сигнала ${\displaystyle \mathbf {s} }$ является супергауссовым, мы можем ожидать

1. коэффициент эксцесса извлечённого сигнала ${\displaystyle \mathbf {y} }$ максимален в точности тогда, когда ${\displaystyle \mathbf {y} =\mathbf {s} }$.
2. коэффициент эксцесса извлечённого сигнала ${\displaystyle \mathbf {y} }$ максимален, когда ${\displaystyle \mathbf {w} }$ ортогонален проекциям осей ${\displaystyle S_{1}}$ или ${\displaystyle S_{2}}$, поскольку мы знаем, что вектор оптимального веса должен быть ортогонален преобразованным осям ${\displaystyle S_{1}}$ и ${\displaystyle S_{2}}$.

Для смеси сигналов от разных источников мы можем использовать коэффициент эксцесса ортогонализации Грама ― Шмидта (ОГШ) для извлечения сигналов. Если дана смесь M сигналов в M-мерном пространстве, ОГШ проектирует эти точки данных в (M-1)-мерное пространство с помощью вектора весов. Мы можем гарантировать независимость выделенных сигналов с помощью ОГШ.

С целью поиска правильного значения ${\displaystyle \mathbf {w} }$ мы можем использовать метод градиентного спуска. Прежде всего, мы избавляемся от корреляции и преобразуем ${\displaystyle \mathbf {x} }$ в новую смесь ${\displaystyle \mathbf {z} }$, которая имеет единичную дисперсию и ${\displaystyle \mathbf {z} =(z_{1},z_{2},\ldots ,z_{M})^{T}}$. Этот процесс может быть выполнен путём применения сингулярного разложения к ${\displaystyle \mathbf {x} }$,

${\displaystyle \mathbf {x} =\mathbf {U} \mathbf {D} \mathbf {V} ^{T}}$

Масштабируем каждый вектор ${\displaystyle U_{i}=U_{i}/\operatorname {E} (U_{i}^{2})}$ и положим ${\displaystyle \mathbf {z} =\mathbf {U} }$. Сигнал, выделенный взвешенным вектором ${\displaystyle \mathbf {w} }$, равен ${\displaystyle \mathbf {y} =\mathbf {w} ^{T}\mathbf {z} }$. Если вектор весов w имеет единичную длину, то есть ${\displaystyle \operatorname {E} [(\mathbf {w} ^{T}\mathbf {z} )^{2}]=1}$, тогда коэффициент эксцесса можно переписать как:

${\displaystyle K={\frac {\operatorname {E} [\mathbf {y} ^{4}]}{(\operatorname {E} [\mathbf {y} ^{2}])^{2}}}-3=\operatorname {E} [(\mathbf {w} ^{T}\mathbf {z} )^{4}]-3.}$

Процесс обновления для ${\displaystyle \mathbf {w} }$:

${\displaystyle \mathbf {w} _{new}=\mathbf {w} _{old}-\eta \operatorname {E} [\mathbf {z} (\mathbf {w} _{old}^{T}\mathbf {z} )^{3}].}$

где ${\displaystyle \eta }$ является малой константой для гарантирования, что ${\displaystyle \mathbf {w} }$ сходится к оптимальному решению. После каждого обновления мы нормализуем ${\displaystyle \mathbf {w} _{new}={\frac {\mathbf {w} _{new}}{|\mathbf {w} _{new}|}}}$ и множество ${\displaystyle \mathbf {w} _{old}=\mathbf {w} _{new}}$ и повторяем процесс обновления пока он не сойдётся. Мы можем использовать также другой алгоритм для обновления вектора весов ${\displaystyle \mathbf {w} }$.

Другим подходом является использование негэнтропии[8] вместо коэффициента эксцесса. Негэнтропия является устойчивым методом по отношению коэффициента эксцесса, поскольку коэффициент эксцесса очень чувствителен к выбросам. Метод негэнтропии основывается на важном свойстве распределения Гаусса — нормальная случайная величина имеет наибольшую энтропию среди всех непрерывных случайных переменных с одинаковой дисперсией. Это также является причиной, почему мы хотим найти наиболее негауссовые переменные. Простое доказательство можно найти в статье дифференциальной энтропии.

${\displaystyle J(x)=S(y)-S(x)\,}$

y являются гауссовой случайной переменной некоторой ковариантной матрицы,

${\displaystyle S(x)=-\int p_{x}(u)\log p_{x}(u)du}$

Аппроксимация для негэнтропии равна

${\displaystyle J(x)={\frac {1}{12}}(E(x^{3}))^{2}+{\frac {1}{48}}(kurt(x))^{2}}$

Доказательство можно найти на странице 131 книги «Анализ независимых компонент», которую написали Аапо Хювяринен, Юха Кархунен и Эркки Ойя[3]. Эта аппроксимация также страдает теми же проблемами, что и коэффициент эксцесса (чувствительность к выбросам). Разрабатывались и другие подходы[9]

${\displaystyle J(y)=k_{1}(E(G_{1}(y)))^{2}+k_{2}(E(G_{2}(y))-E(G_{2}(v))^{2}}$

Выбор ${\displaystyle G_{1}}$ и ${\displaystyle G_{2}}$

${\displaystyle G_{1}={\frac {1}{a_{1}}}\log(\cosh(a_{1}u))}$ и ${\displaystyle G_{2}=-\exp(-{\frac {u^{2}}{2}})}$

АНК, по существу, является многомерной параллельной версией поиска наилучшей проекции. В то время как поиск наилучшей проекции выделяет серию сигналов по одному из смеси M сигналов, АНК выделяет M сигналов параллельно. Это приводит к большей устойчивости АНК по сравнению с поиском наилучшей проекции[2].

Метод поиска наилучшей проекции, чтобы обеспечить независимость выделяемых сигналов, использует ортогонализацию Грама ― Шмидта, в то время как АНК использует метод infomax и оценку максимального правдоподобия для обеспечения независимости выделяемого сигнала. Ненормальность выделяемого сигнала достигается с помощью соответствующей модели.

Процесс АНК, основанный на infomax, коротко: если дана смесь сигналов ${\displaystyle \mathbf {x} }$ и набор одинаковых независимых функций распределения ${\displaystyle g}$, мы ищем несмешивающую матрицу ${\displaystyle \mathbf {W} }$, которая максимизирует совместную энтропию сигналов ${\displaystyle \mathbf {Y} =g(\mathbf {y} )}$, где ${\displaystyle \mathbf {y} =\mathbf {Wx} }$ являются сигналами, отобранными по ${\displaystyle \mathbf {W} }$. Если дана оптимальная ${\displaystyle \mathbf {W} }$, сигналы ${\displaystyle \mathbf {Y} }$ имеют максимальную энтропию и, потому, независимы, что гарантирует, что выделенные сигналы ${\displaystyle \mathbf {y} =g^{-1}(\mathbf {Y} )}$ также независимы. Функция ${\displaystyle g}$ обратима и является моделью сигнала. Заметим, что если плотность вероятности модели источника сигнала ${\displaystyle p_{s}}$ соответствует плотности вероятности выделенного сигнала ${\displaystyle p_{\mathbf {y} }}$, то максимизация совместной энтропии ${\displaystyle Y}$ также максимизирует количество взаимной информации между ${\displaystyle \mathbf {x} }$ и ${\displaystyle \mathbf {Y} }$. По этой причине, использование энтропии для выделения независимых сигналов известно как infomax.

Рассмотрим энтропию векторной переменной ${\displaystyle \mathbf {Y} =g(\mathbf {y} )}$, где ${\displaystyle \mathbf {y} =\mathbf {Wx} }$ является набором сигналов, выделенных несмешивающей матрицей ${\displaystyle \mathbf {W} }$. Для конечного набора значений, выбранных из распределения с плотностью вероятности ${\displaystyle p_{\mathbf {y} }}$, энтропия ${\displaystyle \mathbf {Y} }$ может быть оценена как:

${\displaystyle H(\mathbf {Y} )=-{\frac {1}{N}}\sum _{t=1}^{N}\ln p_{\mathbf {Y} }(\mathbf {Y} ^{t})}$

Совместная плотность вероятности ${\displaystyle p_{\mathbf {Y} }}$, как можно показать, связана с совместной плотностью вероятности ${\displaystyle p_{\mathbf {y} }}$ извлечённых сигналов с помощью многомерной формы:

${\displaystyle p_{\mathbf {Y} }(Y)={\frac {p_{\mathbf {y} }(\mathbf {y} )}{|{\frac {\partial \mathbf {Y} }{\partial \mathbf {y} }}|}}}$

где ${\displaystyle \mathbf {J} ={\frac {\partial \mathbf {Y} }{\partial \mathbf {y} }}}$ является матрицей Якоби. Мы имеем ${\displaystyle |\mathbf {J} |=g'(\mathbf {y} )}$, и ${\displaystyle g'}$ является плотностью вероятности, принятых для источников сигналов ${\displaystyle g'=p_{s}}$, поэтому,

${\displaystyle p_{\mathbf {Y} }(Y)={\frac {p_{\mathbf {y} }(\mathbf {y} )}{|{\frac {\partial \mathbf {Y} }{\partial \mathbf {y} }}|}}={\frac {p_{\mathbf {y} }(\mathbf {y} )}{p_{\mathbf {s} }(\mathbf {y} )}}}$

поэтому,

${\displaystyle H(\mathbf {Y} )=-{\frac {1}{N}}\sum _{t=1}^{N}\ln {\frac {p_{\mathbf {y} }(\mathbf {y} )}{p_{\mathbf {s} }(\mathbf {y} )}}}$

Мы знаем, что когда ${\displaystyle p_{\mathbf {y} }=p_{s}}$, ${\displaystyle p_{\mathbf {Y} }}$ является однородным распределением, а ${\displaystyle H({\mathbf {Y} })}$ максимизирована. Поскольку

${\displaystyle p_{\mathbf {y} }(\mathbf {y} )={\frac {p_{\mathbf {x} }(\mathbf {x} )}{|{\frac {\partial \mathbf {y} }{\partial \mathbf {x} }}|}}={\frac {p_{\mathbf {x} }(\mathbf {x} )}{|\mathbf {W} |}}}$

где ${\displaystyle |\mathbf {W} |}$ является абсолютным значением определителя несмешивающей матрицы ${\displaystyle \mathbf {W} }$. Поэтому,

${\displaystyle H(\mathbf {Y} )=-{\frac {1}{N}}\sum _{t=1}^{N}\ln {\frac {p_{\mathbf {x} }(\mathbf {x} ^{t})}{|\mathbf {W} |p_{\mathbf {s} }(\mathbf {y} ^{t})}}}$

так что,

${\displaystyle H(\mathbf {Y} )={\frac {1}{N}}\sum _{t=1}^{N}\ln p_{\mathbf {s} }(\mathbf {y} ^{t})+\ln |\mathbf {W} |+H(\mathbf {x} )}$

поскольку ${\displaystyle H(\mathbf {x} )=-{\frac {1}{N}}\sum _{t=1}^{N}\ln p_{\mathbf {x} }(\mathbf {x} ^{t})}$, и максимизация ${\displaystyle \mathbf {W} }$ не влияет на ${\displaystyle H_{\mathbf {x} }}$, мы можем максимизировать функцию

${\displaystyle h(\mathbf {Y} )={\frac {1}{N}}\sum _{t=1}^{N}\ln p_{\mathbf {s} }(\mathbf {y} ^{t})+\ln |\mathbf {W} |}$

чтобы получить независимость извлечённого сигнала.

Если имеется M маргинальных плотностей вероятности модели, совместные плотности вероятности ${\displaystyle p_{\mathbf {s} }}$ независимы и используют супергауссову модель плотности вероятности для источников сигналов ${\displaystyle p_{\mathbf {s} }=(1-\tanh(\mathbf {s} )^{2})}$, то мы получаем

${\displaystyle h(\mathbf {Y} )={\frac {1}{N}}\sum _{i=1}^{M}\sum _{t=1}^{N}\ln(1-\tanh(\mathbf {w_{i}^{T}x^{t}} )^{2})+\ln |\mathbf {W} |}$

В сумме, если задана наблюдаемая смесь сигнала ${\displaystyle \mathbf {x} }$, соответствующий набор извлечённых сигналов ${\displaystyle \mathbf {y} }$ и модель источника сигнала ${\displaystyle p_{\mathbf {s} }=g'}$, мы можем найти оптимальную несмешивающую матрицу ${\displaystyle \mathbf {W} }$ и сделать извлечённые сигналы независимыми и негауссовыми. Подобно ситуации с поиском наилучшей проекции, мы можем использовать метод градиентного спуска для поиска оптимального решения несмешивающей матрицы.

Оценка максимального правдоподобия (англ. Maximum likelihood estimation, MLE) является стандартным статистическим средством для нахождения значений параметров (например, несмешивающей матрицы ${\displaystyle \mathbf {W} }$), которые обеспечивают лучшее соответствие некоторых данных (например, извлечённых сигналов ${\displaystyle y}$) для данной модели (например, совместной плотности вероятности (ПВ) ${\displaystyle p_{s}}$ источников сигналов)[2].

Модель максимального правдоподобия включает спецификацию плотности вероятности, которая в этом случае является плотностью вероятности ${\displaystyle p_{s}}$ сигналов неизвестного источника ${\displaystyle s}$. При использовании максимального правдоподобия целью является нахождение несмешивающей матрицы, которая даёт извлечённые сигналы ${\displaystyle y=\mathbf {W} x}$ с совместной плотностью вероятности, которые максимально подобны совместной плотностью вероятности ${\displaystyle p_{s}}$ сигналов неизвестного источника ${\displaystyle s}$.

Оценка максимального правдоподобия основывается на предположении, что если модель плотности вероятности ${\displaystyle p_{s}}$ и модель параметров ${\displaystyle \mathbf {A} }$ правильны, то должна быть получена высокая вероятность для ${\displaystyle x}$, что эти данные действительно наблюдаемы. Обратно, если ${\displaystyle \mathbf {A} }$ далёк от верных значений параметров, то следует ожидать низкую вероятность наблюдения данных.

При оценке максимального правдоподобия мы называем вероятность наблюдаемых данных для данного набора значений параметров модели (например, плотности вероятности ${\displaystyle p_{s}}$ и матрицы ${\displaystyle \mathbf {A} }$) правдоподобностью значений параметров модели, заданной наблюдаемыми данными.

Мы определяем функцию правдоподобия ${\displaystyle \mathbf {L(W)} }$ матрицы ${\displaystyle \mathbf {W} }$:

${\displaystyle \mathbf {L(W)} =p_{s}(\mathbf {W} x)|\det \mathbf {W} |.}$

Это равно плотности вероятности в ${\displaystyle x}$, поскольку ${\displaystyle s=\mathbf {W} x}$.

Тогда, если мы хотим найти ${\displaystyle \mathbf {W} }$, то наиболее вероятно иметь сгенерированные наблюдаемые смеси ${\displaystyle x}$ из неизвестных источников сигналов ${\displaystyle s}$ с плотностью вероятности ${\displaystyle p_{s}}$, то нам нужно лишь найти ${\displaystyle \mathbf {W} }$, которая максимизирует правдоподобность ${\displaystyle \mathbf {L(W)} }$. Несмешивающая матрица, которая максимизирует равенство, известна как оценка максимального правдоподобия оптимальной несмешивающей матрицей.

Распространённой практикой является использование логарифма правдоподобия, поскольку его проще всего вычислить. Так как логарифм является монотонной функцией, матрица ${\displaystyle \mathbf {W} }$, которая максимизирует функция ${\displaystyle \mathbf {L(W)} }$, также максимизирует его логарифм ${\displaystyle \ln \mathbf {L(W)} }$. Это позволяет взять логарифм в равенстве выше, что даёт логарифм функции правдоподобия

${\displaystyle \ln \mathbf {L(W)} =\sum _{i}\sum _{t}\ln p_{s}(w_{i}^{T}x_{t})+N\ln |\det \mathbf {W} |}$

Если мы подставим широко используемую модель плотности вероятности с высоким коэффициентом эксцесса для источников сигналов ${\displaystyle p_{s}=(1-\tanh(s)^{2})}$, мы получим

${\displaystyle \ln \mathbf {L(W)} ={1 \over N}\sum _{i}^{M}\sum _{t}^{N}\ln(1-\tanh(w_{i}^{T}x_{t})^{2})+\ln |\det \mathbf {W} |}$

Матрица ${\displaystyle \mathbf {W} }$, максимизирующая эту функцию, является оценкой максимального правдоподобия.