testo

Laboratorio 5
Ottimizzazione
Prof. E. Amaldi
Problema della distanza molecolare e addestramento di reti
neurali
Esercizio 1: geometria molecolare con distanze euclidee
Il problema noto in letteratura come Molecular Distance Geometry Problem, o MDGP, consiste
nel ricavare la forma tri-dimensionale di una molecola in base a informazioni relative alle distanze
tra le coppie di atomi che la compongono. Consideriamo una molecola costituita da m atomi.
Supponiamo che siano note (con precisione infinita) le distanze tra tutte le coppie di atomi. Sia
dunque D ∈ Rm×m la matrice (simmetrica) delle distanze euclidee dij per ogni coppia {i, j}.1
Siano x1 , . . . , xm ∈ R3 le posizioni incognite degli atomi. Senza perdita di generalit`a, possiamo supporre che il sistema di coordinate sia tale che xm = (0, 0, 0) (eventualmente a valle di una
traslazione). Il problema consiste nel determinare le posizioni (coordinate) x1 , . . . , xm−1 ∈ R3
degli m − 1 altri atomi cos`ı da rispettare le distanze dij tra tutte le coppie i e j. I vincoli
corrispondenti sono:
||xi − xj ||2 = dij , i, j = 0, 1, . . . , m.
a) Si formuli il problema in termini di ottimizzazione non lineare non vincolata.
b) Data la non convessit`
a del problema di ottimizzazione ottenuto, implementare un metodo
Multistart basato su uno dei due metodi di ottimizzazione non lineare implementati nel
laboratorio precedente. Siano:
• ε > 0 la tolleranza dell’algoritmo Multistart
• f¯ il valore ottimo esatto (noto) della funzione obiettivo
• M il numero massimo di iterazioni del metodo Multistart.
L’algoritmo Multistart procede come segue:
(a) Sia x un qualsiasi punto in R3(m−1) . Sia x∗ ← x.
(b) Se f (x∗ ) < f¯ + ε oppure se sono state fatte pi`
u di M iterazioni, l’algoritmo termina;
altrimenti, si va al passo c).
(c) Trovare un minimo locale x0 = (x1 , . . . , xm−1 ) a partire dal punto iniziale x, con un
metodo di ottimizzazione nonlineare (il quale, a sua volta, avr`a una tolleranza ε0 > 0
e un numero massimo di iterazioni M 0 ).
(d) Se f (x0 ) < f (x∗ ) aggiornare x∗ ← x0 .
(e) Scegliere un nuovo punto iniziale x ∈ R3(m−1) , a caso.
(f) Tornare al passo b).
c) Applicare l’algoritmo alla seguente istanza con m = 4 atomi (con dii = 0 per ogni i):
1
Un riferimento bibliografico per questo esercizio `e: J.M. Yoon, Y. Gad, Z. Wu, Mathematical modeling of
protein structure using distance geometry, Technical report TR00-24, DCAM, Rice University, Houston, 2000,
scaricabile da: http://www.caam.rice.edu/caam/trs/tr00.html#TR00-24.
Edizione 2014 a cura di L. Taccari. Edizione originale di L. Liberti.
1
Laboratorio 5
Prof. E. Amaldi
Ottimizzazione
dij
1
2
3
2
1.526
0
-
3
2.491389536
1.526
0
4
3.837572036
2.491389535
1.526
Per la costruzione della soluzione aleatoria iniziale, si usi la funzione rnd.m, che genera un
vettore casuale di n elementi compresi tra -bound e bound, secondo una distribuzione di
probabilit`
a uniforme.
function xrnd = rnd(n, bound)
xrnd = zeros(n,1);
for i = 1:size(xrnd,1)
xrnd(i) = (rand()-0.5)*2*bound;
end
end %of function
d) Si applichi il Metodo di Gauss-Newton per la soluzione di problemi di minimi quadrati
non lineari. A tale scopo, forniamo il file jac.m, che calcola lo Jacobiano nel punto x di
una funzione vettoriale r con m componenti. Il secondo parametro rappresenta il numero
di componenti della funzione.
% Jacobian of a vector function at a point
function J = jac(r, m, x)
n = length(x);
J = zeros(m,n);
h = 0.0001;
for i = 1:m
for j = 1:n
delta = zeros(n, 1); delta(j) = h;
rd=r(x+delta);
rx=r(x);
J(i,j) = (rd(i) - rx(i)) / h;
end
end
end %end of function
Si utilizzi inoltre la funzione implementata in Matlab pinv, che, data una matrice A, ne
−1
calcola la pseudoinversa AT A AT .
Il metodo di Gauss-Newton Si tratta di una variante del metodo di Newton per
risolvere problemi di minimi quadrati non lineari. Consideriamo il seguente problema
generale:
f (x) =
m
X
(ri (x))2 ,
(1)
i=1
dove x ∈ Rn ed r(x) = (r1 (x), . . . , rm (x))T `e il vettore dei residui. Assumiamo che tutte
le funzioni ri (x) siano non lineari.
Derivando (1) si ottiene
∇x f (x) =
m
X
2ri (x)∇x ri (x).
i=1
Edizione 2014 a cura di L. Taccari. Edizione originale di L. Liberti.
2
Laboratorio 5
Prof. E. Amaldi
Ottimizzazione
n
∂ri
∂xk
o
lo Jacobiano di r valutato nel punto x: per brevit`a, indichiamo
 T

∇x r1 (x)


..
. Possiamo quindi riscrivere
Jx r(x) con J (x). Si osservi che J (x) = 
.


T
∇x rm (x)
l’espressione come
Sia Jx r(x) =
ik
∇f (x) = 2Jx r(x)T r(x).
La matrice Hessiana di f valutata nel punto x, che per brevit`a indichiamo con H (x), `e
quindi
T
H (x) = 2J (x) J (x) + 2
m
X
(ri (x))∇2 ri (x).
(2)
i=1
Se i residui sono piccoli, si pu`
o trascurare l’ultimo termine di (2), che porta a
H (x) ≈ 2J (x)T J (x) .
Si noti che, se i residui sono funzioni lineari di x, allora ∇2 ri (x) = 0 per ogni i e l’approssimazione al secondo ordine di Taylor `e esatta. In questo caso il metodo viene a coincidere
con quello dei minimi quadrati lineari.
Il passo del metodo di Newton
xk+1 = xk − [H (xk )]−1 ∇f (xk )
corrisponde dunque, usando l’approssimazione dell’Hessiano prima derivata, a:
h
i−1
xk+1 = xk − J (xk )T J (xk )
J (xk )T r(xk ).
Affinch´e il metodo di Gauss-Newton converga, la soluzione iniziale deve essere sufficientemente vicina ad un punto stazionario di f e i termini trascurati dell’Hessiano (2) devono
essere piccoli. Si noti che il metodo non richiede richiede il calcolo dell’Hessiano H (x) ad
ogni iterazione, che implicherebbe il calcolo degli Hessiani ∇2 ri per ogni ri .
Figura 1: Configurazione ottimale per l’istanza del problema.
Edizione 2014 a cura di L. Taccari. Edizione originale di L. Liberti.
3
Laboratorio 5
Prof. E. Amaldi
Ottimizzazione
Esercizio 2 (extra): addestramento di reti neurali
A una piccola ditta di assemblaggio di circuiti integrati `e stato spedito un chip logico con due
piedini di input e uno di output, di cui purtroppo si `e persa la bolla di accompagnamento. Non
si sa n´e da dove provenga, n´e che funzione abbia. Un tecnico di laboratorio conduce alcuni
esperimenti applicando agli ingressi le tensioni y1 , y2 riportate nella Tabella 1, e registra le
tensioni osservate in uscita z: Raccolti i dati, il tecnico decide di utilizzarli per addestrare una
Tabella 1: Dati dell’esercizio 2.
y1
y2
z
0.1
0.1
0.1
0.9
0.9
0.05
0.0
0.95 0.98
0.95
0.1
0.95
rete neurale che modellizzi la funzione logica implementata nel chip.
ξ
−ξ
3
a) Sia ϕ(ξ) = eeξ −e
la funzione di attivazione, il cui andamento `e rappresentato nella Figu+e−ξ
ra 2b. Si addestri la rete neurale rappresentata nella Figura 2a mediante il set di dati
(training set) sperimentali precedentemente riportati, individuando un insieme di pesi
u associati agli archi della rete tali da ben approssimare i dati. Si utilizzi l’algoritmo
Multistart per trovare un buon ottimo locale.
0
1
2
ϕ(ξ)
−2
−1
ξ
−3
(a) Schema della rete neurale considerata.
−2
−1
0
1
2
3
(b) Grafico della funzione ϕ (ξ).
Figura 2: Caratteristiche della rete neurale dell’esercizio 2.
Si usino il file f hmap.m, che rappresenta la funzione implementata dalla rete neurale in
Edizione 2014 a cura di L. Taccari. Edizione originale di L. Liberti.
4
Laboratorio 5
Ottimizzazione
Prof. E. Amaldi
Figura 2a:
function f = f_hmap(u, y)
f = f_act(u(5)* f_act(u(1)*y(1) + u(2)*y(2)) ...
+ u(6)* f_act(u(3)*y(1) + u(4)*y(2)));
end
function f = f_act(xi)
f = (exp(xi) - exp(-xi)) / (exp(xi) + exp(-xi));
end
e il file f nnet.m, che codifica la funzione d’errore da minimizzare:
function f = f_nnet(u)
y = [ 0.1 0.1 ;
0.9 0.9 ;
0.0 0.95 ;
0.95 0.1 ];
z = [ 0.1 0.05 0.98 0.95 ]’;
f = 0;
for i = 1:length(z)
f = f + ( z(i) - f_hmap(u, [y(i,1) y(i,2)]’) )^2;
end
end
b) Si calcolino le uscite della rete con le coppie di input (0,0), (0,1), (1,0), (1,1) e si stabilisca
che tipo di chip `e stato consegnato alla ditta.
Edizione 2014 a cura di L. Taccari. Edizione originale di L. Liberti.
5

Download Report