Optim.jl：负反Hessian

Question

我正在使用Optim.jl库来最小化Julia中的函数，使用BFGS算法。今天，我问了一个question关于同一个库，但为了避免混淆，我决定将它分成两部分。

我还想在优化后得到负反Hessian的估计值，以便进一步计算。

在库中的Optim的GitHub的网站，我发现下面的工作例如：

using Optim 
rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2 
result  = optimize(rosenbrock, zeros(2), BFGS()) 

我怎样才能得到结果优化后的负逆海森？是否有任何字段可以识别Hessian，Inverse Hessian或负Hessian？

编辑

感谢您的意见。你认为编辑“optimize.jl”会使函数返回逆Hessian更有效吗？请参阅下面的工作示例 - 编辑已经在管线226被引入：

if state.method_string == "BFGS" 
     return MultivariateOptimizationResults(state.method_string, 
               initial_x, 
               f_increased ? state.x_previous : state.x, 
               f_increased ? state.f_x_previous : state.f_x, 
               iteration, 
               iteration == options.iterations, 
               x_converged, 
               options.x_tol, 
               f_converged, 
               options.f_tol, 
               g_converged, 
               options.g_tol, 
               f_increased, 
               tr, 
               state.f_calls, 
               state.g_calls, 
               state.h_calls), state.invH 
    else 
     return MultivariateOptimizationResults(state.method_string, 
               initial_x, 
               f_increased ? state.x_previous : state.x, 
               f_increased ? state.f_x_previous : state.f_x, 
               iteration, 
               iteration == options.iterations, 
               x_converged, 
               options.x_tol, 
               f_converged, 
               options.f_tol, 
               g_converged, 
               options.g_tol, 
               f_increased, 
               tr, 
               state.f_calls, 
               state.g_calls, 
               state.h_calls) 
    end 

或者只是：

return MultivariateOptimizationResults(state.method_string, 
             initial_x, 
             f_increased ? state.x_previous : state.x, 
             f_increased ? state.f_x_previous : state.f_x, 
             iteration, 
             iteration == options.iterations, 
             x_converged, 
             options.x_tol, 
             f_converged, 
             options.f_tol, 
             g_converged, 
             options.g_tol, 
             f_increased, 
             tr, 
             state.f_calls, 
             state.g_calls, 
             state.h_calls), state 

要使优化后全面进入“状态”。

EDIT 2

因为这种变化将在Optim.jl库的新版本推出，就没有必要继续讨论。就目前而言，extended_trace和after_while！技巧工作。就我个人而言，我更喜欢后者，所以我将结束讨论，给予丹茨茨正确的答案。

Answer 1

另一个不是最佳的方法是连接到内部Optim函数after_while!，该函数目前什么都不做，并用它从最后一个状态中提取信息。

在朱莉娅的代码，这看起来像：

julia> using Optim 

julia> rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2 
rosenbrock (generic function with 1 method) 

julia> Optim.after_while!{T}(d, state::Optim.BFGSState{T}, method::BFGS, options) 
    = global invH = state.invH 

julia> result  = optimize(rosenbrock, zeros(2), BFGS()) 
Results of Optimization Algorithm 
* Algorithm: BFGS 
* Starting Point: [0.0,0.0] 
* Minimizer: [0.9999999926033423,0.9999999852005353] 
* Minimum: 5.471433e-17 
* Iterations: 16 
* Convergence: true 
    * |x - x'| < 1.0e-32: false 
    * |f(x) - f(x')|/|f(x)| < 1.0e-32: false 
    * |g(x)| < 1.0e-08: true 
    * f(x) > f(x'): false 
    * Reached Maximum Number of Iterations: false 
* Objective Function Calls: 69 
* Gradient Calls: 69 

julia> invH 
2×2 Array{Float64,2}: 
0.498092 0.996422 
0.996422 1.9983 

这是没有吸引力的使用全局变量和因运行/编译optimize之前定义after_while!（但也许V0.6，这是已经解决） 。

正如@ DSM在他的回答中指出的那样，想要访问优化器的最后一个状态是很自然的。如果追踪不是答案，那可能就是这样。

Answer 2

我知道一种方法，但是否值得，而不是自己估计逆Hessian，我不确定。如果您通过Optim.Options(store_trace=true, extended_trace=true)，则可以获得包含最后一个invH的优化路径的记录。例如，在

result = optimize(rosenbrock, zeros(2), BFGS(), 
        Optim.Options(store_trace=true, extended_trace=true)); 

我们可以得到

julia> result.trace[end] 
    16  5.471433e-17  2.333740e-09 
* Current step size: 1.0091356566200325 
* g(x): [2.33374e-9,-1.22984e-9] 
* ~inv(H): [0.498092 0.996422; 0.996422 1.9983] 
* x: [1.0,1.0] 


julia> result.trace[end].metadata["~inv(H)"] 
2×2 Array{Float64,2}: 
0.498092 0.996422 
0.996422 1.9983 

至少有两件事情，我不喜欢这种做法，虽然：

首先是打开extended_trace=true似乎强制show_trace=true - 你会注意到我没有显示计算的输出！感觉像一个错误。这可以通过将show_every设置为很大的值来减轻（尽管不会被删除），或者通过完全重定向stdout来避免。

第二个是我们应该能够访问最后一个状态而不存储整个路径，但是否这实际上是一个问题将取决于问题的大小。

Answer 3

目前在Optim.jl中做到这一点的最简单的方法是执行以下操作。

首先，负载的Optim和OptimTestProblems（有一个例子断的工作）

julia> using Optim, OptimTestProblems 

julia> prob = OptimTestProblems.UnconstrainedProblems.examples["Himmelblau"] 
OptimTestProblems.MultivariateProblems.OptimizationProblem{Void,Void,Float64,String,Void}("Himmelblau", OptimTestProblems.MultivariateProblems.UnconstrainedProblems.himmelblau, OptimTestProblems.MultivariateProblems.UnconstrainedProblems.himmelblau_gradient!, nothing, OptimTestProblems.MultivariateProblems.UnconstrainedProblems.himmelblau_hessian!, nothing, [2.0, 2.0], [3.0, 2.0], 0.0, true, true, nothing) 

然后指定部分optimize需求都在正确的顺序输入：

julia> x0 = prob.initial_x 
2-element Array{Float64,1}: 
2.0 
2.0 

julia> obj = OnceDifferentiable(prob.f, prob.g!, x0) 
NLSolversBase.OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1},Val{false}}(OptimTestProblems.MultivariateProblems.UnconstrainedProblems.himmelblau, OptimTestProblems.MultivariateProblems.UnconstrainedProblems.himmelblau_gradient!, NLSolversBase.fg!, 0.0, [NaN, NaN], [NaN, NaN], [NaN, NaN], [0], [0]) 

julia> m = BFGS() 
Optim.BFGS{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64},Optim.##67#69}(LineSearches.InitialStatic{Float64} 
    alpha: Float64 1.0 
    scaled: Bool false 
, LineSearches.HagerZhang{Float64} 
    delta: Float64 0.1 
    sigma: Float64 0.9 
    alphamax: Float64 Inf 
    rho: Float64 5.0 
    epsilon: Float64 1.0e-6 
    gamma: Float64 0.66 
    linesearchmax: Int64 50 
    psi3: Float64 0.1 
    display: Int64 0 
, Optim.#67, Optim.Flat()) 

julia> options = Optim.Options() 
Optim.Options{Float64,Void}(1.0e-32, 1.0e-32, 1.0e-8, 0, 0, 0, false, 0, 1000, false, false, false, 1, nothing, NaN) 

julia> bfgsstate = Optim.initial_state(m, options, obj, x0) 
Optim.BFGSState{Array{Float64,1},Array{Float64,2},Float64,Array{Float64,1}}([2.0, 2.0], [6.91751e-310, 6.9175e-310], [-42.0, -18.0], NaN, [6.9175e-310, 0.0], [6.91751e-310, 0.0], [6.91751e-310, 0.0], [1.0 0.0; 0.0 1.0], [6.91751e-310, 0.0], NaN, [6.9175e-310, 6.91751e-310], 1.0, false, LineSearches.LineSearchResults{Float64}(Float64[], Float64[], Float64[], 0)) 

julia> res = optimize(obj, x0, m, options, bfgsstate) 
Results of Optimization Algorithm 
* Algorithm: BFGS 
* Starting Point: [2.0,2.0] 
* Minimizer: [2.9999999999998894,2.000000000000162] 
* Minimum: 5.406316e-25 
* Iterations: 7 
* Convergence: true 
    * |x - x'| ≤ 1.0e-32: false 
    |x - x'| = 5.81e-09 
    * |f(x) - f(x')| ≤ 1.0e-32 |f(x)|: false 
    |f(x) - f(x')| = 2.93e+09 |f(x)| 
    * |g(x)| ≤ 1.0e-08: true 
    |g(x)| = 4.95e-12 
    * Stopped by an increasing objective: false 
    * Reached Maximum Number of Iterations: false 
* Objective Calls: 42 
* Gradient Calls: 42 

然后我们可以从optimize中突变的状态访问逆Hessian。

julia> bfgsstate.invH 
2×2 Array{Float64,2}: 
    0.0160654 -0.00945561 
-0.00945561 0.034967 

并将其与通过计算实际Hessian的逆得到的逆Hessian进行比较。

julia> H=similar(bfgsstate.invH) 
2×2 Array{Float64,2}: 
6.91751e-310 6.91751e-310 
6.91751e-310 6.91751e-310 

julia> prob.h!(H, Optim.minimizer(res)) 
34.00000000000832 

julia> H 
2×2 Array{Float64,2}: 
74.0 20.0 
20.0 34.0 

julia> inv(H) 
2×2 Array{Float64,2}: 
    0.0160681 -0.0094518 
-0.0094518 0.0349716 

这与在BFGS运行的最后一步获得的相似。