Formalismo ed applicabilità del metodo ICA (Independent Component Analysis)

1. Formalismo ed applicabilità del metodo ICA (Independent Component Analysis) Francesca Marcucci Università di Perugia e INFN Udine 31 gennaio 2003

2. ICA: Independent Component Analysis ICA è una tecnica statistica per la decomposizione di un complesso dataset nelle sue sottoparti indipendenti ed è particolarmente utile nella soluzione di problemi di Blind Source Separation (BSS) • Modello: • Supponiamo che  M segnali di media nullas1 s2 s3.......sMma siano osservabili solo • N combinazioni lineari delle variabili si • xj= aij · si x=AsA={aij} • Statistical “latent variables” model Se la matrice A non è nota il problema puo’ essere risolto facendo alcune assunzioni sulle proprietà statistiche delle sorgenti si Dovrebbe essere considerato anche un termine aggiuntivo n per il rumore x=As + n i=1,..,M j=1,..,N

3. Ipotesi per l’ applicabilità di ICA: • N >= M (nel seguito assumiamo M=N senza perdita di generalità) • al massimo una delle sorgenti e’ gaussiana • le sorgenti si sono statisticamente indipendenti • la matrice A ha rango massimo • per ora n=0 , ma il modello puo’ essere esteso anche se di più difficile risoluzione Ambiguità del metodo: Il metodo fornisce una misura dell’ indipendenza delle componenti ma non da informazioni sull’Energia (varianza) e sull’ ordine in cui si ottengono, ovvero la matrice A puo’ essere scritta (dopo la convergenza) come: A=PD P=permutazione D=matrice diagonale Soluzione: Whitening or Sphering

4. PREPROCESSAMENTO DEI DATI: Se i dati si non hanno media nulla allora si sottrae il valor medio (ad xi) Whitening o sphering: Serve ad ottenere dei nuovi dati x con varianza unitaria x=Vx dove E{ xxT }=1 (=I) Se E{ x xT }=C allora V=C-1/2 infatti E{ xxT }= E{V xxT VT}=C-1/2 C C-1/2=I

5. Illustrazione del metodo: Supponiamo di avere due variabili indipendenti uniformemente distribuite nella regione illustrata, con media nulla e varianza unitaria 1/23 |si|<3 Ad es. P(si)= 0 altrove Applichiamo 2 3 A= 2 1 Le direzioni ci danno informazione sulle colonne di A s2 s1 x2 x1

6. … cerchiamo un metodo piu’ generale Procedimento: Si basa sul teorema del “limite centrale”: La distribuzione della somma di variabili random indipendenti tende ad una distribuzione gaussiana Per stimare una delle componenti indipendenti consideriamo y = wT x = i wi xi 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 W e’ una matrice ortogonale WWT=I

7. 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 wTx Misure di non-gaussianità: KURTOSIS kurt(y)=E{y4} – 3(E{y2})2 è nullo per variabili gaussiane quindi si cerca il max di |Kurt(y)| NEGENTROPY 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 INFORMATION I(y1,…,yM) = i H(yi) – H(y) È nulla per variabili indipendenti e non negativa  va minimizzata

8. W Q  x x Modello di rete neurale: y • y è una stima del vettore s y = W x • Q è una stima della matrice A x = Q y • “Apprende” una matrice W tale che y=Wx sono indipendenti • “Apprende” una matrice Q tale da minimizzare • E{||n’||2}=E{||x-Qy||2} 

9. V Q BT x  x x y W=BTV Con pre-withening:

10. 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) RxxL ’’w2wk= -L’wk wk+1 = wk -  Rxx-1L’wk/ r(w) 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

11. FastICA : • Caso semplice “one-unit” • (una sola unità computazionale  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 • Sceglie un iniziale vettore w random • calcola w+ = E{xg (wTx)} – E{g ’(wTx)} w g derivata di una funzione non quadratica • controlla se w = w+ / ||w+|| • se non converge (w  w+ = 1, hanno la stessa direzione) ritorna al punto 2 • Tale algoritmo one-unit permette di determinare solo 1 componente ma può essere facilmente esteso per la stima di più componenti indipendenti improntando una rete “several-unit” con neuroni di pesi w1,…,wn Converge più rapidamente del metodo ICA; non necessita della stima di funzioni g o di parametri di altri parametri, è gratuito e disponibile sul web.

12. Recente applicazione: (Baccigalupi et al. 2002) 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: • 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: • funzione separabile sj(r,) = fj() sj(r) • B(r-r’,)=B(r) indipendente dalla frequenza • aij= t(’)f (’) d’  x(r)=A s(r) * B(r) + n • n è un rumore bianco, indipendente dal segnale , Gaussiano e stazionario

13. Synchrotron angular power spectra Input output

14. CMB angular power spectra Input output

15. 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

16. 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:N pixel Timmagini M sorgenti X = AS • X matrice TxN riga Xt: singola immagine al tempo t colonna Xn: serie temporali (curve di luce) del pixel n • S matrice MxN righe: immagini delle componenti indipendenti per il singolo pixel n Xn = am Smn • A 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 Smncaratterizza il comportamento spaziale

17. T = 35 e N=100x100 pixel dopo whitening 3510 componenti indip. Immagine originaria 1° e 2° autovettori: Raggi cosmici

18. 5° autovettore: Sorgente puntiforme

19. 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