NMath/SolverFoundation 声称非凸模型，尽管可以通过 Matlab "quadprog" 求解
NMath/SolverFoundation claims non-convex model, despite being solvable by Matlab's "quadprog"
我有一个由 Hessian H 和向量 c 定义的二次规划问题。这个问题在使用 Matlab 的 quadprog 函数解决的参考文献中
result = quadprog(H, c, [], [], [], [], lb, ub, init, opts);
lb 和 ub 是与 c 长度相同的下限和上限向量，并且 init 无论如何都会被忽略，因为 opts 定义算法为 interior-point-convex.
我的理解是每个参数都有下限和上限（0 和 Inf），但它们之间没有进一步的线性约束。
我现在想在 .NET 中解决这个问题。我首先尝试在 MS Solver Foundation 中实现它，但他们对所有内容的特定符号让我很头疼，所以经过更多搜索后我发现 NMath 有一个实现，所以我试了一下。
我先定义了问题
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
这里的 n 是 H 中的行数，也就是要求解的变量数。我为每个变量添加了界限。
H 和 c 包含与 Matlab 中相同的值。
我现在定义了求解器和求解器选项（基于 NMath 教程中的一些值）
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
并使用
调用求解器
Solver.Solve(Problem, SolverParams)
这会引发异常
Microsoft.SolverFoundation.Common.InvalidModelDataException: 'The model is not convex.'
有内部 Division by Zero 异常。
如果这是唯一的事情，我会假设我在 H 或 c 中有问题，但事实上 Matlab 可以毫无问题地解决这个确切的模型让我有点不知所措。
这两个模型有什么不同？事实上，在 Matlab 中，参数 A、b、Aeq 和 beq 被设置为 []，这意味着不等式和等式约束都不存在，因此不为 NMath 模型设置任何参数是有意义的。
有人指出我哪里出错了吗？
示例值（NMath 表示法）
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给出解值：
0
99.9802
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.1653
是的，事实证明提出的问题确实是非凸的，因此给定的求解器无法求解。
用于非凸、框约束优化的特殊算法更加复杂且计算困难（NP-hard）。
我可以使用 alglib 来解决这些问题。他们甚至有一个示例，我只需要进行最少的更改即可包含额外的线性组件：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
我有一个由 Hessian H 和向量 c 定义的二次规划问题。这个问题在使用 Matlab 的 quadprog 函数解决的参考文献中
result = quadprog(H, c, [], [], [], [], lb, ub, init, opts);
lb 和 ub 是与 c 长度相同的下限和上限向量，并且 init 无论如何都会被忽略，因为 opts 定义算法为 interior-point-convex.
我的理解是每个参数都有下限和上限（0 和 Inf），但它们之间没有进一步的线性约束。
我现在想在 .NET 中解决这个问题。我首先尝试在 MS Solver Foundation 中实现它，但他们对所有内容的特定符号让我很头疼，所以经过更多搜索后我发现 NMath 有一个实现，所以我试了一下。
我先定义了问题
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
这里的 n 是 H 中的行数，也就是要求解的变量数。我为每个变量添加了界限。
H 和 c 包含与 Matlab 中相同的值。
我现在定义了求解器和求解器选项（基于 NMath 教程中的一些值）
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
并使用
调用求解器Solver.Solve(Problem, SolverParams)
这会引发异常
Microsoft.SolverFoundation.Common.InvalidModelDataException: 'The model is not convex.'
有内部 Division by Zero 异常。
如果这是唯一的事情，我会假设我在 H 或 c 中有问题，但事实上 Matlab 可以毫无问题地解决这个确切的模型让我有点不知所措。
这两个模型有什么不同？事实上，在 Matlab 中，参数 A、b、Aeq 和 beq 被设置为 []，这意味着不等式和等式约束都不存在，因此不为 NMath 模型设置任何参数是有意义的。
有人指出我哪里出错了吗？
示例值（NMath 表示法）
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给出解值：
0
99.9802
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.1653
是的，事实证明提出的问题确实是非凸的，因此给定的求解器无法求解。
用于非凸、框约束优化的特殊算法更加复杂且计算困难（NP-hard）。
我可以使用 alglib 来解决这些问题。他们甚至有一个示例，我只需要进行最少的更改即可包含额外的线性组件：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