UNIVERSITÀ DI TRIESTE

Dipartimento di FisicaCorso di Laurea in Fisica

Tesi di Laurea

Simulazione numerica delproblema di Thomson

generalizzato con applicazioni allereti neurali

Relatoreprof. Marco Budinich

LaureandaSofia Benedetta Rosati

Marzo 2016

Anno accademico 2014/2015

Indice

1 Teoria 31.1 Reti neurali . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Problema di Thomson . . . . . . . . . . . . . . . . . . . . . . . . . 71.3 Dal problema fisico al Neural Relax . . . . . . . . . . . . . . . . . . 91.4 Approfondimenti . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2 Ottimizzazione programma 152.1 Tempo di esecuzione . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.1.1 Costante di apprendimento . . . . . . . . . . . . . . . . . . 152.1.2 Criterio di stop . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.2 Altre osservazioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.2.1 Inizializzazione casuale . . . . . . . . . . . . . . . . . . . . . 222.2.2 Simulazione cariche . . . . . . . . . . . . . . . . . . . . . . 23

3 Analisi dei risultati 253.1 Indipendenza dei risultati . . . . . . . . . . . . . . . . . . . . . . . 25

3.1.1 Binarietà . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.1.2 Distribuzione uniforme . . . . . . . . . . . . . . . . . . . . . 273.1.3 Indipendenza statistica . . . . . . . . . . . . . . . . . . . . 34

3.2 W e trasformazione lineare . . . . . . . . . . . . . . . . . . . . . . 353.3 Conclusioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

A Listato del codice 41

Bibliografia 47

ii

Introduzione

Questa tesi inizia il suo studio dove è finito quello di Elisa Benedetti nel suo articoloNeural Relax [1], che rappresenta, inoltre, insieme alla sua tesi, anche la primaessenziale fonte di informazioni teoriche. Infatti, considerato l’algoritmo elaboratonel suo studio, il Neural Relax per l’appunto, se n’è qui riproposta una provacomputazionale, con lo scopo di ottimizzare la simulazione numerica in prospettivadi applicazioni, rendendola più rapida, e di approfondire l’analisi dei dati ottenuti.

La tesi si articolata in tre capitoli. Nel primo sono introdotte le reti neurali e ilconcetto di apprendimento, il problema di Thomson nella formulazione classica egeneralizzata e la sua utilità per l’algoritmo del Neural Relax, e come quest’ultimosia strutturato, infine vengono accennati alcuni risultati e problematiche teoricheda risolvere che riguardano l’algoritmo stesso. Nel secondo capitolo viene spiegatoil metodo di ottimizzazione della simulazione numerica che prevede sostanzialmen-te un’apprendimento accelerato per la rete, e altre simulazioni effettuate per megliostudiare il problema che hanno avuto, fra gli altri scopi, quello di rispondere alladomanda: ci sono minimi relativi nella funzione di cui si discende il gradiente nel-l’apprendimento? Nel terzo capitolo vengono presentati nel dettaglio gli studi fattiper l’analisi dei dati ottenuti dall’algoritmo, che avevano l’obiettivo di verificareche questi avessero le caratteristiche aspettate dalla teoria.

1

2

Capitolo 1

Teoria

1.1 Reti neuraliIl Perceptron è sostanzialmente il modello matematico di un neurone (McCulloche Pitts, 1943 [2]) con un algoritmo di apprendimento per la classificazione binaria.

Il neurone implementa ovvero una funzione che mappa il proprio input ~x, unvettore reale ∈ Rn, in un output y, a valori in prima istanza binari

y = Θ(n∑i=0

wixi + b) (1.1)

dove ~w è un vettore reale di dimensione n come ~x, e b è la soglia, un numero realeche non dipende da nessun valore di input. La funzione Θ(x) è detta funzione ditrasferimento, e in questo caso è la funzione di Heaviside, che vale 0 se x < 0 e 1per x ≥ 0; altrimenti, più in generale, può essere una qualsiasi funzione che limitil’output fra [0,1] o [-1,1] e che accetti qualunque valore reale in entrata, come latangente iperbolica, che verrà usata successivamente.

Questa modello si può intuitivamente comprendere se si pensa che per ogniinput presentato al neurone, esso ne pesa ogni elemento (attraverso gli elementi delvettore ~w) e deciderà se darà output positivo se in totale questa somma supereràla soglia b (figura 1.1).

Per semplificare la notazione si aggiungerà al vettore ~x un input fittizio chevale sempre 1 (includendolo nelle n dimensioni) e la soglia sarà invece compresanel vettore ~w. Dunque, non esplicitando il prodotto scalare fra vettore dei pesi einput, si può riscrivere l’eq 1.1 come

y = Θ(~w · ~x).

3

1 – Teoria

∑y

x1

w1x2w2

x3w3

xn

wn

1

b

Figura 1.1: Rappresentazione schematica di un neurone.

Sfruttando la rappresentazione vettoriale, ~w è il vettore ortogonale al piano cheripartisce lo spazio degli input fra quelli che danno output positivo e quelli negati-vo, di equazione ~w · ~x = 0 (figura 1.2).

x1

x2

Figura 1.2: Esempio di soluzione nel caso semplice di un input in due dimensioni: laretta rossa rappresenta il piano definito dal peso ~w che divide lo spazio degli inputfra quelli a output positivo e quelli negativo. I cerchi rappresentano un output diun segno, i quadrati di quello opposto.

Dato un insieme di esempi di input ~xν , ν = 1, . . .m, per permettere ad unneurone di imparare una funzione classificatrice, esso viene sottoposto ad una regoladi apprendimento, applicata dopo aver presentato al neurone uno degli esempi, così

4

1.1 – Reti neurali

definita:w′i = wi + ε∆wi (1.2)

dove ogni valore wi viene modificato di una quantità ∆wi al passo successivo;vettorialmente si ha ~w′ = ~w + ε∆~w.1

Conoscendo l’output desiderato ξν per ogni ν-esimo esempio, si può sceglierecome passo di apprendimento

∆wi = (1− ξνyν)ξν , (1.3)

che si annulla per yν = ξν ; in termini vettoriali questa modifica tende ad inclinare ilpiano separante ~w ·~x = 0 in modo tale da includere od escludere l’esempio ν-esimonon appreso (cioè dove vale yν /= ξν).

Il Perceptron è un sistema potente ma non può matematicamente risolvere ogniproblema (definito come il trovare quel peso ~w∗ tale per cui yν = ξν ∀ν), si veda ilcaso semplice dello XOR rappresentato nel caso di input a due dimensioni di figura1.3.

x1

x2

?

Figura 1.3: Esempio di problema di tipo XOR: come si può capire subito non èpossibile trovare un piano che possa dividere gli input in modo tale da averne diun solo tipo da ogni lato.

1Un apprendimento si dice Hebbiano (dallo studio di Hebb, 1949 [4]) se soddisfa la regolamatematica (di ispirazione biologica) secondo la quale il collegamento fra due neuroni si rafforza(o indebolisce) se questi hanno stati simili allo stesso tempo: matematicamente, definiti dueneuroni i e j con stati yi e yj , la modificazione ∆w del peso che li collega w, deve essere del tipo

∆w = εyiyj .

5

1 – Teoria

Rilassando la pretesa di ottenere alla fine dell’apprendimento la giusta rispostaper ogni esempio, si può scegliere di passare ad una funzione di trasferimentocontinua e derivabile (come la tangente iperbolica) e definire la funzione deglierrori

E(~w) = 12

m∑ν=1

(ξν − yν)2 = 12

m∑ν=1

(ξν − f(~xν · ~w))2. (1.4)

Equivalentemente alla definizione di apprendimento 1.3, questa funzione si an-nulla quando tutti gli esempi sono appresi, e se ne può quindi sfruttare la derivataper arrivare al suo minimo

wi = wi + ε∆wi = wi = wi − ε∂E

∂wi

da cui

∂E

∂wi= −

m∑ν=1

(ξν − f(~xν · ~w)) f ′((~xν · ~w)xiν . (1.5)

Il pro di usare questo metodo è che, indipendentemente dalla complessità delsistema e della funzione da apprendere, si arriva ad un minimo che corrispondead un parziale apprendimento da parte della rete degli esempi. Il contro è chesoffre della cosiddetta curse of dimensionality: maggiore è il numero di pesi, moltomaggiori sono i minimi relativi in cui è possibile cadere senza possibilità di uscita.

Si può costruire una rete di neuroni (una cosidetta rete Feed Forward, dove l’in-put di ogni neurone di un certo strato, proviene solo dall’output di strati precedential proprio e mai da quelli successivi), defininendo in modo analogo alla equazione1.4 una funzione che dipende dai pesi di ogni strato (del tipo E( ~W1, . . . , ~Wn) ) ericavare equivalentemente una regola di aggiornamento dei pesi che ne discendela derivata (come in 1.5), ma l’apprendimento si rivelerà sempre più difficile alcrescere della rete perché sempre più grande sarà il numero di minimi relativi.

Però, più in generale, non conoscendo le risposte giuste ξν ma solo una fun-zione che, opportunamente minimizzata, può portarvici, è un metodo standardquello di utilizzare la derivata di questa funzione per portarsi nel minimo (“steepe-st descent”). Sarà proprio questo il metodo usato nell’algoritmo del Neural Relax(sezione 1.3) dove la funzione da minimizzare è il potenziale elettrostatico: grazieall’armonicità di questa funzione, sarà possibile ovviare al problema dei minimirelativi della derivata della funzione da discendere (il potenziale non ha minimirelativi in un insieme aperto e convesso), ma solo in parte, come verrà precisatoalla fine della sezione 1.4.

6

1.2 – Problema di Thomson

1.2 Problema di Thomson

L’obiettivo del problema di Thomson [9] è quello di determinare il minimo del-l’energia potenziale elettrostatica in una configurazione di m cariche uguali vin-colate sulla superficie di una sfera unitaria, sottoposte alla forza di repulsionecoulombiana.

L’interazione elettrostatica presente fra ogni coppia di cariche uguali (Qi =Qj = Q) è data da

Ui,j(~yi, ~yj) = keQ2

|~yi − ~yj |

dove ke è la costante di Coulomb, e |~yi − ~yj | è la distanza fra due cariche i, j, lacui posizione sulla superficie della sfera è definita dai vettori ~y ∈ R3.

Il potenziale complessivo del sistema è

U(~y1, . . . , ~ym) =∑i<j

keq2

|~yi − ~yj |

Semplificheremo d’ora in poi imponendo ke = 1 e che Qi = 1 per tutte le cariche

U(~y1, . . . , ~ym) =∑i<j

1|~yi − ~yj |

. (1.6)

Le soluzioni analitiche di questo problema sono conosciute solo per pochi valoridi m (Schwartz, 2010 [5]), si guardi la figura 1.4 per un esempio.

Figura 1.4: Configurazioni di equilibrio per alcune cariche su una sfera.

7

1 – Teoria

U(~y1, . . . , ~ym) è una funzione armonica 2 (Axler, Bourdon e Ramey, 2001 [7])per cui, per il principio del minimo, non ha minimi in un insieme convesso e apertocome una sfera: il valore di minimo in questo insieme viene assunto dalla funzionesolo sulla sua frontiera. 3

Questo significa che, conseguentemente, la configurazione di riposo delle caricheè con quest’ultime sulla superficie del contenitore: anche guardando un conduttorereale, le cariche possono massimizzare la propria distanza reciproca solo standosulla superficie.

Si generalizza ora il problema di Thomson ad uno spazio a h ≥ 3 dimensioni,e dalla sfera si passa ad un conduttore cubico, di lato due, centrato nell’origine,definito in h dimensioni come

Hh ≡ ~y ∈ Rh : |yk| < 1 k = 1,2, . . . , h (1.7)

con la distanza Euclidea estesa a

|~z − ~t| =

√√√√ h∑k=1

(zk − tk)2 con ~z, ~t ∈ Rh.

In questo nuovo insieme, sempre aperto e convesso, deduciamo che la soluzionecontinui ad essere con le cariche sulla superficie, e che inoltre, se il numero dellecariche è minore o uguale al numero di vertici, allora queste stanno nei vertici.

Per m → ∞ la densità spaziale delle cariche deve essere costante, e necessa-riamente deve essere costante anche la densità di carica superficiale. Perciò, datol’ipercubo Hh, ogni iperpiano passante per l’origine di Rh e che non passa pernessun vertice (per evitare complicazioni) taglia Hh in due parti che contengonolo stesso numero di vertici.4 Con la densità di carica costante, ne consegue che inognuno dei semispazi ci deve essere sempre metà della carica totale.

Siccome questo fatto è valido per ogni iperpiano che attraversi l’orgine di Rh, èvero anche per gli h iperpiani yk = 0, ovvero per tutti e h gli iperpiani normali ad

2Una funzione f : A→ R, dove A è un sottospazio aperto di Rn, è detta armonica, se soddisfal’equazione di Laplace ∇2f = 0.

3Se f ∈ C2(Ω) è armonica, con Ω ⊂ Rn insieme aperto non vuoto con Ω compatto, allora valeche

minΩ

f = min∂Ω

f.

4Infatti ogni vertice appartiene ad uno dei semispazi, mentre quello a coordinate opposteappartiene all’altro. La stessa argomentazione, con le dovute differenze, varrebbe anche avendoun’ipersfera.

8

1.3 – Dal problema fisico al Neural Relax

ognuna delle h coordinate prese singolarmente. Per cui, ricordando che le coordi-nate della cariche saranno binarie, si avranno m

2 cariche con la k-esima coordinatanegativa e altrettante l’avranno positiva.

Quindi, chiamando q(yi) il rapporto fra il numero di cariche con la coordinatai-esima positiva e il totale delle cariche, se m→∞, si avrà che

q(yi) = 12 ∀ i. (1.8)

Definendo l’estensione di q(yi) a t coordinate con q(y1, . . . , yt), per induzione siprova qui che l’equazione 1.8 implica che

q(~y) =h∏i=1

q(yi) con q(yi) = 12 ∀ i. (1.9)

Consideriamo l’iperpiano yi = 0 che seca la distribuzione in due parti ugualie il piano ortogonale a questo yj = 0. Nei quattro sottospazi delimitati dai dueiperpiani, per i ragionamenti precedenti, deve esserci un quarto della carica totalem4 . Per ogni scelta dei due piani definiti da yi e yj si deve ottenere q(yi, yj) = 1

4 ,perciò q(yi, yj) = q(yi)q(yj).

Per continuare il ragionamento per induzione, supponiamo che lo stesso valgacon p iperpiani ortogonali diversi, e cioè che q(y1, y2, . . . , yp) =

∏pj=1 q(yj) = 1

2p .Sfruttando la struttura diHh e (il fatto che le cariche hanno coordinate binarie),

se si aggiunge una (p + 1)-esima coordinata e l’iperpiano di equazione yp+1 = 0,quest’ultimo taglierà tutte le distribuzioni di cariche precedenti in due metà, equindi q(y1, y2, . . . , yp, yp+1) =

∏p+1j=1 q(yj) = 1

2p+1 , completando la dimostrazione.5

1.3 Dal problema fisico al Neural RelaxOra abbiamo tutti gli strumenti per introdurre l’algoritmo del Neural Relax, lacui idea centrale è quella di sfruttare il sistema fisico del problema di Thomsongeneralizzato e una “buona” funzione come quella del potenziale elettrostatico, perottenere un algoritmo di apprendimento per uno strato di Perceptron continui. Ilfine ultimo di questo metodo è quello di avere un processo di pretrattamento di

5C’è però da precisare che solo per m = 2h si può continuare la prova per induzione fino alpasso k = h, dando q(~y) = 2−h per ogni ~y per completare la fattorizzazione della distribuzioneq(~y). Se m < 2h, si può provare solo che tutti i momenti di ordine k di q(~y) sono zero fino ak = log2 m.

9

1 – Teoria

immagini (e più in generale di dati binari), adatto a renderne più efficace l’elabo-razione da parte di memorie associative.

La rete è strutturata con uno strato degli input di dimensione n, e uno stratodi h neuroni per l’output, come illustrato in figura 1.5.

Figura 1.5: Schema della rete

Con un set di m esempi binari, ~xν ∈ ±1n, ν = 1, . . . ,m, l’output della rete,~yν , rappresenta altrettanti vettori che appartengono allo spazio Rh; in particolare,quest’ultimi vengono interpretati come posizioni di cariche, che rimarranno vinco-late all’interno dell’ipercubo h-dimensionale Hh, definito in 1.7. Nella simulazionenumerica gli input sono rappresentati da immagini binarie (figura 3.1). L’outputdi ogni neurone è calcolato con

yi = tanh(~wi · ~x).6 (1.10)L’energia di ogni carica ν-esima, in virtù della repulsione coulombiana, è

Uν =∑

µ=1, µ/=ν

1|~yi − ~yj |

per cui l’energia totale sarà

U(~y1, . . . , ~ym) =m∑ν=1

Uν =∑i<j

1|~yi − ~yj |

.

6Comprendendo in un’unica equazione tutti gli output ottenuti da ogni neurone si userà

~y = tanh(WT ~x)

dove W è una matrice h× n e ~y ∈ Rh.

10

1.3 – Dal problema fisico al Neural Relax

Questa è l’energia che l’algoritmo del Neural Relax minimizza, modificando glielementi della matrice dei pesi W con la discesa del potenziale

w′ij = wij − ε∆wij = wij − ε∂U(~y1, . . . , ~ym)

∂wij

con ε costante piccola; si può quindi calcolare l’aggiornamento

∆wij = − ∂U

∂wij= − ∂

∂wij

(∑µ<ν

1|~yµ − ~yν |

)=

=∑µ<ν

(~yνi − ~yµi)|~yµ − ~yν |2

[xµi(1− y2

µi)− xνj(1− y2νi)].

(1.11)

Lasciando evolvere il sistema grazie alla repulsione coulombiana in Rh, si mi-nimizza l’energia fino a che non si raggiunge l’equilibrio. Nell’ipotesi che m ≤ 2h(cioè avendo meno cariche che vertici), le cariche occuperanno i vertici di Hh. Per-ciò, raggiunto il minimo, si otterranno coordinate binarie in output senza ulterioririchieste e vincoli per l’algoritmo.

Si può notare infatti che in questa regola di apprendimento compaiono i termi-ni del tipo (1− y2

µi): questi tendono a fermare l’apprendimento quando |yµi| ' 1,ovvero alla fine dell’algoritmo le coordinate tendono ad essere binarie.

Siccome la regola di apprendimento sfrutta una funzione armonica, il problemadei minimi relativi (che è tipico dell’usare la discesa del gradiente), è in parte risolto,appunto perché la funzione del potenziale non ne ha in un insieme convesso e chiuso.È risolto però in parte perché le coordinate delle cariche non sono totalmentelibere, visto che è possibile muoverle solo in funzione dei pesi, ma ci si aspetta checomunque questa tecnica possa rendere la funzione da discendere decisamente piùsmooth che con altri metodi (fine della sezione 1.4).

Sapendo che raggiunto il minimo del potenziale la distribuzione delle cariche èuniforme sulla superficie dell’ipercubo (per m → ∞), nella sezione precedente siè dimostrato che il rapporto di cariche rispetto al totale che si troveranno ad unacerta data posizione ~y, sarà fattorizzato rispetto alle coordinate prese singolarmente(eq. 1.9).

Si capisce che quelle che in quel momento si chiamavano solo coordinate, adessosi interpretano anche come output dei vari neuroni. Per cui, quello che abbiamovisto valere per la distribuzione delle cariche, vale anche per la distribuzione diprobabilità associata ai possibili output dei neuroni: la probabilità che si abbia uncerto output finale formato dall’insieme di tutti i neuroni è fattorizzata rispettoalla probabilità associata ad ogni neurone preso singolarmente. E si può utilizzareanche la stessa definizione usata nella sezione precedente, per riferirsi adesso aglioutput piuttosto che alle coordinate.

11

1 – Teoria

Avere una distribuzione fattorizzata implica che le variabili ad essa associatesiano indipendenti: i dati ottenuti da questa rete sono statisticamente indipendenti,per cui

E[yi1yi2 . . . yir ] = 0 ∀ i1 /= i2 /= · · · /= ir, ∀ 1 ≥ r ≥ h. (1.12)

Inoltre da questo deriva che l’informazione associata a questi dati sia massi-mizzata (si guardi la sezione seguente per una dimostrazione che riguarda proprioquesto).

Riassumendo, questo algoritmo è teoricamente capace di

1. dare output praticamente binario, |yµi| ' 1

2. dare una distribuzione di probabilità dell’output fattorizzata (equazione 1.9)

3. dare dati statisticamente indipendenti

4. massimizzare l’entropia dell’output

5. dare una riduzione dimensionale dell’input.

1.4 ApprofondimentiPer completezza, si riporterà qui un risultato teorico, una dimostrazione per sommicapi che tratta la capacità di massimizzare l’informazione dell’output, e un limitedell’algoritmo dovuto al tipo di funzione di trasfermento utilizzata.

Si considera il problema di trovare il minimo dell’energia potenziale U , nel casosemplice di un solo neurone, per cui l’equazione 1.10 diventa

y = tanh(wx+ w0),

con input (per non rendere banale la trattazione del problema) e output continui,con distribuzione di probabilità p(x) per x, e di carica ρ(y) per y.

L’energia del sistema è

U =∫∫

ρ(y)ρ(y′)|y − y′|

dydy′;

chiamando φ(y) ≡∫ ρ(y′)|y−y′|dy

′, il potenziale totale del punto y diventa

U =∫∫

ρ(y)φ(y)|y − y′|

dy ≡∫q(y)dy,

dove q(y) è la densità lineare dell’energia, che è per definizione positiva e propor-zionale al quadrato del campo elettrico (se appropriamente normalizzata q(y) dà

12

1.4 – Approfondimenti

la distribuzione della probabilità di y). Il problema, data x e p(x), è quello dideterminare i parametri w, w0, che minimizzano U .

Utilizzando la regola di trasformazione di una distribuzione di una variabile,p(x) a quella di un’altra variabile rispetto alla quale è legata da un relazione, comey = f(x), in questo caso si ha

q(y) = p(x)|∂f(x)∂x |

.

Sapendo che la soluzione del problema fisico equivalente della distribuzione dicarica che minimizza l’energia potenziale in un filo finito unidimensionale, sarebbedi avere una q(y) costante, ne consegue che necessariamente |∂f(x)

∂x | ∝ p(x) e chela funzione y = f(x) deve soddisfare la seguente relazione di proporzionalità

f(x) ∝∫p(x)dx. (1.13)

Si sa che questo porta alla soluzione per un neurone singolo che massimizzal’entropia (Atick (1992) [3]). Se modificando w e w0 si riesce a mantenere veral’equazione 1.13, allora il sistema minimizza l’energia e massimizza l’informazione.

Nel nostro caso si può riscrivere

tanh′(wx+ w0)|w| ∝ p(x),

sapendo che tanh′(x) > 0 ∀ x.Più in generale se non è possibile soddisfare sempre la 1.13, si può cercare per

quali w si minimizza U .Per cui, partendo da

∂U

∂w= ∂

∂w

∫q(y)dw =

∫∂

∂w

p(x)|∂f(x)∂x |

dx

si può dimostrare che −∂U∂w

= EU [ 1w− 2yw]

− ∂U

∂w0= EU [−2y]

dove EU [x] rappresenta il valore medio dell’argomento sulla distribuzione dell’ener-gia potenziale. Confrontandole con le relazioni di apprendimento che implementanoil celebre algoritmo dell’Analisi delle Componenti Indipendenti, (Bell e Sejnowski,1995 [6]), con le dovute differenze, si può vedere che sono equivalenti, per cui sicapisce che l’algoritmo del Neural Relax è strettamente legato alla ricerca dellecomponenti indipendenti e tende, come nel ICA, a massimizzare l’informazione.

13

1 – Teoria

******

Abbiamo visto che il potenziale U(~y1, . . . , ~ym) è una funzione dellemh coordina-te della cariche e che è dimostrabilmente armonica, ma nel nostro caso, attraversol’equazione 1.10, possiamo cambiare le coordinate delle cariche solo attraverso i(n+ 1)h pesi wij . È semplice verificare che U(wij) non è più armonica:

∂U

∂wij= ∂U

∂yi

∂yi∂wij

,

∂2U

∂w2ij

= ∂2U

∂y2i

(∂yi∂wij

)2+ ∂U

∂yi

∂2yi∂w2

ij

,

e in generale, ∇2U(wij) =∑i,j

∂2U∂w2

ij/= 0. Questo significa che la restrizione

imposta sulla posizione delle cariche ~yν , per il fatto di essere definite con la tangenteiperbolica come ~yν = tanh(WT~xν) (che inoltre rafforza il vincolo di ~yν ∈ Hh) rendela funzione dell’energia non più armonica nelle coordinate “libere” wij .

Questo implica che non possiamo provare formalmente che la funzione U(wij)sia senza minimi locali e che la discesa del gradiente porti sempre alla soluzione cheabbiamo descritto, questo perché possiamo muovere le cariche ~yν solo attraversola variazione dei pesi wij , quindi non liberamente.

Si potrebbe argomentare che è ragionevole aspettarsi che le caratteristiche dellasoluzione non cambino sostanzialmente visto che le cariche sono molto lontane unadall’altra in Hh, e specialmente se m 2h. Ma la forza di una prova formale èpersa. Questo argomento sicuramente merita ulteriori approfondimenti.

14

Capitolo 2

Ottimizzazione programma

Come spiegato nell’introduzione, uno degli obiettivi primari di questa tesi era quellodi ottimizzare la simulazione numerica del Neural Relax [1], perciò si è consideratonecessario per la confrontabilità dei risultati partire dai medesimi dati iniziali edallo stesso algoritmo implementato da Elisa Benedetti.

In tabella 2.1 sono riportate le dimensioni delle principali matrici dell’algoritmo.

Tabella 2.1

matrici costanti~W (pesi) h× n h (neuroni) 64~y (output) h×m m (esempi) 7 †~x (input) n×m n (entrate) 1089+1 ‡

2.1 Tempo di esecuzione2.1.1 Costante di apprendimentoUna volta scritto il programma in Fortran 90, ciò che è stato possibile subito osser-vare nella sua esecuzione, è l’andamento discendente del potenziale U . Nonostantequesto sia fisicamente realistico, è anche poco efficiente: esso tende velocemente

†Nella capitolo 3 verranno usati 20 input.‡Il più uno rappresenta la soglia.

15

2 – Ottimizzazione programma

ad avere una variazione pressoché nulla, rendendo poco effettivo anche un numerogrande (>104) di cicli di apprendimento (figure 2.1 e 2.2).

Figura 2.1: Discesa dell’energia potenziale U durante l’algoritmo in 105 cicli.

Figura 2.2: Variazione del potenziale rispetto al ciclo precedente: superati i primi5000 cicli la variazione non supera 10−6 (le tacche e i gradini indicano che si èraggiunto il limite macchina nella precisione del calcolo).

16

2.1 – Tempo di esecuzione

Questo fatto è cruciale nella misura in cui il tempo di esecuzione del codice èproporzionale al numero di iterazioni: con 106 passi, e circa un decimo di secondoper eseguire un solo ciclo, significa avere un tempo totale di attesa di quasi 40 ore(figura 2.3).

Figura 2.3: Potenziale finale in funzione del logaritmo del numero di cicli.

Cercando quindi un metodo per rendere più rapido l’algoritmo si è rivolta l’at-tenzione verso la matrice dei pesi W , perché essa rappresenta il fulcro della reteneurale -in quanto determina l’attivazione dei neuroni in corrispondenza di ognidato input e conseguentemente decide la discesa o meno del potenziale-, e se n’èstudiato l’andamento e la direzione.

Affinché il potenziale sia minimizzato le m cariche-immigini ridotte rappresen-tate dall’output a h dimensioni, ~yν , devono non solo essere dirette verso i verticigiusti (ovvero dove viene massimizzata la distanza relativa fra le cariche) ma de-vono anche avere ogni componente il più possibile vicina ad uno, in altre parole, lecariche devono stare il più vicino possibile ai vertici (si guardi la sezione 1.2).

Più il prodotto ~wi · ~xν cresce, più rallenta la crescita dell’output verso |1|, pervia della tangente iperbolica che vale uno solo per l’argomento che tende ad infinito(eq. 1.10).

Allo stesso tempo anche l’aggiornamento |∆wij | decresce per via del termine(1 − y2

ν) nella regola dell’apprendimento (eq. 1.11), per cui si ha un feedback po-sitivo nel rallentamento nella discesa del potenziale.

Dunque si è deciso di monitorare la direzione diW , seguendo l’intuizione secon-do la quale essa si assesti all’inizio dell’esecuzione del programma e che la discesadel potenziale sia poi una conseguenza dell’allungamento della matrice dei pesi (equindi un semplice avvicinamento dell’output ai vertici).

Una prima conferma di questo fatto è che, una volta conclusa l’esecuzione delprogramma con un certo U finale, moltiplicando la matrice W finale per un fattore

17

2 – Ottimizzazione programma

qualsiasi maggiore di uno, se si ricalcola l’output, e quindi il potenziale, si ricavaun U ′ che è sempre minore o uguale a U ; maggiore è il fattore moltiplicativo, mag-giore è la diminuizione del potenziale - almeno fino ad un limite oltre il quale nonsi hanno miglioramenti.

Per meglio comprendere perché è importante studiare come si modifica la matri-ce dei pesi, si consideri un vettore ~wi, ovvero il vettore dei pesi associato all’i-esimoneurone, e si scelga una carica, per esempio la ν-esima: il prodotto scalare ~wi · ~xνdetermina l’output dell’i-esimo neurone per quanto riguarda il ν-esimo esempio, eallo stesso tempo questo rappresenta il valore dell’i-esima coordinata della ν-esimacarica, secondo l’equazione

yiν = tanh(~wi · ~xν).

Separiamo ora le direzioni dai moduli, distinguendo i versori (indicati con ilsimbolo ˆ ) dei vettori

yiν = tanh(|~wi||~xν | wi · xν),

siccome vale che sgn(tanh(x)) = sgn(x) allora

sgn(yiν) = sgn(wi · xν)

perché ovviamente |~wi||~xν | > 0.Perciò un cambio di direzione del vettore ~wi tende a cambiare il segno dell’i-

esima coordinata, mentre un cambio di modulo varia solo il valore in modulo dellacoordinata.

Questo significa che seguire la direzione degli h vettori ~wi è molto più signifi-cativo che seguire quella degli m vettori ~yν che indicano la posizione delle cariche.Ammettiamo infatti che si sia raggiunto un passo dell’apprendimento in cui tuttii wi non cambiano più, mentre i moduli |~wi| continuano a variare: i segni dellecoordinate di una carica sono decisi (e quindi è deciso il vertice in cui andrà allafine dell’apprendimento)1, ma il vettore che indica la posizione della carica stessarisulterà di star ruotando e cambiando direzione, non permettendo di capire sel’apprendimento sta procedendo in modo sostanziale o meno. La figura 2.4 è unesempio di questa situazione, semplificata ad un caso di cariche in sole due dimen-sioni: w1 e w2 non stanno variando più, per cui il quadrante in cui la carica si

1Non si sta considerando il caso in cui wi ⊥ xν perché con la condizione che l’algoritmo portile cariche nella configurazione di minimo, si è già illustrato come questa sia raggiungibile solo sele cariche si trovano nei vertici, e se wi · xν = 0 significherebbe che la carica si trova a metà diuno spigolo.

18

2.1 – Tempo di esecuzione

sposta è ben determinato, mentre i moduli |~w1| e|~w2| stanno continuando a modi-ficare la direzione del vettore ~yν .

y1

y2

0

(1,1)

12

3

Figura 2.4: Variazione della posizione di una carica in due dimensioni all’internodi un solo quadrante; i numeri indicano passi successivi dell’apprendimento.

Ad ogni modo, per misurare la variazione degli angoli degli h vettori ~wi a ncomponenti, relativi ad ognuno dei neuroni, se n’è trovato il vettore normalizzato,wi (senza poi effettivamente normalizzare ~wi). Al passo successivo del programmasi è fatto lo stesso con i vettori aggiornati ~w′

i → w′

i, e per la regola dei coseni siè ricavato l’angolo θi compreso fra ogni coppia i-esima dei vettori a norma uno(figura 2.5):

|w′

i − wi|2 = |wi|2 + |w′

i|2 − 2|wi| · |w′

i| cos θi = 2(1− cos θi).

~wi

~w′i∆~wi

wi

w′i |w′i − wi|θi

Figura 2.5

19

2 – Ottimizzazione programma

Trovato perciò un modo per monitorare la variazione della direzione e dellalunghezza degli ~wi, se n’è considerato per semplicità uno solo (avendo precedente-mente osservato che la variazione era simile per tutti) e si è graficato il risultatoper un numero alto di clici (figure 2.6a e 2.6b)

Figura 2.6

(a) Tipico grafico della variazione totale del-l’angolo (in radianti) dall’inizio del program-ma, in funzione del ciclo. Dopo le primemigliaia di cicli l’angolo non cambia più.

(b) Particolare dei primi 2000 cicli con la va-riazione dell’angolo (in radianti) ad ogni ciclorispetto al precedente.

Come si può vedere si ha avuto una conferma dell’intuizione e si è dunquecercato un metodo per consentire ai vettori dei pesi di orientarsi inizialmente con lamaggiore libertà possibile, e, in un secondo momento, permettendogli di allungarsivelocemente.

La soluzione si è rivelata essere più semplice del previsto, visto che ciò che deter-mina proprio il passo di apprendimento è la costante ε: rendendola una variabile,crescente, si è ottenuto esattamente questo risultato (figure 2.8a e 2.8b).

Numericamente è risultato sufficiente per l’assestamento del vettore dei pesi unacrescita quadratica dell’ε (figura 2.7)- od anche inizialmente quadratica e, dopo uncerto numero di cicli, cubica.

Con questo metodo è bastato un numero di cicli dell’ordine 102-103 per ottenererisultati sostanzialmente migliori (da 1.82 a 1.76), guadagnando un fattore 103 intermini di tempo risparmiato (ovvero da alcuni giorni a circa un minuto).

2.1.2 Criterio di stopUn altro modo per rendere più veloce il codice è stato quello di ridefinire il criterioper fermare il programma.

Un metodo standard per arrestare gli algoritmi di apprendimento delle retineurali è quello di monitorare il più grande elemento nel vettore dell’aggiornamento

20

2.1 – Tempo di esecuzione

Figura 2.7: Discesa del potenziale in 500 cicli con ε crescente quadraticamente.

Figura 2.8

(a) Variazione totale dell’angolo (in radianti)dall’inizio del programma di un ~wi. Dopo circa400 cicli l’angolo non cambia.

(b) Variazione dell’angolo di un ciclo rispettoal precedente.

dei pesi ∆wij : se questo diventa minore di una certa soglia sufficientemente piccola,si deduce che l’apprendimento non apporta più modificazioni sostanziali.

Bisogna però notare che l’elemento ∆wij indica in generale la modificazioneche viene fatta al vettore ~wi, senza però specificare quanto questo stia ruotando oallungandosi.

Come già detto, in questo caso il fattore fondamentale era l’orientamento di ~wi,piuttosto che il suo modulo: sebbene questo debba essere sufficientemente grandeper avere un output vicino ad uno, si deve prima di tutto sapere in che vertici

21

2 – Ottimizzazione programma

dovranno stare le cariche.Questo suggerisce che si può scegliere come regola di stop una variazione del-

l’angolo nulla: ulteriori iterazioni tendono solo ad allungare gli ~wi. Nel caso già mo-strato in figura 2.7, sarebbero bastati i primi 400 cicli, permettendo di guadagnareun altro 25% sul tempo di esecuzione.

2.2 Altre osservazioni

2.2.1 Inizializzazione casuale

Avendo sostanzialmente ridotto il tempo di esecuzione del programma è stato pos-sibile usare tante inizializzazioni diverse della matrice W . Con valori troppo alti(> 0.05) o troppo bassi (< 10−5) per i vari elementi di W , l’apprendimento erapressoché impossibile, ma rimanendo entro l’intervallo adatto, quello che si è os-servato è che il potenziale ottenibile dipende in certa misura dagli elementi inizialidella matrice dei pesi.

Il programma è stato scritto in Fortran 90 e la funzione per la generazio-ne dei numeri (pseudo)casuali con distribuzione uniforme era quella di libreria,random_number(), deterministica una volta scelto il vettore di numeri interi perinizializzarla, per cui, per cambiare la matrice dei pesi iniziale da un’esecuzioneall’altra, si è dunque cambiato questo vettore.

Guardando il grafico 2.9, che raccoglie alcune migliaia di discese del potenziale,si può vedere che esse sono tutte molto simili. Sembrano anzi quasi la traslazioneverticale di un solo andamento, indicando che la capacità dell’algoritmo NR diminimizzare la funzione del potenziale elettrostatico, sia non sono limitata maanche sempre la stessa, per cui, dato il punto di inizio, si può praticamente giàprevedere il punto finale. Sicuramente questo mette in luce il fatto che ε è grande(almeno dopo un certo numero di cicli), perché l’algoritmo, avvicinandosi ad unminimo, tende a non uscirne. Ad ogni modo nella scrittura del programma si è dataprimaria importanza al rapporto qualità della soluzione/tempi di esecuzione e vistoche il potenziale finale che si otteneva prima era peggiore di più della metà di quelliraccolti nel grafico 2.9 e che inoltre i tempi di esecuzione erano enormemente piùlunghi, si pensa che si possa considerare un male minore una certa limitatezza nellacapacità di minimizzare il potenziale. Saranno di sicuro necessarie più accuratecalibrazioni dei parametri dell’algoritmo per ottenere il miglior risultato nel minortempo, e probabilmente dovrà essere presa in esame la possibilità di aumentare laprecisione macchina e il numero di cifre significative delle variabili.

22

2.2 – Altre osservazioni

Figura 2.9: Discesa del potenziale in 60 cicli con 4000 diverse inizializzazioni di W .

2.2.2 Simulazione carichePer avere un riscontro sui risultati ottenuti con la rete neurale si è pensato auna simulazione esclusivamente fisica: cariche uguali confinate in un ipercubo ah dimensioni (definizione 1.7) sottoposte solo alla forza repulsiva elettromagnetica(lo stesso potenziale di equazione 1.6).

Per quanto riguarda il potenziale si è visto che i valori finali tipici stavanonegli stessi intervalli ottenuti con la rete neurale, e così anche la discesa stessa delpotenziale ha lo stesso tipo di andamento.

Invece, per quanto riguarda la disposizione finale delle cariche, sono stati usatigli stessi test illustrati nella prossima sezione e non si sono riscontrate particolaridifferenze (se non giusto risultati migliori in corrispondenza di potenziali minori),per cui si ritiene di aver avuto la conferma da questa verifica che la simulazionecon le reti neurali implementi un algoritmo “fisicamente” realistico.

23

24

Capitolo 3

Analisi dei risultati

3.1 Indipendenza dei risultatiL’obiettivo dell’algoritmo del Neural Relax è quello di ottenere una trasformazio-ne lineare per dati binari che li renda indipendenti e scorrelati (massimizzandol’informazione) con il vincolo però che questa li lasci binari.

Per la simulazione numerica si è partiti da una matrice, x, che rappresenta unset di m = 20 immagini binarie di 33x33 pixel,1 da cui se ne ottiene un’altra, y, cheinvece ha solo h = 64 dimensioni per ognuna delle venti immagini (figura 3.1); adogni modo in tabella 2.1, sono ricordate tutte le dimensioni delle matrici principali.

I controlli che sono stati effettuati sull’output sono finalizzati a capire se i datiottenuti rispettano tre principali caratteristiche:

1. binarietà

2. distribuzione uniforme delle cariche nello spazio

3. indipendenza statistica dell’output.

3.1.1 BinarietàPer la binarietà basta un semplice istogramma dell’output, in figura 3.2a ce n’è unesempio ottenuto dopo 200 cicli di apprendimento. Come si può vedere anche unnumero così basso (relativamente all’ordine di 106 cicli che erano necessari prima)basta per avere output pressoché binario. In figura 3.2b c’è l’istogramma del valore

1Ogni immagine viene resa una stringa di 332 + 1 = 1089 + 1 bit; il +1 è per la soglia.

25

3 – Analisi dei risultati

Figura 3.1: Immagini usate per l’apprendimento e poi l’analisi statisticadell’output.

assoluto dell’output, dove è chiaro quanto si distanzino i vari bit da 1, il più lontanoè a 0.9976.

Figura 3.2

(a) Istogramma dell’output dopo un ap-prendimento di 200 cicli.

(b) Istogramma dell’output in valoreassoluto.

26

3.1 – Indipendenza dei risultati

3.1.2 Distribuzione uniforme

Per verificare che le cariche fossero in distribuzione uniforme sono stati fatti dueprincipali tipi di analisi diverse.

Matrice Q

Come prima cosa si è guardata la matrice degli overlap Qνµ ≡ yν · yµ (figura3.3a). Interpretandola in senso “fisico” ogni suo elemento rappresenta la distanzareciproca fra ogni coppia di cariche: chiamando k il numero di coordinate in comunefra una certa immagine ~yµ e un’altra ~yν , allora ~yµ · ~yν = 2k − h e la distanza chele separa è d = 2

√h− k. 2

Per figurarsi cosa si dovrebbe trovare nel caso ottimale si possono fare alcuneosservazioni preliminari: se si avesse un numero di cariche uguale al numero divertici (m = 2h) si dovrebbero trovare nell’istogramma tutte le possibili distanzedi un vertice rispetto a tutti gli altri, in una distribuzione simmetrica rispetto azero (in particolare, si avrebbe la distribuzione binomiale di B(64,1/2)); se m = 2,a prescindere dalla dimensione dell’ipercubo, queste dovrebbero stare in verticidiametralmente opposti, dando un istogramma (con i vettori normalizzati) con duecolonne alte due centrate una su −1 e una su 1; se m . 2h, è ragionevole pensareche le cariche non abbiano sufficiente spazio per allontanarsi molto una dall’altra,dando perciò un istogramma con molti valori positivi; via via che si riduce questonumero di cariche allo stesso modo ci si può aspettare che (escluso il picco dei valorisulla diagonale di Q che valgono sempre 1) la media della distribuzione tenda aspostarsi sempre più verso valori negativi, questo perché le cariche sono in mediasempre più distanti una dall’altra.

Quello che comunque è possibile immaginare per la configurazione di ottimoè che la distanza che intercorre fra una carica e la sua prima vicina dovrà esserela stessa di tutte le altre cariche rispetto alla rispettiva prima vicina (e così viaper la seconda vicina, la terza ecc.), affinché la figura geometrica che si forma siail più possibile regolare. Per cui guardando la figura 3.3b (che è semplicementela matrice Q dove sono state ordinate le righe dal valore più basso al più alto)si immaginerebbe di trovare colonne con gli stessi numeri: la prima vicina è allastessa distanza per tutte (e così per la seconda vicina, terza ecc.); allo stesso modosi dovrebbero trovare righe tutte uguali fra loro, mentre come si può vedere infigura 3.4 si hanno vari andamenti.

2Con una sola coordinata non in comune, d = 2, ovvero avranno una distanza che le separaequivalente ad un lato dell’ipercubo; con una differenza di due bit, d = 2

√2: avranno una

diagonale di una faccia come distanza di separazione ecc.

27

3 – Analisi dei risultati

Figura 3.3

(a) Matrice Q degli overlap.(b) Matrice Q con righe ordinate dalvalore più basso al più alto.

(c) Legen-da dei co-lori.

(d) Istogramma di Q.

(e) Valori medi (cerchi) e deviazione standard(quadratini) di ogni colonna di Q dopo esserestata ordinata.

Dal momento che queste distanze non sono uguali questa prova non è moltosignificativa (visto che si perdono le associazioni e le cariche diventano irriconosci-bili), se non nella misura in cui non si notano particolari anomalie. In figura 3.6si ha l’istogramma di alcune delle colonne della matrice Q ordinata (mentre in 3.5si hanno i contributi di ogni colonna separata): si può notare, guardando i numeriche si riferiscono alla colonna, che gli elementi delle colonne si distribuiscono inmodo circa gaussiano, con un picco centrale; nel grafico 3.3e si può trovare invecela devizione standard e la media calcolate sui valori di tutte le colonne: con rife-rimento alla deviazione standard, si può vedere che l’estensione dei valori è circa

28

3.1 – Indipendenza dei risultati

Figura 3.4: Profili delle righe della matrice Q ordinata.

Figura 3.5: Istogramma delle colonne della matrice Q dopo essere stata ordinata,dove sono evidenziati i contributi di ogni colonna (indicate dai numeri). Si notiche il profilo del grafico è lo stesso di figura 3.3d.

29

3 – Analisi dei risultati

Figura 3.6: Istogrammi di alcune colonne della matrice Q.

(a) IV colonna. (b) VII colonna.

(c) X colonna. (d) XIII colonna.

(e) XVI colonna. (f) XIX colonna.

sempre la stessa.

Considerando però la configurazione di ottimo per un numero di cariche m ≤h+ 1, sembra ammissibile pensare che le cariche abbiano tutto lo spazio necessarioper disporsi con distanze reciproche quanto più simili, equilibrando così il piùpossibile le forze di repulsione presenti fra di loro.

Questo ragionamento è motivato dal fatto che in h dimensioni si possono co-struire simplessi regolari (politopi regolari con il minor numero di vertici, che sonola generalizzazione del tetraedro, dove ogni distanza reciproca fra i vertici è la

30

3.1 – Indipendenza dei risultati

stessa) fino a quello h dimensionale, che ha h+ 1 vertici. Queste figure sono peròinscrivibili in ipersfere e non è detto che si possa far combaciare i loro vertici conquelli di un ipercubo, ad ogni modo si ritiene che il lato del simplessom ≤ h dimen-sionale (che ha m+ 1 vertici) dia un suggerimento per la distanza media ottimalefra le cariche. In particolare è utile sapere che per un m-simplesso regolare si hache il prodotto scalare fra qualunque coppia di vettori normalizzati che indicano laposizione dei sui vertici vale sempre − 1

m .Per cui si guardino i grafici di figura 3.7, che provengono da due esecuzioni del

programma diversi, dove è segnato il valore medio della distribuzione (in blu) e− 1

20 = −0.05 (in rosso): nel grafico 3.7a il valore finale del potenziale era maggioredi quello nel caso della figura 3.7b, e si può notare che

1. la distribuzione del secondo è molto più concentrata del primo

2. il valore medio della distribuzione e − 1m sono molto più vicini nel secondo

caso (quasi sovrapposti)

3. in ogni caso il valore medio è più alto di − 1m .

Ciò significa che allo scendere del potenziale si ha che:

1. la distanza reciproca fra le cariche tende ad uniformarsi verso un valore unico

2. la distanza media si avvicina ad un valore ottimale, che si ritiene essereragionevolmente stimato dal lato dell’m-simplesso

3. la distanza reciproca media delle cariche tende ad aumentare, e quindi intermini di prodotto scalare, questi tendono a valori più negativi.

Figura 3.7: Istogrammi di Q in due esecuzioni del programma diverse. La freccia inblu indica il valore medio dell’istogramma (esclusa la colonna in 1), mentre quellain rosso indica − 1

20 .

(a) Uf = 17.14, valore medio= −0.0116 (b) Uf = 16.48, valore medio= −0.0463

31

3 – Analisi dei risultati

Tutto questo può dare un’idea di come si distribuiscono le cariche, ma non sene ha ancora una visione chiara e precisa.

Piani secanti

Una prova che ci pare più significativa è quella di testare la soluzione con una seriedi piani passanti per l’origine di Hh generati casualmente: se le cariche fosseronella configurazione ottima metà dei piani dovrebbero trovare ogni carica solo daun lato.

Per fare questo si sono generati semplicemente dei vettori h-dimensionali, peridentificare i piani, e si è considerato il prodotto scalare fra questi e la posizionedelle cariche per distinguere il lato di appartenenza.

Inizialmente si è guardata la somma totale delle cariche che si sono trovate daun lato dei piani: con 105 piani, e svariati set di piani diversi, si è trovato che ilrapporto è sempre rimasto in 0.498, 0.501.

Per vedere se fosse possibile essere più accurati in queste verifiche si è osservatoquante volte in totale ogni carica presa singolarmente si era trovata da un lato di unpiano piuttosto che dall’altro. Si è osservato (figura 3.8) che, come ci si aspettava,ogni carica si è trovata all’incirca metà delle volte da un lato (si noti che la scalain ordinata va da 0.48 a 0.52).

Figura 3.8: Rapporto di quante volte ogni carica presa singolarmente si è trovatada un lato dei piani, sul totale dei piani.

Ad ogni modo, la significatività di queste verifiche ci pare limitata. Una cheinvece crediamo possa essere la più determinante è quella di accertare, preso ognipiano, quante cariche esso si è trovato da un lato, piuttosto che dall’altro, eistogrammando questo valore.

32

3.1 – Indipendenza dei risultati

Figura 3.9: In ascissa ci sono il numero di cariche che si sono trovate da un la-to dei piani, in ordinata invece c’è il numero di volte che si è trovata una certaconfigurazione.

Qualunque sia il numero di cariche, esistono in totale 2m diverse configurazioni(fatte dividendo le cariche in due gruppi di appartenenza), ma solo nel caso parti-colare di avere m ≤ h tutte queste diverse configurazioni possono essere implemen-tate da iperpiani nello spazio a h dimensioni con il vincolo che questi passino perl’origine.3 Queste 2m configurazioni seguono una distribuzione binomiale B(m, 1

2 ).Perciò, avendo m ≤ h, qualunque sia la distribuzione delle cariche nello spazio,

esistono sempre quei 2m piani (trovarli poi è più difficile) che dividono le cariche intutte le 2m configurazioni possibili, dando un istogramma esattamente binomiale[8, p. 41-44.]. Ma se le cariche sono concentrate in un sottospazio dell’ipercubo,i piani che implementano tutte le configurazioni possibili saranno uno specifico eparticolarissimo sottoset di tutti i piani che passano per l’origine di Hh.

Dunque, per assicurarsi che m ≤ h cariche siano ben distribuite nello spazio,bisogna prendere piani distribuiti uniformemente sull’ipersfera: se con questi sitrova un istogramma delle configurazioni il più possibile simile ad una binomialeallora è ragionevole pensare che le cariche siano in distribuite uniformemente.

3Per m ≥ h certe configurazioni delle 2m totali non sarebbero possibili da ottenere con unpiano, perché potrebbero esistere problemi di tipo XOR (figura 1.3).

33

3 – Analisi dei risultati

Essendo 20 le immagini, si può già approssimare la distribuzione binomiale conuna gaussiana N(m2 ,

√m2 ), e in figura 3.9, ottenuta dall’analisi dell’output di un

apprendimento di 200 passi, si può osservare che la corrispondenza in questo casoè buona.

Si è poi memorizzata una serie di output ottenuti durante l’apprendimento, esi è utilizzato lo stesso set di piani per verificare la distribuzione delle cariche neltempo. Questa prova è stata fatta per un apprendimento di 500 passi, ma si èvisto che dopo circa 200 passi i valori non si modificavano, per cui, per rendereapprezzabile la modificazione, a titolo di esempio si mostrano qui, in figura 3.10,i primi 30 passi di un apprendimento, con intervalli di 5 passi fra uno e l’altro.Come era ragionevole aspettarsi, la distribuzione si stringe e si alza nel tempo,assomigliando sempre più ad una binomiale, per cui si capisce che in corrisponden-za della discesa del potenziale si ha una sempre più uniforme distribuzione dellecariche nello spazio. L’osservazione che non ci siano particolari modificazioni (senon fluttuazioni numeriche intorno a dei valori fissi) dopo un certo numero di cicli,suggerisce che questo potrebbe essere un altro criterio di arresto del programma,sebbene chiedersi se la posizione delle cariche sia rimasta la stessa da un passodell’apprendimento all’altro (come si verifica qua), sia equivalente a chiedersi se ladirezione dei vettori dei pesi non si siano modificate (come nel metodo spiegatonella sezione 2.1.2); sono necessarie ulteriori ricerche per verificare quale metodosia il più performante e affidabile.

3.1.3 Indipendenza statisticaPer l’indipendenza statistica, si è verificato che

E[yi1yi2 . . . yir ] = 0 ∀ i1 /= i2 /= · · · /= ir, ∀ 1 ≥ r ≥ h (3.1)

dove

E[yir ] = 1m

m∑ν=1

yνir

dove si assume che la probabilità associata ad ogni bit sia 1m . Per farlo si sono

calcolati i valori medi di tutte le combinazioni possibili dei prodotti yi1yi2 . . . yirfino ad un massimo di tre elementi, poi solo un campione di 105 per le combinazionidi 4,5 e 6 elementi e si sono istogrammate le distribuzioni dei valori. Si è calcolataanche la covarianza, che dovrebbe essere nulla come tutti gli altri momenti: ciòsignifica che gli output non solo sono indipendenti, ma anche scorrelati. Ci siè fermati all’ordine sei, sebbene fosse necessario per la verifica dell’indipendenzastatistica arrivare fino al grado massimo, ovvero 64, perché si è visto che non c’eranoparticolari cambiamenti dal quarto in poi.

34

3.2 – W e trasformazione lineare

Figura 3.10: I numeri nella legenda indicano il ciclo a cui è stato fatto il grafico. Ipunti sono stati uniti solo per maggiore chiarezza.

Quello che ci interessa osservare dai grafici raccolti in figura 3.11 è che tuttequeste distribuzioni sono centrate sullo zero, come si voleva in equazione 3.1.

Si può considerare ogni bit di ogni singolo vettore ~yν dell’output come unavariabile aleatoria con valore medio uguale a zero (perché la probabilità che valga1 è uguale a quella di ottenere -1). Di conseguenza, fare la somma di alcunedi queste variabili, ricorda il problema del random walk unidimensionale, dove ilnumero di immagini è il numero di passi fatti, e la probabilità associata ad unvalore negativo o positivo è proprio 1

2 . Per cui, la variabile aleatoria associata allasomma di uno stesso bit (o del prodotto di due, tre ecc.) di tutte le immagini, hadistribuzione binomiale. Come già detto con 20 immagini inizia ad essere già buonal’approssimazione della binomiale con una gaussiana, per questo motivo, nei grafici3.11, sono stati sovrapposte agli istogrammi delle gaussiane con µ = 0 e σ = 1√

20 :come si vede, dal terzo compreso in poi (dove si inizia ad avere abbastanza numeriper una buona statistica), gli istogrammi hanno una ottima corrispondenza con legaussiane.

3.2 W e trasformazione lineareAlla fine si riportano qui delle osservazioni interessanti che ancora non è possibilespiegare.

35

3 – Analisi dei risultati

Se si considera la matrice dei pesi W e quella dell’output y e se ne fa il pro-dotto WT y quello che si trova, riarrangiando i vettori in matrici quadrate, sono leimmagini di input, modificate, ma perfettamente riconoscibili (figura 3.12), fattonon aspettato.

Sono poi notevoli gli istogrammi degli elementi di W (3.13a), delle norme deipesi associati ad ogni neurone (3.13b) e quello dei valori che si ottengono se si fa ilprodotto della matrice W con la trasposta dell’input: il primo perché la distribu-zione che si trova ha un profilo molto particolare, poi ci dice che la maggior partedei valori di W sono nulli e che nel complesso sono relativamente “piccoli”; il se-condo perché si trova che tutte le norme dei vettori dei pesi ~wi stanno in un piccolointervallo in una distribuzione molto piccata; e il terzo in particolare è interessanteperché si ottiene una distribuzione pressoché con due picchi a valori opposti (comesi cerca di avere nel caso di y). Quest’ultimo grafico sembra suggerire l’intuizioneche sia in qualche modo possibile prescindere dalla tangente iperbolica della defini-zione dell’output 1.10, perché anche senza si ritrova una trasformazione che lascial’input binario. Senza tangente iperbolica si avrebbe una trasformazione per l’out-put perfettamente lineare che verosimilmente renderebbe la funzione minimizzataarmonica, cancellando così il problema dei minimi relativi.

3.3 Conclusioni

Per concludere si riassume il percorso fatto e i risultati ottenuti.Si sono inizialmente esposti i componenti che costituiscono il Neural Relax: le

reti neurali e il problema di Thomson generalizzato; poi come e cosa si prefiggedi raggiungere: ottenere una trasformazione per dati binari che li renda scorrelati,indipendenti statisticamente e che li lasci binari. Infine si è mostrato che il NeuralRelax massimizza l’informazione e poi come non è possibile escludere teoricamen-te che la funzione da minimizzare nelle variabili “libere” wij sia priva di minimirelativi.

Nella seconda parte si è studiata l’evoluzione delle variabili dell’algoritmo du-rante l’apprendimento e si è proposto un metodo per ridurre i tempi d’esecuzioneche consiste nel far crescere la costante di apprendimento ε; inoltre, nelle simula-zioni, non si è trovata l’evidenza di minimi relativi. Ad ogni modo, è sicuro che leottimizzazioni numeriche hanno ancora molto margine per migliorare l’algoritmo,ma si ritiene che sia già un buon risultato quello di aver ridotto i tempi d’attesa diun fattore 103.

Nella terza parte si è illustrato tutta la serie di ragionamenti che hanno appro-fondito la nostra comprensione dei dati che si ottengono da questo algoritmo. Inprimo luogo infatti, si sono elaborati dei test che permettessero di indagare le pro-prietà che ci interessavano dei dati, e con essi si è cercato di definire la loro giusta

36

3.3 – Conclusioni

interpretazione, e, in secondo luogo, nessuno di essi ha dato risultati in contraddi-zione con le caratteristiche aspettate, per cui non si ha avuto motivo di smentirenessuno dei presupposti teorici.

Ad ogni modo è sicuro che ci siano molti argomenti che meritano maggioriapprofondimenti, e alcuni che sembrano avere grande potenzialità, perciò si ritieneche ci sia ancora molto spazio per nuove analisi e si auspica che queste si possanoconcretizzare in future ricerche.

37

3 – Analisi dei risultati

Figura 3.11

(a) E[yi1 ] (b) E[yi1 yi2 ]

(c) E[yi1 yi2 yi3 ] (d) E[yi1 yi2 yi3 yi4 ]

(e) E[yi1 . . . yi5 ] (f) E[yi1 . . . yi6 ]

(g) cov[yi, yj ]

38

3.3 – Conclusioni

Figura 3.12: Immagini ottenute dal prodotto WT · y.

Figura 3.13

(a) Istogramma degli elementi diW .

(b) Istogramma delle norme dei 64vettori ~wi.

(c) Istogramma dei valori ottenuti dalprodotto W · xT .

39

40

Appendice A

Listato del codice

Di seguito si riporta il codice scritto in Fortran 90 per implementare l’algoritmodel Neural Relax.

1 module m12 implicit none3

4 ! Costanti globali5

6 integer , parameter :: n =10907 integer , parameter :: h=648 integer , parameter :: m=209 integer , parameter :: maxloops =100

10 integer , parameter :: seed_size =1211

12 real , parameter :: halt =0.00000113 real , parameter :: giga =100.14

15 ! Variabili globali16 real :: start1 , stop1 , nowtime17 real :: eps18

19 real , dimension (n,m) :: x20 real , dimension (h,n) :: w21 real , dimension (h,n) :: dw22 real , dimension (n,h) :: dwT23 real , dimension (h,m) :: y24 real , dimension (h) :: dwStop25

26 integer , dimension ( seed_size ) :: vector =(/ 651 , 138 , 982 , &27 554 , 570 , 89, 13, 892 , 758 , 788 , 868 , 951/)28

29

30 contains31

41

A – Listato del codice

32 ! Inizializzazione vettore degli input , dei pesi e degli output33 subroutine Initialize ()34 integer :: i, j35 real , dimension (h,n) :: t36 real :: r37

38 do i=1,m39 read(unit =10 , fmt =*) x(:,i)40 end do41

42 call random_seed (put= vector )43

44 do i=1,h45 do j=1,n46 call random_number (r)47 w(i,j)=(r -0.5)/ giga48 end do49 end do50

51 w(:,n)=0.52 x(n ,:)=1.53 y=0.54 U1 =0.55 eps=epsi56

57 end subroutine Initialize58

59

60 ! Calcolo dell ’ output61 subroutine OutputUpdate ()62 integer :: i,j63

64 y=0.65

66 do i=1,h67 do j=1,m68 y(i,j)= tanh( dot_product (w(i ,:) ,x(:,j)))69 end do70 end do71

72 end subroutine OutputUpdate73

74

75 ! Calcolo del potenziale76 function UCalc () result (U)77 real :: U78 integer :: indmu , indnu79

80 U=0.81

82 do indmu =1,m -183 do indnu = indmu +1,m

42

84 U=U+1./ Dist(indmu , indnu )85 end do86 end do87

88 end function UCalc89

90

91 ! Calcolo della distanza92 function Dist(indmu , indnu ) result (d)93 integer :: indmu , indnu94 integer :: i95 real :: d96

97 d=0.98 d=sqrt( sum( (y(:, indmu )-y(:, indnu ))**2 ) )99

100 end function Dist101

102

103 ! Regola di aggiornamento104 subroutine LR(indmu , indnu )105 integer , intent (in) :: indmu , indnu106 integer :: i,j107

108 do i=1,h109 do j=1,n110 dw(i,j)= dw(i,j)- (y(i, indmu )-y(i, indnu )) * &111 ((1. -y(i, indmu )**2)* x(j, indmu ) - &112 (1. -y(i, indnu )**2)* x(j, indnu ) ) / &113 (Dist(indnu , indmu ))**3114 end do115 end do116

117 end subroutine LR118

119

120 end module m1121

122

123 program NR124 use m1125 implicit none126

127 ! Variabili locali128 integer :: i, j, k129 integer :: indnu , indmu130 integer :: progress =1131 real :: U1 ,U2 , Uvar132 real :: maxvar133 real :: step =0.00001134 real :: epsi =0.0000001135 real :: tot =0

43

A – Listato del codice

136

137 real , dimension (h,n) :: wprec =0.138 real , dimension (h,n) :: wvar =0.139 real , dimension (h) :: wnorm140 real , dimension (h) :: wprecnorm141

142

143 !File di scrittura / lettura144 open(unit =10 , file="images -20. txt", action ="read")145 open(unit =11 , file="wI", action ="read")146 open(unit =12 , file="wF", action =" write ")147 open(unit =13 , file="dw", action =" write ")148 open(unit =14 , file="y", action =" write ")149 open(unit =15 , file="UW", action =" write ")150 open(unit =23 , file=" angoloW ", action =" write ")151

152 call cpu_time ( start1 )153 call Initialize ()154

155 print *, " ********* ******** ********* "156 print *, " ********* NR ********* "157 print *, " avanzamento tempo (s) loop U attuale "158

159 ! Algoritmo di apprendimento160 do i=1, maxloops161

162 eps=epsi+step*i**2163

164 if (mod(i *100./ maxloops ,5.)==0) then165 call cpu_time ( nowtime )166 print *, 5* progress , "% ", nowtime -start1 , i, Ucalc ()167 progress = progress +1168 end if169

170 U1=U2171 wprec =w172

173 call OutputUpdate ()174

175 dw =0.176 dwT =0.177 dwStop =0.178

179 U2= Ucalc ()180 Uvar=U2 -U1181

182 do indmu =1,m -1183 do indnu = indmu +1,m184 call LR(indmu , indnu )185 end do186 end do187

44

188 w=w-eps*dw189

190 do j=1,h191 wnorm (j)= sqrt( dot_product (w(j ,:) ,w(j ,:)))192 wprecnorm (j)= sqrt( dot_product ( wprec (j ,:) , wprec (j ,:)))193 wvar(j ,:)=w(j ,:)/ wnorm (j)- wprec (j ,:)/ wprecnorm (j)194 end do195

196 do j=1,h197 dwStop (j)= sqrt( dot_product (wvar(j ,:) , wvar(j ,:)))198 end do199

200 maxvar = maxval ( dwStop )201 tot=acos (1. -( maxvar **2)/2.)+ tot202

203 call cpu_time ( nowtime )204 write (unit =15 , fmt =*) i, Ucalc (),Uvar , maxval (dw), maxval (w), &205 eps , nowtime - start1206 write (unit =23 , fmt =*) i, 1. -( maxvar **2)/2. , &207 acos (1. -( maxvar **2)/2.) , tot208

209 end do210

211 ! Scrittura sui file212 do j=1,h213 write (unit =12 , fmt =*) w(j ,:)214 write (unit =13 , fmt =*) dw(j ,:)215 write (unit =14 , fmt =*) y(j ,:)216 end do217

218 print *, " "219 print *, " totale cicli ", maxloops220 print *, "U fine programma =", Ucalc ()221

222 w=w *10**8223 call OutputUpdate ()224 print *, "U con W_fin *10^8= ", Ucalc ()225 y=int(y)226 print *, "U con int(y)=", Ucalc ()227 print *, " "228

229 call cpu_time ( stop1 )230

231 print *, " tempo ", stop1 -start1 , "s; ", (stop1 - start1 )/60. , " min; ", &232 (stop1 - start1 )/3600. , " h"233

234

235 end program NR

45

46

Bibliografia

[1] E. Benedetti, M. Budinich, Neural Relax, Neural Computation, 24(11), 3091-3110, (2012)

[2] W.S. McCulloch, W.H. Pitts, A logical calculus of the ideas immanent innervous activity, Bulletin of Mathematical Biophysics, 5, 115-133, (1943)

[3] J. Atick, Could information theory provide an ecological theory of sensoryprocessing?, Network: Computation in Neural Systems, 3(2), 213-251, (1992)

[4] D.O. Hebb, The organization of behavior, New York: Wiley & Sons, (1949)[5] R.E. Schwartz, The five electron case of Thomson’s problem,

http://www.math.brown.edu/res/papers/electron.pdf[6] A.J. Bell e T.J. Sejnowski, An information-maximization approach to blind

separation and blind deconvolution, Neural Computation, 7(6), 1129-1159,(1995)

[7] S. Axler, P. Bourdon e W. Ramey, Harmonic function theory (2nd ed.), NewYork: Springer, (2001)

[8] M. Budinich, Introduzione alla teoria delle reti neurali, Università degli studidi Trieste, (2015)

[9] J.J. Thomson, On the Structure of the Atom: an Investigation of the Stailityand Periods of Oscillation of a number of Corpuscles arranged at equal intervalsaround the Circumference of a Circle; with Application of the Results to theTheory of Atomic Structure, Philosophical Magazine Series 6, 7, 39, 237-265,(1904)

47