Formalismo ed applicabilità

del metodo ICA

(Independent Component Analysis)

Francesca MarcucciFrancesca MarcucciUniversità di Perugia e INFNUniversità di Perugia e INFN

Udine 31 gennaio 2003Udine 31 gennaio 2003

ICA: Independent Component Analysis

i=1,..,Mj=1,..,N

Modello:

Supponiamo che MM segnali di media nulla ss11 s s2 2 ss33 .............. ssMM ma siano osservabili solo NN combinazioni lineari delle variabili si

xj = aij · si x=As A={aij}

Statistical “latent variables” model

Se la matrice AA non è nota il problema puo’ essere risolto facendo alcune assunzioni sulle proprietà statistiche delle sorgenti ssii

Dovrebbe essere considerato anche un termine aggiuntivo nn per il rumorex=As + n

ICA è una tecnica statistica per la decomposizione di un complesso datasetnelle sue sottoparti indipendenti ed è particolarmente utile nella soluzione

di problemi di Blind Source Separation (BSS)

Ipotesi per l’ applicabilità di ICA:

N >= MN >= M (nel seguito assumiamo M=NM=N senza perdita di generalità) al massimo una delle sorgenti e’ gaussiana le sorgenti ssii sono statisticamente indipendenti la matrice AA ha rango massimo per ora n=0n=0 , ma il modello puo’ essere esteso anche se di più difficile risoluzione

Ambiguità del metodo:Ambiguità del metodo:Il metodo fornisce una misura dell’ indipendenza delle componenti ma Il metodo fornisce una misura dell’ indipendenza delle componenti ma non da informazioni sull’Energia (non da informazioni sull’Energia (varianzavarianza) e sull’ ) e sull’ ordineordine in cui si in cui si ottengono, ovvero la matrice A puo’ essere scritta (dopo la ottengono, ovvero la matrice A puo’ essere scritta (dopo la convergenza) come:convergenza) come:

A=PD A=PD P P=permutazione =permutazione DD=matrice =matrice diagonalediagonale

Soluzione:Soluzione: Whitening or Sphering Whitening or Sphering

PREPROCESSAMENTO DEI DATI:PREPROCESSAMENTO DEI DATI:

Se i dati si non hanno media nulla allora si sottrae il valor medio (ad xi)

Whitening o sphering:Whitening o sphering:

Serve ad ottenere dei nuovi dati x con varianza unitaria

x=Vx dove E{ x xT }=1 (=I)

Se E{ x xT }=C allora V=C-1/2

infatti E{ x xT }= E{V x xT VT}=C-1/2 C C-1/2=I

Illustrazione del metodo:Illustrazione del metodo:Supponiamo di avere due variabili indipendenti uniformemente distribuite nella regione illustrata, con media nulla e varianza unitaria

1/23 |si|<3

Ad es. P(si)=

0 altrove

Applichiamo

2 3

A= 2 1

Le direzioni ci danno

informazione sulle colonne di A

s1

s2

x2

x1

Per stimare una delle componenti indipendenti consideriamo

y = wy = wT T x = x = ii w wii x xii

se w fosse l’ i-ima riga di A-1 allora y= si

z = AT w y = wT x = wT As = zT s

Se i dati hanno varianza unitaria WW e’ una matrice ortogonale WWWWTT=I=I

… cerchiamo un metodo piu’ generaleProcedimento:

Si basa sul teorema del “limite centrale”:

La distribuzione della somma di variabili random indipendenti tende ad una distribuzione gaussiana

Come usare il teorema del limite centrale?Come usare il teorema del limite centrale?

Ora abbiamo y = zTs ossia una combinazione lineare delle sorgenti indipendenti. Tale somma è “piu’ gaussiana” delle componenti originarie e lo diventa “al minimo” quando y=si ossia z ha solo l’i-imo elemento non nullo.

w scelto in modo da massimizzare la non-gaussianità di wTxMisure di non-gaussianità:Misure di non-gaussianità:

KURTOSIS KURTOSIS kurt(y)=E{y4} – 3(E{y2})2 è nullo per variabili gaussiane

quindi si cerca il max di |Kurt(y)|

NEGENTROPYNEGENTROPY J(y)=H(ygauss) – H(y) con H(y)= f(y) log f(y) dy

è nulla per variabili gaussiane (quelle con la max entropia H)

MUTUAL INFORMATIONMUTUAL INFORMATION I(y1,…,yM) = i H(yi) – H(y)

È nulla per variabili indipendenti e non negativa va minimizzata

Modello di rete neurale:Modello di rete neurale:

y

x x

W Q

y è una stima del vettore s y = W x

Q è una stima della matrice A x = Q y

1. “Apprende” una matrice W tale che y=Wx sono indipendenti

2. “Apprende” una matrice Q tale da minimizzare

E{||n’||2}=E{||x-Qy||2}

x

x

x

V Q

y

BT

W=BTVCon pre-withening:Con pre-withening:

LEARNING:LEARNING:

Massimizzare/minimizzare rispetto a w una delle funzioni F(w) precedenti imponendo dei vincoli ad esempio E{y2}=1 e E{y}=0 , ad esempio utilizzando i moltiplicatori di Lagrange :

gradient-ascendent method: wk+1 = wk + L ’wk

Newton-Like method:

L ’’w2 = r(w) Rxx L ’’w

2 wk= -L ’wk wk+1 = wk - Rxx-1 L ’wk/ r(w)

ALGORITMI:ALGORITMI:

Herault-Jutten: fallisce per piu’ di 2 sorgenti

EASI: performance uniforme

Bell’s and Seinowsky’s: performance uniforme e non richiede pre-withening

Chicocki and Amari: per feedforward e recurrent network

BIGRADIENT: necessario prewhitening, molto flessibile

NONLINEAR PCA: senza prewhitening separa solo componenti sinusoidali.

Adatto principalmente per funzioni sub-gaussiane

FastICA :FastICA :Caso semplice “one-unit”

(una sola unità computazionale (una sola unità computazionale 1 neurone con peso w) 1 neurone con peso w)

FastICA trova un vettore unitario w tale che massimizzi la non-gaussianità di wTx (utilizzando la Negentropy) con il metodo Newton-Like

1. Sceglie un iniziale vettore w random

2. calcola w+ = E{xg (wTx)} – E{g ’(wTx)} w g derivata di una funzione non quadratica

3. controlla se w = w+ / ||w+||

4. se non converge (w w+ = 1, hanno la stessa direzione) ritorna al punto 2

Tale algoritmo one-unit permette di determinare solo 1 componente Tale algoritmo one-unit permette di determinare solo 1 componente ma può essere facilmente esteso per la stima di più componenti ma può essere facilmente esteso per la stima di più componenti indipendenti improntando una rete “several-unit” con neuroni di indipendenti improntando una rete “several-unit” con neuroni di

pesi wpesi w11,…,w,…,wnn

Converge più rapidamente del metodo ICA; non necessita della stima di funzioni g o di parametri di altri parametri, è gratuito e disponibile sul web.

Recente applicazione:Recente applicazione:(Baccigalupi et al. 2002)

SCOPO:SCOPO: Separazione di componenti astrofisiche sovrapposte, ricostruendone sia le caratteristiche spaziali che spettrali, senza assunzioni a priori se non l’indipendenza e l’assenza di componenti gaussiane

MODELLO:MODELLO:

xi(r,)=ij sj(r,) (N differenti processi fisici) x=vettore M-dim

M canali di misura (diverse bande di frequenza) e strumento caratterizzato da una PSF B(r,) e funzione di risposta t(’)

x(r)= B(r-r’,’)j t(’)sj(r’,’) dr’ d’ + n

Ipotesi:Ipotesi:

1. funzione separabile sj(r,) = fj() sj(r)

2. B(r-r’,)=B(r) indipendente dalla frequenza

3. aij= t(’)f (’) d’ x(r)=A s(r) * B(r) + n

4. n è un rumore bianco, indipendente dal segnale , Gaussiano e stazionario

Synchrotron

angular

power

spectra

Input output

Input output

CMB

angular

power

spectra

Limiti:Limiti:

La ricostruzione della matrice di separazione peggiora nellì ipotesi in cui il rapporto tra due componenti è fortemente variabile lungo la

skymap

ES: le polveri dominano sul piano galattico mentre CMB domina ad alte latitudini

La ricostruzione è ottenuta con un’ errore migliore dell’ 1% nelle regioni in cui S/N 1.5 ,

l’errore cresce fino al 10% per S/N 1

Ancora un’ applicazione in astrofisica:Ancora un’ applicazione in astrofisica:

(Maria Funaro, Erkki Oja e Harri Valpola,2002)

Scopo: Individuare e rimuovere gli “artefacts” che influenzano le immagini (fluttuazioni,stelle della nostra galassia, rumore strumentale) basandosi sull’ analisi del profilo temporale della luminosità dei pixel e sull’indipendenza delle componenti dell’ immagine.

Dati: Immagini della Galassia M31

Modello: NN pixel T T immagini MM sorgenti X = ASX = AS

XX matrice TxN riga Xt: singola immagine al tempo t

colonna Xn: serie temporali (curve di luce) del pixel n

SS matrice MxN righe: immagini delle componenti indipendenti

per il singolo pixel n Xn = am Smn

AA matrice TxM : le M colonne di A (mixing vectors am) sono delle “curve di

luce virtuali” le cui combinazioni lineari danno quelle reali Xn

am caratterizza il comportamento temporale della sorgente m

Smn caratterizza il comportamento spaziale

T = 35 e N=100x100 pixel dopo whitening 3510 componenti indip.

Immagine originaria

1° e 2° autovettori: Raggi cosmici

5° autovettore: Sorgente puntiforme5° autovettore: Sorgente puntiforme

Conclusioni Ci sono pochi casi in letteratura di applicazioni ICA in

astrofisica, ma in questo campo l’indipendenza delle componenti assicura l’applicabilità del metodo.

La bontà statistica del metodo e’ legata principalmente alla minimizzazione della funzione di costo nella rete neurale implementata .

E’ necessario verificarne l’accuratezza con modelli simulati piu’ vicini alla realtà osservativa

ICA è sicuramente piu’ rapido dei metodi tradizionali … ma è ugualmente attendibile?

PROPOSTA:

Utilizzare in un primo momento FastICA ,Likelihood e Wavelet con gli output del light simulator e confrontare i due metodi