Gradientoptimering med nloptr i R

Dette indlæg introducerer gradientnedstigningsoptimering i R ved hjælp af nloptr-pakken.

Til løsning af transportproblemer eller netværksmodelleringsproblemer er lineær programmering ofte tilstrækkelig.

Ikke-lineær programmering kan ikke desto mindre blive relevant når man overvejer yderligere begrænsninger eller mål som ikke er lineære.

R giver en pakke til løsning af ikke-lineære problemer: nloptr

I dette indlæg vil jeg anvende nloptr-pakken til at løse nedenstående ikke-lineære optimeringsproblem ved anvendelse af gradientafstamningsmetode. Gradientafstamningsalgoritmer ser efter retningen for den stejleste ændring, dvs. retningen for maksimum eller minimum første afledte. Logikken er, at hvis jeg fortsætter med at bevæge mig i retning af den stejleste nedstigning, vil jeg hurtigt bevæge mig mod mit optimale.

Problemet ved hånden:

Vigtig note: Når man bruger nloptr til at modellere ikke-lineære problemer skal man angive problemet som et minimeringsproblem, og alle begrænsninger i ligestilling skal være af typen <= 0 (mindre end). Funktionen kan håndtere ligestillingsbegrænsninger, så du ikke behøver at erstatte dem ved at overlappe begrænsninger i ligestilling.

Gradientoptimering kræver, at man vælger et startpunkt for optimeringen. Dette startpunkt skal altså defineres af analysten.

Jeg modellerer og løser dette problem med gradientnedstigning, ved hjælp af nloptr-pakken i R (med nloptr-funktionen!):

# importer nloptr pakken
library(nloptr)

# definer funktionen som skal optimeres 
eval_f <- function(x){
  return(
    list(
      "objective"=x[1]^2+x[2]^2,
      "gradient"=c(2*x[1],
                   2*x[2])
       )
  )
}

# definer funktionen som skal repræsentere restriktioner i form af uligheder
eval_g_ineq <- function(x){
  return(
    list(
      "constraints"=x[1]+x[2]-100,
      "jacobian"=c(1,
                   1)
    )
  )
}

# definer startpunkt for optimeringen
x_0 <- c(10,10)

# yderligere indstillinger som anbefalet i nloptr dokumentationen
local_opts <- list( "algorithm" = "NLOPT_LD_MMA",
                    "xtol_rel"  = 1.0e-7 )

opts <- list( "algorithm" = "NLOPT_LD_AUGLAG",
              "xtol_rel"  = 1.0e-7,
              "maxeval"   = 1000,
              "local_opts" = local_opts )

# modeller og løs problemet
solution <- nloptr(x0=x_0,
                   eval_f=eval_f,
                   lb=NULL,
                   ub=NULL,
                   eval_g_ineq = eval_g_ineq,
                   eval_g_eq = NULL,
                   opts=opts)

Vi kan nu tage et nærmere kig på resultatet af optimeringen:

print(solution)

## 
## Call:
## 
## nloptr(x0 = x_0, eval_f = eval_f, lb = NULL, ub = NULL, eval_g_ineq = eval_g_ineq, 
##     eval_g_eq = NULL, opts = opts)
## 
## 
## Minimization using NLopt version 2.4.2 
## 
## NLopt solver status: 3 ( NLOPT_FTOL_REACHED: Optimization stopped because 
## ftol_rel or ftol_abs (above) was reached. )
## 
## Number of Iterations....: 102 
## Termination conditions:  xtol_rel: 1e-07 maxeval: 1000 
## Number of inequality constraints:  1 
## Number of equality constraints:    0 
## Optimal value of objective function:  1.09462903864043e-94 
## Optimal value of controls: -7.398071e-48 -7.398071e-48

Der findes også en god dokumentation om gradientoptimering på Wikipedia: https://en.wikipedia.org/wiki/Gradient_descent