Светът наоколо е динамична смесица от сигнали от различни източници. Точно както цветовете в горната снимка преливат един в друг, пораждайки нови нюанси и тонове, всичко, което възприемаме, е сливане на по-прости компоненти. През повечето време дори не сме наясно, че светът около нас е такава хаотична смесица от независими процеси. Само в ситуации, в които различни стимули, които не се смесват добре, се конкурират за нашето внимание, ние осъзнаваме тази бъркотия. Типичен пример е сценарият на коктейл, където човек слуша гласа на друг човек, докато филтрира гласовете на всички останали гости. В зависимост от силата на звука в стаята това може да бъде проста или трудна задача, но по някакъв начин нашите мозъци са способни да отделят сигнала от шума. Въпреки че не е разбрано как мозъците ни правят това разделяне, има няколко изчислителни техники, които имат за цел да разделят сигнала на основните му компоненти. Един от тези методи се нарича Iнезависим Cкомпонентен Aанализ (ICA) и тук ще имаме по-близък вижте как работи този алгоритъм и как да го запишете в код на Python. Ако се интересувате повече от кода, отколкото от обяснението, можете директно да разгледате Бележника на Jupyter за тази публикация в Github.
Какво представлява анализът на независимите компоненти?
Нека засега останем с примера с коктейла. Представете си, че двама души говорят, можете да чуете и двамата, но единият е по-близо до вас от другия. Звуковите вълни от двата източника ще се смесят и ще достигнат до ушите ви като комбиниран сигнал. Вашият мозък ще размеси двата източника и ще възприемете гласовете и на двамата гости поотделно, като този, който стои по-близо до вас, е по-силен. Сега нека опишем това по по-абстрактен и опростен начин. Всеки източник е синусоида с постоянна честота. И двата източника се смесват в зависимост от това къде стоите. Това означава, че източникът, който е по-близо до вас, ще бъде по-доминиращ в смесения сигнал от този, който е по-далеч. Можем да запишем това по следния начин във векторно-матрица:
Където x е наблюдаваният сигнал, s са сигналите на източника и A е смесителната матрица. С други думи, нашият модел предполага, че сигналите x се генерират чрез линейна комбинация от изходните сигнали. В кода на Python нашият пример ще изглежда така:
>> import numpy as np
>>> # Number of samples
>>> ns = np.linspace(0, 200, 1000)
>>> # Sources with (1) sine wave, (2) saw tooth and (3) random noise
>>> S = np.array([np.sin(ns * 1),
signal.sawtooth(ns * 1.9),
np.random.random(len(ns))]).T
>>> # Quadratic mixing matrix
>>> A = np.array([[0.5, 1, 0.2],
[1, 0.5, 0.4],
[0.5, 0.8, 1]])
>>> # Mixed signal matrix
>>> X = S.dot(A).T
Както може да се види от графиките на Фигура 1 по-долу, кодът генерира един синусоидален сигнал, един трионен сигнал и някакъв случаен шум. Тези три сигнала са наши независими източници. В графиката по-долу можем да видим и трите линейни комбинации на сигналите на източника. По-нататък виждаме, че първият смесен сигнал е доминиран от компонента на триона, вторият смесен сигнал е повлиян повече от компонента на синусоида и последният смесен сигнал е доминиран от компонента на шума.
Сега, според нашия модел, можем да извлечем изходните сигнали отново от смесените сигнали, като умножим x с обратното на A:
Това означава, че за да намерим източника на сигнали, трябва да изчислим W. Така че задачата за останалата част от тази публикация ще бъде да се намери W и да се извлекат трите независими източника на сигнала от трите смесени сигнала.
Предпоставки за работа на ICA
Сега, преди да продължим, трябва да помислим малко повече за това какви свойства трябва да притежават нашите сигнали на източника, така че ICA да оцени успешно W. Първото предварително условие за работата на алгоритъма е, че смесените сигнали са линейна комбинация от произволен брой сигнали на източника. Второто предварително условие е сигналите на източника да са независими. И така, какво означава независимост? Два сигнала са независими, ако информацията в сигнал s1 не дава информация за сигнал s2. Това означава, че те не са корелирани, което означава, че тяхната ковариация е 0. Тук обаче трябва да се внимава, тъй като некорелацията не означава автоматично независимост. Третото предварително условие е независимите компоненти да не са гаусови. Защо така? Съвместното разпределение на плътността на два независими негаусови сигнала ще бъде равномерно на квадрат; вижте горния ляв график на Фигура 2 по-долу. Смесването на тези два сигнала с ортогонална матрица ще доведе до два сигнала, които сега вече не са независими и имат равномерно разпределение върху паралелограма; вижте долния ляв график на Фигура 2. Което означава, че ако сме на минималната или максималната стойност на един от нашите смесени сигнали, ние знаем стойността на другия сигнал. Следователно те вече не са независими. Правенето на същото с два гаусови сигнала ще доведе до нещо друго (вижте десния панел на Фигура 2). Съвместното разпределение на сигналите на източника е напълно симетрично, както и съвместното разпределение на смесените сигнали. Следователно той не съдържа информация за матрицата на смесване, чиято обратна стойност искаме да изчислим. От това следва, че в този случай ICA алгоритъмът ще се провали.
И така, в обобщение, за да работи алгоритъмът ICA, трябва да бъдат изпълнени следните предварителни условия: Нашите източници са (1) линейна смес от (2) независими, ( 3) негаусови сигнали.
Така че нека бързо проверим дали нашите тестови сигнали отгоре отговарят на тези предварителни условия. В лявата графика по-долу виждаме графиката на сигнала на синусоида спрямо сигнала на триона, докато цветът на всяка точка представлява компонента на шума. Сигналите се разпределят на квадрат, както се очаква за негаусови случайни променливи. По същия начин смесените сигнали образуват успоредник в десния график на Фигура 3, което показва, че смесените сигнали вече не са независими.
Стъпки на предварителна обработка
Сега приемането на смесените сигнали и подаването им директно в ICA не е добра идея. За да получите оптимална оценка на независимите компоненти, препоръчително е да направите известна предварителна обработка на данните. По-долу са обяснени по-подробно двете най-важни техники за предварителна обработка.
Центриране
Първата стъпка на предварителна обработка, която ще обсъдим тук, е центриране. Това е просто изваждане на средната стойност от нашия вход X. В резултат на това центрираните смесени сигнали ще имат нулева средна стойност, което означава, че и нашите сигнали на източника s са с нулева средна стойност. Това опростява изчислението на ICA и средната стойност може по-късно да бъде добавена обратно. Функцията за центриране в Python изглежда по следния начин.
>>> def center(x): >>> return x - np.mean(x, axis=1, keepdims=True)
Избелване
Втората стъпка на предварителна обработка, от която се нуждаем, е избелване на нашите сигнали X. Целта тук е да се трансформира линейно X, така че потенциалните корелации между сигналите да бъдат премахнати и техните дисперсии да са равни на единица. В резултат ковариационната матрица на избелените сигнали ще бъде равна на идентичната матрица:
Където I е матрицата на идентичността. Тъй като също трябва да изчислим ковариацията по време на процедурата за избелване, ще напишем малка функция на Python за нея.
>>> def covariance(x): >>> mean = np.mean(x, axis=1, keepdims=True) >>> n = np.shape(x)[1] - 1 >>> m = x - mean >>> return (m.dot(m.T))/n
Кодът за стъпката на избелване е показан по-долу. Базира се на разлагането на сингулярна стойност (SVD) на ковариационната матрица на X. Ако се интересувате от подробности за тази процедура, препоръчвам тази статия.
>>> def whiten(x):
>>> # Calculate the covariance matrix
>>> coVarM = covariance(X)
>>> # Singular value decoposition
>>> U, S, V = np.linalg.svd(coVarM)
>>> # Calculate diagonal matrix of eigenvalues
>>> d = np.diag(1.0 / np.sqrt(S))
>>> # Calculate whitening matrix
>>> whiteM = np.dot(U, np.dot(d, U.T))
>>> # Project onto whitening matrix
>>> Xw = np.dot(whiteM, X)
>>> return Xw, whiteM
Внедряване на алгоритъма FastICA
Добре, сега, когато разполагаме с нашите функции за предварителна обработка, най-накрая можем да започнем да прилагаме алгоритъма ICA. Има няколко начина за прилагане на ICA въз основа на контрастната функция, която измерва независимостта. Тук ще използваме приближение на негентропиявъв версия на ICA, наречена FastICA.
И така, как работи? Както беше обсъдено по-горе, едно от предпоставките за работа на ICA е сигналите на нашия източник да не са гаусови. Интересно нещо за два независими, не-гаусови сигнала е, че тяхната сума е по-гаусова от който и да е от сигналите на източника. Следователно трябва да оптимизираме W по начин, че получените сигнали на Wx да са възможно най-различни от Гаус. За да направим това, се нуждаем от мярка за гауситност. Най-простата мярка би била ексцесът, който е четвъртият момент от данните и измерва „опашатостта“ на разпределението. Нормалното разпределение има стойност 3, равномерното разпределение като това, което използвахме на Фигура 2, има ексцес ‹ 3. Имплементацията в Python е права, както може да се види от кода по-долу, който също изчислява другите моменти от данните. Първият момент е средната стойност, вторият е дисперсията, третият е асимметрията и четвъртият е ексцесът. Тук 3 се изважда от четвъртия момент, така че нормалното разпределение има ексцес 0.
>>> def kurtosis(x): >>> n = np.shape(x)[0] >>> mean = np.sum((x**1)/n) # Calculate the mean >>> var = np.sum((x-mean)**2)/n # Calculate the variance >>> skew = np.sum((x-mean)**3)/n # Calculate the skewness >>> kurt = np.sum((x-mean)**4)/n # Calculate the kurtosis >>> kurt = kurt/(var**2)-3 >>> return kurt, skew, var, mean
За нашето внедряване на ICA обаче няма да използваме ексцес като контрастна функция, но можем да го използваме по-късно, за да проверим нашите резултати. Вместо това ще използваме следната контрастна функция g(u) и нейната първа производна g’(u):
Алгоритъмът FastICA използва горните две функции по следния начин витерационна схема с фиксирана точка:
Така че според горното това, което трябва да направим, е да направим случайно предположение за теглата на всеки компонент. Точковият продукт на произволните тегла и смесените сигнали се предава в двете функции g и g’. След това изваждаме резултата от g’ от g и изчисляваме средната стойност. Резултатът е нашият нов вектор на теглата. След това можем директно да разделим новия вектор на теглата на неговата норма и да повтаряме горното, докато теглата престанат да се променят. Няма да има нищо лошо в това. Въпреки това проблемът, пред който сме изправени тук, е, че в итерацията за втория компонент можем да идентифицираме същия компонент като в първата итерация. За да разрешим този проблем, трябва да декорираме новите тегла от предварително идентифицираните тегла. Това се случва в стъпката между актуализирането на теглата и разделянето на тяхната норма. Тогава в Python имплементацията изглежда по следния начин:
>>> def fastIca(signals, alpha = 1, thresh=1e-8, iterations=5000):
>>> m, n = signals.shape
>>> # Initialize random weights
>>> W = np.random.rand(m, m)
>>> for c in range(m):
>>> w = W[c, :].copy().reshape(m, 1)
>>> w = w/ np.sqrt((w ** 2).sum())
>>> i = 0
>>> lim = 100
>>> while ((lim > thresh) & (i < iterations)):
>>> # Dot product of weight and signal
>>> ws = np.dot(w.T, signals)
>>> # Pass w*s into contrast function g
>>> wg = np.tanh(ws * alpha).T
>>> # Pass w*s into g'
>>> wg_ = (1 - np.square(np.tanh(ws))) * alpha
>>> # Update weights
wNew = (signals * wg.T).mean(axis=1) -
>>> wg_.mean() * w.squeeze()
>>> # Decorrelate weights
>>> wNew = wNew -
np.dot(np.dot(wNew, W[:c].T), W[:c])
>>> wNew = wNew / np.sqrt((wNew ** 2).sum())
>>> # Calculate limit condition
>>> lim = np.abs(np.abs((wNew * w).sum()) - 1)
>>> # Update weights
>>> w = wNew
>>> # Update counter
>>> i += 1
>>> W[c, :] = w.T
>>> return W
Така че сега, след като имаме написан целия код, нека стартираме всичко!
>>> # Center signals >>> Xc, meanX = center(X) >>> # Whiten mixed signals >>> Xw, whiteM = whiten(Xc) >>> # Run the ICA to estimate W >>> W = fastIca(Xw, alpha=1) >>> #Un-mix signals using W >>> unMixed = Xw.T.dot(W.T) >>> # Subtract mean from the unmixed signals >>> unMixed = (unMixed.T - meanX).T
Резултатите от ICA са показани на Фигура 4 по-долу, където горният панел представлява оригиналните сигнали на източника, а долният панел - независимите компоненти, извлечени от нашата реализация на ICA. И резултатът изглежда много добър. Върнахме и трите източника!
Така че най-накрая нека проверим едно последно нещо: ексцесът на сигналите. Както можем да видим на Фигура 5, всички наши смесени сигнали имат ексцес ≤ 1, докато всички възстановени независими компоненти имат ексцес 1,5, което означава, че са по-малко гаусови от техните източници. Това трябва да е така, тъй като ICA се опитва да максимизира негаусианството. Също така добре илюстрира факта, споменат по-горе, че сместа от не-гаусови сигнали ще бъде по-гаусова от източниците.
И така, за да обобщим: видяхме как работи ICA и как да го внедрим от нулата в Python. Разбира се, има много налични реализации на Python, които могат да се използват директно. Въпреки това винаги е препоръчително да разберете основния принцип на метода, за да знаете кога и как да го използвате. Ако се интересувате да се потопите по-дълбоко в ICA и да научите за подробностите, препоръчвам този документ от Aapo Hyvärinen и Erkki Oja, 2000 г..
В противен случай можете да проверите пълния „код тук“, да ме последвате в „Twitter“ или да се свържете чрез „LinkedIn“.
Кодът за този проект може да бъде намерен в Github.
