# esempio gradient descent:
# -------------------------------------------------------------
# 

# LISTATO ############################### (non eseguire in R)

1  # definiamo l'area di disegno
2  op <- par(pty = "s", mar = c(2, 1, 2, 0), cex.axis = 0.75,
    cex.main = 0.85)
3
4  # 1. generiamo i vettori delle variabili x1 ed x2 (i valori coprono i
    medesimi
    intervalli con lo stesso valore del passo) 
5  x2 <- x1 <- seq(-2 * pi, 2 * pi, 0.1)
6
7  # 2. definiamo la funzione da minimizzare
8  f<-function(x, y) x^2 + 3 * sin(y)
9 
10 # contour plot della funzione
11 contour(x1, x2, outer(x1, x2, f), col = "blue", main = "")
12
13 # titolo 
14 mtext(side = 3, paste("minimizzare : ", expression(z == x[1]^2 + 3 *
   sin(x[2]))))
15
16
17 # scegliamo un punto con il mouse (potete sceglierne 3 ... non tutti in
    una volta
    ...)
18 for (i in 1:3) {
19 
20 ######### GRADIENT DESCENT (DISCESA A GRADIENTE)
##########################################
21 # Razionale:
22 # La discesa a gradiente è basata sull'osservazione che se una
    funzione
23 # f(x) a valori reali è definita e differenziabile in un punto a,
    allora
24 # f *decresce più velocemente* se ci muoviamo da a nella direzione del 
25 # gradiente *NEGATIVO* di f in corrispondenza del punto a.
26 #
27 # Ne consegue che:
28 # se b = a - s * gradiente(f(a)) per s > 0 e sufficientemente piccolo 
29 # allora f(a) >= f(b)
30 # NBB: notare che usiamo il NEGATIVO del gradiente : - s *
    gradiente(f(a)) ( così
    otteniamo una discesa e non un'ascesa) 
31 #
32 # si parte quindi da un punto che presumiamo essere vicino ad un punto
     di minimo
33 # locale, e consideriamo la sequenza:
34 # a0 (punto di partenza), a1, a2, a3, ... tale che :
35 # an+1 = an - b * gradient(f(an)), con n >= 0 (n è l' "indice" del
    passo della
    ricerca)
36 #
37 # Abbiamo : f(ao) >= f(a1) >= f(a3) >= f(a3) >= ...
38 # e quindi -> SPERIAMO! <- di trovare un minimo locale (... non c'è
    garanzia di
    riuscita)
39 # NB: notate che non è stato detto niente sulla funzione gradiente() . 
40
41
42 # discesa a gradiente ( gradient descent o steepest descent) 
43    # un click del mouse viene letto come 2 valori nel piano (x ed y) e
      ci fornisce un "punto di partenza"
44 # per la ricerca
45
46    x = unlist(locator(1))    # NB: restituisce 2 valori (ricerca in 2D
    ... 2 variabili)
47
48
49    s = 0.05 
50    newx = x - s * c(2 * x[1], 3 * cos(x[2]))      # NB: anche qui 2
      variabili ( an+1 )
51    eps = 0.001                                    # eps serve per
      definire il criterio di stop
52    gap = abs(f(newx[1], newx[2]) - f(x[1], x[2])) # gap: f(an) -
      f(an+1)
53
54    # RICERCA: qui definiamo un criterio di stop (altrimenti si continua
      all'infinito). Continuiamo a cercare
55    #          fino a quando gap > eps 
56    while (gap > eps) {
57
58    #disegniamo una freccia da f(an) a f(an+1)
59        arrows(x[1], x[2], newx[1], newx[2], length = 0.075,
                col = "red")
60
61    # settiamo an+1 come il nuovo punto di "partenza" (nel
		prossimo passo di ricerca sara' an+1 a rivestire il ruolo
62    # che è stato svolto da an
63        x = newx
64
65    # utilizziamo il nuovo punto di partenza (x) per
		approssimare il gradiente
66        newx = x - s * c(2 * x[1], 3 * cos(x[2]))
67
68    # calcoliamo il gap tra f(an corrente) e f(an+1
	corrente)
69        gap = abs(f(newx[1], newx[2]) - f(x[1], x[2]))
70
71    # chiediamo ad R di aspettare un tempo pari a gap/3 
72        Sys.sleep(gap/3)
73
74    # abbiamo fatto l'update di x e di gap ... se gap è >
		di eps ripetiamo il ciclo
75    }
76 }
77
78 par(op)


############ versione senza commenti ... #####################################


op <- par(pty = "s", mar = c(2, 1, 2, 0), cex.axis = 0.75,cex.main = 0.85)
x2 <- x1 <- seq(-2 * pi, 2 * pi, 0.1)
f <- function(x, y) x^2 + 3 * sin(y)
contour(x1, x2, outer(x1, x2, f), col = "blue", main = "")
mtext(side = 3, expression(z == x[1]^2 + 3 * sin(x[2])))

for (i in 1:3) {
    x = unlist(locator(1))
    s = 0.05
    newx = x - s * c(2 * x[1], 3 * cos(x[2]))
    eps = 0.001
    gap = abs(f(newx[1], newx[2]) -f(x[1], x[2]))

          while (gap > eps) {
	          arrows(x[1], x[2], newx[1], newx[2], length = 0.075,
		  col = "red")
		  x = newx
		  newx = x - s * c(2 * x[1], 3 *cos(x[2]))
		  gap = abs(f(newx[1],newx[2]) - f(x[1],x[2]))
		  Sys.sleep(gap/3)
	  }

}

par(op)