NMath/SolverFoundation claims non-convex model, despite being solvable by Matlab's "quadprog"

I have a quadratic programming problem defined by the Hessian H and the vector c. This problem is in a reference solved using Matlab's quadprog function like that

result = quadprog(H, c, [], [], [], [], lb, ub, init, opts);

lb and ub are vectors of lower and upper bounds of the same length as c, and init is ignored anyway since opts defines the algorithm to be interior-point-convex.

I understand this in terms of a problem that has lower and upperbounds (0 and Inf) for each of the parameters, but no further linear constraints between them.

I now want to solve this problem in .NET. I first tried to implement it in MS Solver Foundation but their specific notation for everything gave me a headache, so after some more searching I found that NMath has an implementation, so I gave their trial a go.

I first defined the problem

Dim Problem As New CenterSpace.NMath.Analysis.QuadraticProgrammingProblem(H, c)
For i = 1 To n
    Problem.AddLowerBound(i - 1, 0)
    Problem.AddUpperBound(i - 1, Double.PositiveInfinity)
Next

Here n is the number of rows in the H, aka the number of variables to solve for. For each variable I add the bounds. H and c contain the same values as in Matlab.

I now defined the solver and solver options (based on some values in an NMath tutorial)

Dim Solver As New InteriorPointQPSolver()

Dim SolverParams As New InteriorPointQPSolverParams
SolverParams.KktForm = InteriorPointQPSolverParams.KktFormOption.Blended
SolverParams.Tolerance = "1e-6"
SolverParams.MaxDenseColumnRatio = 0.9
SolverParams.PresolveLevel = InteriorPointQPSolverParams.PresolveLevelOption.Automatic
SolverParams.SymbolicOrdering = InteriorPointQPSolverParams.SymbolicOrderingOption.Automatic

and call the solver using

Solver.Solve(Problem, SolverParams)

This throws an exception

Microsoft.SolverFoundation.Common.InvalidModelDataException: 'The model is not convex.'

with inner Division by Zero exception.

If this was the only thing to go by, I would just assume I have something wrong in H or c, but the fact that Matlab can solve this exact model without problem leaves me a bit at a loss.

Where do the two models differ? The fact that in Matlab the arguments A, b, Aeq and beq are set to [] means that neither inequality nore equality constraints exist, so it makes sense not to set any for the NMath model.

Does someone have a pointer where I am going wrong?

Example values (NMath notation)

H
"11x11 [ 0 0 0 0 0 0 0 0 0 0 0  0 9 0.0362597318967715 0.0589000506405063 0.0949592651496192 0.151356936171003 0.237273856132009 0.363526400836095 0.540694082222912 0.776143076440992 1.07111504928792  0 0.0362597318967715 0.0100090487691427 0.0216132625192289 0.0333985757439813 0.0392023626107542 0.0464624233197008 0.0560189763874193 0.0662001400288286 0.0757654577788433 0.0835891629836335  0 0.0589000506405063 0.0216132625192289 0.0215096597923902 0.0361130413465718 0.0514313528568181 0.0611682501074947 0.0724431110814198 0.0855864103736708 0.0983346261825746 0.108993492590539  0 0.0949592651496192 0.0333985757439813 0.0361130413465718 0.0398513724131197 0.0590171495132719 0.0794702824167896 0.0945317161015367 0.110685477350358 0.127501862369773 0.142107023175787  0 0.151356936171003 0.0392023626107542 0.0514313528568181 0.0590171495132719 0.0685957299907155 0.0944181301705371 0.121897035464389 0.143574330419051 0.164966499696443 0.18482585143016  0 0.237273856132009 0.0464624233197008 0.0611682501074947 0.0794702824167896 0.0944181301705371 0.112508953001228 0.147491705343541 0.183846332946191 0.212902676348244 0.239092382667909  0 0.363526400836095 0.0560189763874193 0.0724431110814198 0.0945317161015367 0.121897035464389 0.147491705343541 0.177279888808415 0.223930012345989 0.270439636810068 0.306727776312255  0 0.540694082222912 0.0662001400288286 0.0855864103736708 0.110685477350358 0.143574330419051 0.183846332946191 0.223930012345989 0.26858976843593 0.328738895388255 0.385593364732688  0 0.776143076440992 0.0757654577788433 0.0983346261825746 0.127501862369773 0.164966499696443 0.212902676348244 0.270439636810068 0.328738895388255 0.390566939588293 0.464670732441943  0 1.07111504928792 0.0835891629836335 0.108993492590539 0.142107023175787 0.18482585143016 0.239092382667909 0.306727776312255 0.385593364732688 0.464670732441943 0.544201571301849 ]"

c
"[ 0 -899.999332235966 -3.63818566995197 -5.9055875689195 -9.5157467986296 -15.160783363483 -23.7589207356301 -36.3918800495691 -54.1175844999643 -77.6724662955247 -107.180324535075 ]"

Matlab gives the solution values:

Solution

Right, it turns out that the problem posed is indeed non-convex and thus not solvable by the given solver.

Special algorithms for non-convex, box-constrained optimization are more complicated and computationally hard (NP-hard).

I was able to solve these using alglib instead. They even had an example I only had to minimally change to include the additional linear components: http://www.alglib.net/translator/man/manual.vbnet.html#example_minqp_d_nonconvex

'Solve the quadratic programming program using alglib
Dim init(H.Rows - 1) As Double
Dim lb(H.Rows - 1) As Double
Dim ub(H.Rows - 1) As Double
For i = 0 To init.Count - 1
    init(i) = 1
    lb(i) = 0
    ub(i) = Double.PositiveInfinity
Next

Dim a(,) As Double = H.ToDoubleMatrix2

Dim x() As Double = New Double() {}
Dim state As alglib.minqpstate = Nothing
Dim rep As alglib.minqpreport = Nothing

'  create solver, set quadratic/linear terms, constraints
alglib.minqpcreate(init.Count, state)
alglib.minqpsetquadraticterm(state, a)
alglib.minqpsetlinearterm(state, c.Row(1))
alglib.minqpsetstartingpoint(state, init)
alglib.minqpsetbc(state, lb, ub)

'  Set scale of the parameters.
'  It is strongly recommended that you set scale of your variables.
'  Knowing their scales is essential for evaluation of stopping criteria
'  and for preconditioning of the algorithm steps.
'  You can find more information on scaling at http://www.alglib.net/optimization/scaling.php
alglib.minqpsetscale(state, init)

'  solve problem with BLEIC-QP solver.
'  default stopping criteria are used.
alglib.minqpsetalgobleic(state, 0.0, 0.0, 0.0, 0)
alglib.minqpoptimize(state)
alglib.minqpresults(state, x, rep)
Return x