Identificazione di sorgenti gamma tramite ICA
Download
Report
Transcript Identificazione di sorgenti gamma tramite ICA
Formalismo ed applicabilità
del metodo ICA
(Independent Component Analysis)
Francesca Marcucci
Università di Perugia e INFN
Udine 31 gennaio 2003
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 nulla s1 s2 s3 ....... sM ma siano osservabili solo
N combinazioni lineari delle variabili si
xj = aij
· si
i=1,..,M
j=1,..,N
x=As
A={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
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
Soluzione: Whitening or Sphering
P=permutazione D=matrice diagonale
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{ 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:
Supponiamo di avere due variabili indipendenti uniformemente distribuite nella
regione illustrata, con media nulla e varianza unitaria
1/23
|si|<3
s2
Ad es. P(si)=
0 altrove
s1
Applichiamo
A=
2
3
2
1
Le direzioni
x2
ci danno
informazione sulle colonne di A
x1
… 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
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
Modello di rete neurale:
W
Q
x
x
y
y è una stima del vettore s
y=Wx
Q è una stima della matrice A
x=Qy
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}
W=BTV
Con pre-withening:
V
Q
BT
x
x
x
y
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 ’’w = r(w) Rxx
2
L ’’w wk= -L ’wk
2
wk+1 = wk - Rxx-1 L ’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
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
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 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.
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:
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
CMB
angular
power
spectra
Input
output
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:
(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
X matrice TxN
T immagini
M sorgenti
X = AS
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
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 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