FSharp.Stats
Optimization
Summary: This tutorial teaches how to use optimization methods within FSharp.Stats
Nelder-Mead
The Nelder-Mead method (also downhill simplex method) can be used to find the minimum or maximum of an objective function. Please check out Mathias' blog post about the nelder mead algorithm.
Quadratic function
Task: Identify the minimum of the following function:
\[f(x)=x^2-0.32x-0.13\]
open FSharp.Stats
open FSharp.Stats.Optimization
open System
open Plotly.NET
open Plotly.NET.TraceObjects
let myFunction (xs: Vector<float>) =
let x = xs.[0]
x**2. + 0.32*x + 0.13
// initial guess for the optimization
let x0 = vector [| 0.95|]
// default solver options
let nmc = NelderMead.NmConfig.defaultInit()
// optimization procedure
let optim =
//let stopCrit =
// { OptimizationStop.defaultStopCriteria with MinFunctionEpsilon = 1e-24 }
//NelderMead.minimizeWithStopCriteria nmc x0 myFunction stopCrit
NelderMead.minimize nmc x0 myFunction
// optimization results as x, y, and z coordinate
let xs,ys =
optim.Vectors.[0..40] |> Array.map (fun x -> x.[0],myFunction x)
|> Array.unzip
let optimizationPathchart =
[
Chart.Line(x=xs,y=ys,ShowMarkers=true,Name="Optimization path")
[-1. .. 0.005 .. 1.] |> List.map (fun x -> x,myFunction (vector [x])) |> Chart.Line |> Chart.withTraceInfo "myFunction"
Chart.Point([|optim.SolutionVector.[0],optim.Solution|],Name="Solution") |> Chart.withMarkerStyle(Size=20,Symbol=StyleParam.MarkerSymbol.ArrowUp)
]
|> Chart.combine
|> Chart.withTemplate ChartTemplates.lightMirrored
|> Chart.withXAxisStyle ("x",ShowGrid=false)
|> Chart.withYAxisStyle ("myFunction(x)",ShowGrid=false)
|> Chart.withSize (800.,800.)
Rosenbrock
In the following chapter the minium of the three-dimensional rosenbrock function is identified.
\[f(x,y) = (\alpha - x)^2 + \beta(y-x^2)^2\]
When \(\alpha = 1\) and \(\beta = 100\) the minimum is at \(\alpha^2=1\).
Lets define the function, and a starting coordinate for the optimization task.
// Rosenbrock's valley or Rosenbrock's banana function
let rosenbrock (xs: Vector<float>) =
let x, y = xs.[0], xs.[1]
pown (1.0 - x) 2 + 100.0 * pown (y - pown x 2) 2
// initial guess for the optimization
let x0_rb = vector [| 1.85; -1.65|]
// rosenbrock visualization
let rosenBrockChart =
let range = [-2. .. 0.05 .. 2.]
range
|> List.map (fun y ->
range
|> List.map (fun x ->
rosenbrock (vector [x;y])
)
)
|> fun z ->
Chart.Surface(zData=z,X=range,Y=range,Opacity = 0.5, Contours = Contours.initXyz (Show = true))
let startConditionsChart =
[
Chart.Point3D([x0_rb.[0],x0_rb.[1],rosenbrock x0_rb])
rosenBrockChart
]
|> Chart.combine
|> Chart.withTemplate ChartTemplates.lightMirrored
|> Chart.withXAxisStyle ("x", Id = StyleParam.SubPlotId.Scene 1,ShowGrid=false)
|> Chart.withYAxisStyle ("y", Id = StyleParam.SubPlotId.Scene 1,ShowGrid=false)
|> Chart.withZAxisStyle ("rosenbrock(x,y)",ShowGrid=false)
|> Chart.withSize (800.,800.)
Now the functions minimum should be identified using the Nelder-Mead method. Default solver options are used for optimizations.
// default solver options
let nmc_rb = NelderMead.NmConfig.defaultInit()
// optimization procedure
let optim_rb = NelderMead.minimize nmc_rb x0_rb rosenbrock
// minimum x and y value
optim_rb.SolutionVector //vector [|0.9999978246; 1.000002057|]
// minimum z value
optim_rb.Solution //4.110573695e-09
// all z values during optimization steps
optim_rb.Values
// all x and y values during optimization steps
optim_rb.Vectors
Plotting of the optimization path
The minimum was correctly identified to be located at \((1,1)\). Lets investigate the path the Nelder-Mead method took to converge to this result.
// optimization results as x, y, and z coordinate
let x3d,y3d,z3d =
optim_rb.Vectors.[0..60] |> Array.map (fun x -> x.[0],x.[1],rosenbrock x)
|> Array.unzip3
let optimizationPathChart =
[
Chart.Line3D(x=x3d,y=y3d,z=z3d,ShowMarkers=true)
rosenBrockChart
]
|> Chart.combine
|> Chart.withTemplate ChartTemplates.lightMirrored
|> Chart.withXAxisStyle ("x", Id = StyleParam.SubPlotId.Scene 1,ShowGrid=false)
|> Chart.withYAxisStyle ("y", Id = StyleParam.SubPlotId.Scene 1,ShowGrid=false)
|> Chart.withZAxisStyle ("rosenbrock(x,y)",ShowGrid=false)
|> Chart.withSize (800.,800.)
Auckley function
The Auckley function has many valleys, with one center and global minimum at \((0,0)\).
// Auckley function
let auckley (xs: Vector<float>) =
let x, y = xs.[0], xs.[1]
-20.*exp(-0.2*sqrt(0.5*(x**2. + y**2))) -
exp(0.5*(cos(2. * Math.PI * x) + cos(2. * Math.PI * y))) +
Math.E + 20.
// auckley visualization
let auckleyChart =
let range = [-2. .. 0.05 .. 2.]
range
|> List.map (fun y ->
range
|> List.map (fun x ->
auckley (vector [x;y])
)
)
|> fun z ->
Chart.Surface(zData=z,X=range,Y=range,Opacity = 0.5)
// initial guess for the optimization
let x0_auckley = vector [| 1.2; -1.25 |]
// default solver options
let nmc_auckley = NelderMead.NmConfig.defaultInit()
// optimization procedure
let optim_auckley = NelderMead.minimize nmc_auckley x0_auckley auckley
// optimization results as x, y, and z coordinate
let x3d_auc,y3d_auc,z3d_auc =
optim_auckley.Vectors.[0..44] |> Array.map (fun x -> x.[0],x.[1],auckley x)
|> Array.unzip3
let optimizationPathChart_auckley =
[
Chart.Line3D(x=x3d_auc,y=y3d_auc,z=z3d_auc,ShowMarkers=true)
auckleyChart
]
|> Chart.combine
|> Chart.withTemplate ChartTemplates.lightMirrored
|> Chart.withXAxisStyle ("x", Id = StyleParam.SubPlotId.Scene 1,ShowGrid=false)
|> Chart.withYAxisStyle ("y", Id = StyleParam.SubPlotId.Scene 1,ShowGrid=false)
|> Chart.withZAxisStyle ("auckley(x,y)",ShowGrid=false)
|> Chart.withSize (800.,800.)
The Nelder-Mead method is able to identiy a local minimum, but misses the global minimum
Beale function
// Beale function function
let beale (xs: Vector<float>) =
let x, y = xs.[0], xs.[1]
(1.5 - x + x*y)**2. +
(2.25 - x + x*y**2)**2. +
(2.625 - x + x*y**3)**2.
// Beale visualization
let bealeChart =
let range = [-4. .. 0.01 .. 4.]
range
|> List.map (fun y ->
range
|> List.map (fun x ->
log10 (beale (vector [x;y]))
)
)
|> fun z ->
Chart.Surface(zData=z,X=range,Y=range,Opacity = 0.5)
// initial guess for the optimization
let x0_beale_1 = vector [| -3.5; -3.5 |]
let x0_beale_2 = vector [| 3. ; 2.5 |]
// default solver options
let nmc_beale = NelderMead.NmConfig.defaultInit()
// optimization procedure
let optim_beale_1 = NelderMead.minimize nmc_beale x0_beale_1 beale
let optim_beale_2 = NelderMead.minimize nmc_beale x0_beale_2 beale
// optimization results as x, y, and z coordinate
let x3d_bea_1,y3d_bea_1,z3d_bea_1 =
optim_beale_1.Vectors.[0..34] |> Array.map (fun x -> x.[0],x.[1],beale x |> log10)
|> Array.unzip3
let x3d_bea_2,y3d_bea_2,z3d_bea_2 =
optim_beale_2.Vectors.[0..34] |> Array.map (fun x -> x.[0],x.[1],beale x |> log10)
|> Array.unzip3
let optimizationPathChart_beale =
[
bealeChart
Chart.Line3D(x=x3d_bea_1,y=y3d_bea_1,z=z3d_bea_1,ShowMarkers=true,Name="local minimum")
Chart.Line3D(x=x3d_bea_2,y=y3d_bea_2,z=z3d_bea_2,ShowMarkers=true,Name="global minimum")
]
|> Chart.combine
|> Chart.withTemplate ChartTemplates.lightMirrored
|> Chart.withXAxisStyle ("x", Id = StyleParam.SubPlotId.Scene 1,ShowGrid=false)
|> Chart.withYAxisStyle ("y", Id = StyleParam.SubPlotId.Scene 1,ShowGrid=false)
|> Chart.withZAxisStyle ("log10(beale(x,y))",ShowGrid=false)
|> Chart.withSize (800.,800.)