Optim.jl：負の逆ヘッセ行列

Question

私はBFGSアルゴリズムを使ってJuliaの関数を最小限にするためにOptim.jlライブラリを使用しています。今日、私はquestionに同じライブラリについて尋ねましたが、混乱を避けるために2つに分割することに決めました。

私はまた、最適化後の負の逆ヘッセ行列の推定値を得て、さらなる計算をしたいと思います。 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()) 

にはどうすれば結果から最適化後の負の逆ヘッセ行列を得ることができますか？ヘッセ行列、逆ヘッセ行列、逆逆ヘッセ行列のいずれかを識別するフィールドがありますか？コメントの

EDIT

感謝。関数が "逆ヘッセ行列"を返すように "optimize.jl"を編集する方が効率的だと思いますか？作業例については以下を参照してください - 編集ライン226で導入されました：最適化の後に「状態」へのフルアクセスを持っているために

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 

：

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 

それとも。

EDIT 2

この変更はOptim.jlライブラリの新しいバージョンで導入されますので、この議論を継続する必要はありません。今のところ、extended_traceとafter_while！トリックが動作します。個人的には後者を好むので、@Dan Getzに正しい答えを与える議論を終わらせます。

Answer 1

現在では何もしていない内部オプティマイザ機能after_while!に接続していて、最後の状態から情報を引き出すのに最適ではありません。ジュリアコードで

これは次のようになります（多分v0.6で、これはすでに解決している）。これは、グローバル変数を使用するとoptimizeをコンパイル/実行する前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 

。

@DSMが彼の答えで述べたように、オプティマイザの最後の状態にアクセスしたいのは当然です。トレースが答えではないかもしれないなら、これは多分です。

Answer 2

私は一方的な方法を知っていますが、逆ヘッセ行列を自分で見積もるのではなく、自分では分かりません。 Optim.Options(store_trace=true, extended_trace=true)を渡すと、最後の〜invHを含む最適化パスのレコードを取得できます。

最初extended_trace=trueをオンにすると思われることである：例えば、

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 

を得ることができますshow_trace=trueを強制する - 私は計算の出力を表示していないことに気付くでしょう！バグのように感じる。これは、show_everyを大きな値に設定するか、stdoutを完全にリダイレクトすることで回避できます（ただし、削除されません）。

2つ目は、パス全体を保存せずに最後の状態にアクセスできるようにする必要がありますが、これが実際に問題になるかどうかは、問題のサイズによって異なります。

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で突然変異した状態から逆ヘッセ行列にアクセスすることができます。

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

そして、実際のヘッセ行列の逆数を計算して得られた逆ヘッセ行列と比較してください。

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実行の最後のステップで得られたものと同様です。