Help comparing KNITRO.jl and Matlab results

[pierremabille]

Hi everyone, my problem is the following. I am solving a nonlinear system of (14) equations without giving knitro the derivatives. Results in Matlab and Julia are inconsistent, could you please help me understand why?

In Matlab I use the knitromatlab_fsolve command, with some options (options=optimset('Display','iter-detailed','TolX',1e-15,'TolFun',1e-12,'MaxIter',100,'MaxFunEvals',20000)). I also print the detailed output of Matlab's resolution of the problem to know exactly which Knitro options it is using (available if someone wants to see them: algorithm, gradopt, hessopt, etc.). It turns out that the underlying function solves a constrained minimization problem in which the objective function is a constant and the constraints are the equations of the system I want to solve. The Matlab code solves the problem without any trouble and returns a solution vector different from my initial guess.

Now in Julia, I set up the latter unconstrained, derivative-free minimization problem, with the exact same options. However, the Julia solver seems to stay stuck at the initial guess. When I compare the verbose output with Matlab's, iteration 0 is the same, but iteration 1 in Matlab has a positive step, more fcounts than Julia, and a lower feaserror -- while Julia has a 0 step, the same feaserror as at the initial guess (iteration 0), and less fcounts. As a result Julia returns the same solution as the initial guess, different from Matlab, and tells me that the problem is locally infeasible (which I know is false since I solve what is supposed to be the exact same problem in Matlab).

Could anyone please help me make the Julia code consistent with Matlab? Here is the way I set up the Julia problem, is there something I am doing wrong?

`function eval_f(x::Vector{Float64}) # objective function, here constant 0 end

function eval_g(x::Vector{Float64}, cons::Vector{Float64}) # constraints function
    temp = solveAtPoint(mobj, x, gvec, point, 0, trans)[1]; # produces vector residuals we want to set to 0
    cons[1] = temp[1]
    cons[2] = temp[2]
    cons[3] = temp[3]
    cons[4] = temp[4]
    cons[5] = temp[5]
    cons[6] = temp[6]
    cons[7] = temp[7]
    cons[8] = temp[8]
    cons[9] = temp[9]
    cons[10] = temp[10]
    cons[11] = temp[11]
    cons[12] = temp[12]
    cons[13] = temp[13]
    cons[14] = temp[14]
    # return cons
end

function eval_grad_f(x::Vector{Float64}, grad::Vector{Float64}) # don't supply any of the derivatives
end

function eval_jac_g(x::Vector{Float64}, jac::Vector{Float64})
end

function eval_h(x::Vector{Float64}, lambda::Vector{Float64},
                sigma::Float64, hess::Vector{Float64})
end

function eval_hv(x::Vector{Float64}, lambda::Vector{Float64},
                 sigma::Float64, hv::Vector{Float64})
end

objType = KTR_OBJTYPE_GENERAL
objGoal = KTR_OBJGOAL_MINIMIZE

n=14 # number variables
x_L = -KTR_INFBOUND*ones(14)
x_U = KTR_INFBOUND*ones(14)

m=14 # number constraints
c_Type = [KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL,
KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL]
c_L = zeros(14)
c_U = zeros(14)

#jacIndexCons[i] and jacIndexVars[i] determine the row numbers and the column numbers, respectively, of the nonzero constraint Jacobian element jac[i].
jac_var = vec(repmat(Int32[0,1,2,3,4,5,6,7,8,9,10,11,12,13],14,1))
jac_con = vcat(round(Int32,ones(14)*0),round(Int32,ones(14)*1),round(Int32,ones(14)*2),round(Int32,ones(14)*3),round(Int32,ones(14)*4),
round(Int32,ones(14)*5),round(Int32,ones(14)*6),round(Int32,ones(14)*7),round(Int32,ones(14)*8),round(Int32,ones(14)*9),
round(Int32,ones(14)*10),round(Int32,ones(14)*11),round(Int32,ones(14)*12),round(Int32,ones(14)*13))

hess_row = Int32[]
hess_col = Int32[]

x = gvec # initial guess
lambda = zeros(n+m)
obj = 0

kp = createProblem()
loadOptionsFile(kp, "C:/some_path/knitro3.opt") # path of option file. change solver options in the latter
initializeProblem(kp, objGoal, objType, x_L, x_U, c_Type, c_L, c_U, jac_var, jac_con, hess_row, hess_col, x, lambda)
setCallbacks(kp, eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h, eval_hv)
solveProblem(kp)`

[pierremabille]

Issue solved. I can post a small example if anyone is interested.

[mlubin]

@pierremabille, you're more than welcome to add an example to the documentation, e.g., here: https://github.com/JuliaOpt/KNITRO.jl/blob/master/doc/example.rst

[yeesian]

Now that we're in the process of switching the docs from readthedocs to gh-pages, can you post the example here instead? I'll update the docs to include it.

[pierremabille]

Yes I will do it as soon as possible!

[pierremabille]

Ok, so here is a short example of how to solve a nonlinear system of equations in KNITRO.jl, without specifying any of the derivatives, when there is no analytical expression for the nonlinear system of equations to solve.

The setting and options replicate those of the Matlab function knitromatlab_fsolve, the knitro equivalent of Matlab's fsolve. Of course in my example those derivatives (gradient, constraints jacobian, hessian, etc.) can be computed analytically (because the nonlinear has an analytical expression here), but I didn't do it to show that the method works for any system for equations for which we have no analytical expression. The code is documented, and I attach the corresponding option file knitro.opt.

using KNITRO # load package

#function defining nonlinear system of equations (can have no analytic specification)
function nleq(x::Vector{Float64}) 
    eq1 = x[1]^2 + x[2]^2 - 25
    eq2 = (x[1]-8)^2 + x[2]^2 - 41
    return [eq1,eq2]
end

function eval_f(x::Vector{Float64}) #objective function, here constant
    0
end

function eval_g(x::Vector{Float64}, cons::Vector{Float64}) #constraints function, calling nleq
    temp = nleq(x)
    cons[1] = temp[1]
    cons[2] = temp[2]
end

function eval_grad_f(x::Vector{Float64}, grad::Vector{Float64}) #gradient evaluation, empty
end

function eval_jac_g(x::Vector{Float64}, jac::Vector{Float64}) #constraints jacobian evaluation, empty
end

function eval_h(x::Vector{Float64}, lambda::Vector{Float64}, #hessian evaluation, empty
                sigma::Float64, hess::Vector{Float64})
end

function eval_hv(x::Vector{Float64}, lambda::Vector{Float64},  #hessian of the Lagrangian vector product
                 sigma::Float64, hv::Vector{Float64})
end

objType = KTR_OBJTYPE_GENERAL #type of ojective function (full list available at: https://github.com/JuliaOpt/KNITRO.jl/blob/master/src/ktr_defines.jl)
objGoal = KTR_OBJGOAL_MINIMIZE #type of operation (full list available at: https://github.com/JuliaOpt/KNITRO.jl/blob/master/src/ktr_defines.jl)

n=2 #number of variables
#vectors of lower and upper bounds of variables
x_L = -KTR_INFBOUND*ones(2)
x_U = KTR_INFBOUND*ones(2)

m=2 #number constraints
c_Type = [KTR_CONTYPE_GENERAL,KTR_CONTYPE_GENERAL] #vector of constraints types (full list available at: https://github.com/JuliaOpt/KNITRO.jl/blob/master/src/ktr_defines.jl)
#vectors of constraints lower and upper bounds, here =0 because nonlinear equality constraints are equations of nonlinear system to solve
c_L = zeros(2)
c_U = zeros(2)

#jacIndexCons[i] and jacIndexVars[i] determine the row numbers and the column numbers, respectively, of the nonzero constraint Jacobian element jac[i]
#assume that the matrix is dense, so all entries should be filled with indices 
jac_con = vec(repmat(Int32[0,1],2,1))
jac_var = vcat(round(Int32,ones(2)*0),round(Int32,ones(2)*1))

#can leave hessian indices vectors empty
hess_row = Int32[]
hess_col = Int32[]

initg = [0.,0.] #initial guess
x = Array(Float64, length(initg)) #solution vector, updated when system solved 
x = initg[1:end]
lambda = zeros(n+m) #initial Lagrange multipliers on constraints (can be initialized to 0 even if constraints bind)
obj = [0] #initial value of objective function 

kp = createProblem()
loadOptionsFile(kp, "some_path/knitro.opt") # path of option file. change solver options in the latter
initializeProblem(kp, objGoal, objType, x_L, x_U, c_Type, c_L, c_U, jac_var, jac_con, hess_row, hess_col; initial_x=x, initial_lambda=lambda)
setCallbacks(kp, eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h, eval_hv)

println("initial guess: ", initg)
solveProblem(kp)
x = kp.x
exfl = kp.status 

initg_new = [1.,1.]

#if solver was unsuccessful with first initial guess, try another one
#full list of return codes (knitromatlab ref.) available at: https://www.artelys.com/tools/knitro_doc/3_referenceManual/knitromatlabReference.html#return-codes
if exfl < -199 
    println("new initial guess: ", initg_new)
    x = initg_new[1:end]
    restartProblem(kp, x, lambda)
    solveProblem(kp)
    x = kp.x
    exfl = kp.status
end

println("solution of nonlinear system: ", x)
println("system equations at solution: ", nleq(x)) `

And here is the option file:

#Artelys Knitro 10.1.1 Options file
#http://www.artelys.com/tools/knitro_doc/

#Which algorithm to use.
#auto   = 0 = let Knitro choose the algorithm
#direct = 1 = use Interior (barrier) Direct algorithm
#cg     = 2 = use Interior (barrier) CG algorithm
#active = 3 = use Active Set SLQP algorithm
#sqp    = 4 = use Active Set SQP algorithm
#multi  = 5 = run multiple algorithms (perhaps in parallel)
algorithm    1

#How to compute/approximate the gradient of the objective
#and constraint functions.
#exact        = 1 = user supplies exact first derivatives
#forward      = 2 = gradients computed by internal forward finite differences
#central      = 3 = gradients computed by internal central finite differences
#user_forward = 4 = gradients computed by user-provided forward finite differences
#user_central = 5 = gradients computed by user-provided central finite differences
gradopt      2

#How to compute/approximate the Hessian of the Lagrangian.
#exact           = 1 = user supplies exact second derivatives
#bfgs            = 2 = Knitro computes a dense quasi-Newton BFGS Hessian
#sr1             = 3 = Knitro computes a dense quasi-Newton SR1 Hessian
#product_findiff = 4 = Knitro computes Hessian-vector products by finite differences
#product         = 5 = user supplies exact Hessian-vector products
#lbfgs           = 6 = Knitro computes a limited-memory quasi-Newton BFGS Hessian
hessopt      2

#Whether to enforce satisfaction of simple bounds at all iterations.
#no      = 0 = allow iterations to violate the bounds
#always  = 1 = enforce bounds satisfaction of all iterates
#initpt  = 2 = enforce bounds satisfaction of initial point
honorbnds    1

#Maximum number of function evaluations to allow
#(a negative number implies no limit is imposed).
maxfevals    20000

#Maximum number of iterations to allow
#(if 0 then Knitro determines the best value).
#Default values are 10000 for NLP and 3000 for MIP.
maxit        100

#Specifies the final relative stopping tolerance for the KKT (optimality)
#error. Smaller values of opttol result in a higher degree of accuracy in
#the solution with respect to optimality.
opttol       1e-12

#Specifies the verbosity of output.
#none         = 0 = nothing
#summary      = 1 = only final summary information
#iter_10      = 2 = information every 10 iterations is printed
#iter         = 3 = information at each iteration is printed
#iter_verbose = 4 = more verbose information at each iteration is printed
#iter_x       = 5 = in addition, values of solution vector (x) are printed
#all          = 6 = in addition, constraints (c) and multipliers (lambda)
outlev       iter_verbose

#Whether to allow simultaneous evaluations in parallel.
#no   = 0 = only one thread can perform an evaluation at a time
#yes  = 1 = allow multi-threaded simultaneous evaluations
par_concurrent_evals  no

#Strategy for setting initial x, lambda and slacks with barrier algorithms.
#This option only affects the initial x value when not provided by user.
#auto   = 0 = let Knitro choose the strategy
#strat1 = 1 = initial point strategy 1 (mainly for LP/QP)
#strat2 = 2 = initial point strategy 2
#strat3 = 3 = initial point strategy 3
bar_initpt   strat3

#Which barrier parameter update strategy.
#auto     = 0 = let Knitro choose the strategy
#monotone = 1
#adaptive = 2
#probing  = 3
#dampmpc  = 4
#fullmpc  = 5
#quality  = 6
bar_murule   dampmpc

#Whether or not to penalize constraints in the barrier algorithms.
#auto       = 0 = let Knitro choose the strategy
#none       = 1 = Do not apply penalty approach to any constraints
#all        = 2 = Apply a penalty approach to all general constraints
#equalities = 3 = Apply a penalty approach to equality constraints only
bar_penaltycons   none

#Which penalty parameter update strategy for barrier algorithms.
#auto     = 0 = let Knitro choose the strategy
#single   = 1 = use single penalty parameter approach
#flex     = 2 = use more tolerant flexible strategy
bar_penaltyrule   flex

#Switching rule strategy for barrier algorithms that controls
#switching between optimality and feasibility phases.
#auto     = 0 = let Knitro choose the strategy
#never    = 1 = never switch
#level1   = 2 = allow moderate switching
#level2   = 3 = more agressive switching
bar_switchrule    never

#Which linear system solver to use.
#auto       = 0 = let Knitro choose the solver
#internal   = 1 = use internal solver provided with Knitro
#(not currently active; reserved for future use)
#hybrid     = 2 = use a mixture of linear solvers depending on the linear systems
#qr         = 3 = use dense QR solver always (only for small problems)
#ma27       = 4 = use sparse HSL solver ma27 always
#ma57       = 5 = use sparse HSL solver ma57 always
#mklpardiso = 6 = use sparse Intel MKL Pardiso solver always
linsolver    hybrid

#Whether to use the Knitro Tuner.
#off    = 0 = Knitro Tuner turned off
#on     = 1 = Knitro Tuner enabled
tuner        off

#Specifies the final relative stopping tolerance for the feasibility
#error. Smaller values of feastol result in a higher degree of accuracy
#in the solution with respect to feasibility.
feastol      1e-12

#Specifies the final absolute stopping tolerance for the feasibility error.
#Smaller values of feastol_abs result in a higher degree of accuracy in the
#solution with respect to feasibility.
#feastol_abs  0.001

#Step size tolerance used for terminating the optimization.
xtol         1e-15

#Whether to perform scaling of the problem.
#no            = 0 = no scaling done
#user_internal = 1 = user, if defined, otherwise internal
#user_none     = 2 = user, if defined, otherwise none
#internal      = 3 = Knitro performs internal scaling
scale        no

#Initial trust region radius scaling factor, used to determine
#the initial trust region size.
delta        1

#When using the Interior/Direct algorithm, this parameter
#controls the maximum number of consecutive CG steps before
#trying to force the algorithm to take a direct step again.
#(only used for alg=1).
bar_directinterval  0

#Whether feasibility is given special emphasis.
#no       = 0 = no emphasis on feasibility
#stay     = 1 = iterates must honor inequalities
#get      = 2 = emphasize first getting feasible before optimizing
#get_stay = 3 = implement both options 1 and 2 above
bar_feasible stay

[yeesian]

For clarity: is the reason why it was stuck in your original due to the following (since gvec doesn't change)?

function eval_g(x::Vector{Float64}, cons::Vector{Float64}) # constraints function
    temp = solveAtPoint(mobj, x, gvec, point, 0, trans)[1]; # produces vector residuals we want to set to 0
    cons[1] = temp[1]
    # ...
    cons[14] = temp[14]
    # return cons
end

[pierremabille]

No, it was due to the specifications of derivatives sparsity patterns, and options for computing derivatives in the knitro.opt files (strangely when you print a detailed output in matlab, it sets derivative computation options to 4, a value that doesn't exist for gradopt and hessopt -- here it works for gradopt 2, forward differences, and 3, central differences, and various values for hessopt). Basically knitro.jl was acting as if the Jacobian was 0 everywhere, which led it to be stuck at the initial guess.

[mlubin]

(@pierremabille, just so you know, for any function that you can compute on a computer you can in principle compute its exact derivatives by using techniques from automatic differentiation. Having analytic expressions is irrelevant, in principle.)