I’m having some issues with optimizing an objective with Newton from Optim. I started a brief discussion on slack which has led me to isolate the problem:
Okay I am getting closer to the issue. So consider gradient and hessian:
st0=[160.2, -54.8, -210.8, -0.8, 106.2, -0.0];
st1=[137.6 -34.4 -34.4 -34.4 -34.4 -1.0; -34.4 137.6 -34.4 -34.4 -34.4 -1.0; -34.4 -34.4 137.6 -34.4 -34.4 -1.0; -34.4 -34.4 -34.4 137.6 -34.4 -1.0; -34.4 -34.4 -34.4 -34.4 137.6 -1.0; -1.0 -1.0 -1.0 -1.0 -1.0 0.0]
the newton step I am expecting to happen is
Main> st1\st0
6-element Array{Float64,1}:
0.931395
-0.318605
-1.22558
-0.00465116
0.617442
4.90803e-14
which (if you do two iterations of this) you get close to my optimum. The Optim package is using PositiveFactorizations approach (because st1 is not PosDef) and does (basically) this:
Main> F=cholfact(Positive, LowerTriangular(st1))
Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:
[11.7303 0.0 … 0.0 0.0; -2.93258 11.3578 … 0.0 0.0; … ; -2.93258 -3.78594 … 1.0 0.0; -0.0852493 -0.110056 … 0.0 0.340997]
Main> F\st0
6-element Array{Float64,1}:
-0.457849
-1.70785
-2.61483
-1.3939
-4.12115e-13
-26.55
(which corresponds to the values I see in the trace of the Newton output). This does not converge to the right solution. What is the problem here?
These two approaches boil down to small differences in the hessian that is used in the Newton step:
Main> F.factors*ctranspose(F.factors)
6×6 Array{Float64,2}:
137.6 -34.4 -34.4 -34.4 -34.4 -1.0
-34.4 137.6 -34.4 -34.4 -34.4 -1.0
-34.4 -34.4 137.6 -34.4 -34.4 -1.0
-34.4 -34.4 -34.4 137.6 -34.4 -1.0
-34.4 -34.4 -34.4 -34.4 138.6 4.0
-1.0 -1.0 -1.0 -1.0 4.0 0.232558
Main> st1
6×6 Array{Float64,2}:
137.6 -34.4 -34.4 -34.4 -34.4 -1.0
-34.4 137.6 -34.4 -34.4 -34.4 -1.0
-34.4 -34.4 137.6 -34.4 -34.4 -1.0
-34.4 -34.4 -34.4 137.6 -34.4 -1.0
-34.4 -34.4 -34.4 -34.4 137.6 -1.0
-1.0 -1.0 -1.0 -1.0 -1.0 0.0
Below is 1 iteration using Optim.Newton (showing values that correspond to the PositiveFactorizations)
Main> res=Optim.optimize(logLC, g!, h!, zeros(m+1), show_trace=true, extended_trace=true, method=Optim.Newton(linesearch=LineSearches.Static(alpha=1.0)), iterations=1)
Iter Function value Gradient norm
0 1.348190e+03 2.108000e+02
* Current step size: 1.0
* g(x): [160.2, -54.8, -210.8, -0.8, 106.2, -0.0]
* h(x): [137.6 -34.4 -34.4 -34.4 -34.4 -1.0; -34.4 137.6 -34.4 -34.4 -34.4 -1.0; -34.4 -34.4 137.6 -34.4 -34.4 -1.0; -34.4 -34.4 -34.4 137.6 -34.4 -1.0; -34.4 -34.4 -34.4 -34.4 137.6 -1.0; -1.0 -1.0 -1.0 -1.0 -1.0 0.0]
* x: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
1 9.554691e+02 7.972837e+01
* Current step size: 1.0
* g(x): [10.9235, -30.7894, -2.10531, -31.0505, -79.7284, -6.17442]
* h(x): [79.4839 -23.51 -23.3 -21.8481 -10.8258 -1.0; -23.51 135.839 -54.0743 -41.7406 -16.5137 -1.0; -23.3 -54.0743 138.641 -45.3491 -15.9179 -1.0; -21.8481 -41.7406 -45.3491 124.487 -15.5493 -1.0; -10.8258 -16.5137 -15.9179 -15.5493 58.8067 -1.0; -1.0 -1.0 -1.0 -1.0 -1.0 0.0]
* x: [0.457849, 1.70785, 2.61483, 1.3939, 4.12115e-13, 26.55]
Results of Optimization Algorithm
* Algorithm: Newton's Method
* Starting Point: [0.0,0.0,0.0,0.0,0.0,0.0]
* Minimizer: [0.45784883720971925,1.7078488372097196, ...]
* Minimum: 9.554691e+02
* Iterations: 1
* Convergence: false
* |x - x'| ≤ 1.0e-32: false
|x - x'| = 2.66e+01
* |f(x) - f(x')| ≤ 1.0e-32 |f(x)|: false
|f(x) - f(x')| = 4.11e-01 |f(x)|
* |g(x)| ≤ 1.0e-08: false
|g(x)| = 7.97e+01
* Stopped by an increasing objective: false
* Reached Maximum Number of Iterations: true
* Objective Calls: 2
* Gradient Calls: 2
* Hessian Calls: 2
