// example: read quadratic program in MPS format from file

// the QP below is the first quadratic program example in the user manual

#include <iostream>

#include <fstream>

#include <CGAL/basic.h>

#include <CGAL/QP_models.h>

#include <CGAL/QP_functions.h>

 

// choose exact integral type

#ifdef CGAL_USE_GMP

#include <CGAL/Gmpz.h>

typedef CGAL::Gmpz ET;

#else

#include <CGAL/MP_Float.h>

typedef CGAL::MP_Float ET;

#endif

 

// program and solution types

typedef CGAL::Quadratic_program_from_mps<int> Program;

typedef CGAL::Quadratic_program_solution<ET> Solution;

 

int main() {

std::ifstream in ("first_qp.mps");

Program qp(in); // read program from file

assert (qp.is_valid()); // we should have a valid mps file

 

// solve the program, using ET as the exact type

Solution s = CGAL::solve_quadratic_program(qp, ET());

 

// output solution

std::cout << s;

return 0;

}

 

3.3 Constructing a Program from Iterators

The following program again solves our running example from above, with the same output, but this time with iterators over data stored in suitable containers. You can see that we also store zero entries here (in D). For this toy problem, the previous two approaches (program from data/stream) are clearly preferable, but Section Working from Iterators shows an example where it makes sense to use the iterator-based approach.

下面的程序与上面解决的问题相同，但通过合适容器的相关迭代器。你可以看到我们也将0项保存在D矩阵中。对于本例，前面的两个方式 （输入数据为data和stream）更加好，但 Working from Iterators节显示了迭代器方法在什么地方更加合理。

 

让我们考虑上面Which Programs can be Solved?节提出的（QP）问题的通用形式。 如果D=0，则二次规划是一个线性规划。当边界向量I是零向量并且u的所有项为 ∞，这样的规划被称为非负的（nonnegative）。本包为此特殊类型提供了专门的模型和方法。

从接口视角，这只是在说法上比较容易（syntactic sugar:）：在模型中，我们很容易通过设置缺省边界值得到一个非负规划（nonnegative program），一个线性规划则通过不插入任何D项得到。即使在基于迭代器的方法中 (见 QP_solver/first_qp_from_iterators.cpp)，线性和非负规划可以轻易地通过合适的Const_oneset_iterator<T>样式迭代器得到。

用专门的方法解决线性和非负规划问题的原因在于效率：如果解析器知道这一规划是线性的，则可以节约一些计算资源，因为不需要对不必要的矩阵进行纠缠。正如上面Robustness 节提到的，可以提前检查D是否为0是不恰当的，因为当稠密接口时这需要Ω(n²)时间。

同样，当解析器知道当前的规划是非负的，则它的效率将比l≤x≤u存在上下界时提高很多。所以我们使用概念NonnegativeLinearProgram 而不使用 LinearProgram。

很多情况下，变量根本就没有界，即 l 所有的项为−∞，u所有的项为∞，这种情况即为*规划（free program）。没有专门的方法用于这种情况（一个*二次规划或线性二次规划被当作一个通用的二次规划或线性规划），但所有的所有的预定义的模型使指定所有类型的缺省界变得容易，也包括了*规划的情况。

 

4.1 The Linear Programming Solver

Let's go back to our first quadratic program from above and change it into a linear program by simply removing the quadratic part of the objective function:

线性规划由二次规划转化过来，只需去掉目标函数中的二次项：

Figure 7.2 shows how this looks like. We will not visualize a linear objective function with contour lines but with arrows instead. The arrow represents the (direction) of the vector −c, and we are looking for a feasible solution that is "extreme" in the direction of the arrow. In our small example, this is the unique point "on" the two constraints x1+x2≤7 and −x1+x2≤4, the point (10/3,11/3) marked with a black dot. The optimal objective function value is −32(11/3)+64=−160/3.

图中给出了线性规划的形态。我们采用一个箭标表示线性目标函数。箭头表示向量−c的（方向），我们正在寻找一个可行的解，这个解在箭头方向上是一个“极值”。在我们的例子中，这个“极值”为图上黑点标识的点，它处于两个约束 x1+x2≤7和−x1+x2≤4的交点 (10/3,11/3) 上。其最优目标函数值是−32(11/3)+64=−160/3。

Here is CGAL code for solving it, using the dedicated LP solver, and according to the three ways for constructing a program that we have already discussed in Section How to Enter and Solve a Program.

QP_solver/first_lp.cpp

QP_solver/first_lp_from_mps.cpp

QP_solver/first_lp_from_iterators.cpp

In all cases, the output is

status: OPTIMAL

objective value: -160/3

variable values:

0: 10/3

1: 11/3

这里是CGAL解的代码，使用了LP解析器，使用了How to Enter and Solve a Program节的三种方法

4.2 The Nonnegative Quadratic Programming Solver

If we go back to our first quadratic program and remove the constraint y≤4, we arrive at a nonnegative quadratic program:

在上面的例子中去除 y≤4,就得到一个非负二次规划。

Figure 7.3 contains the illustration; since the constraint y≤4 was redundant, the feasible region and the optimal solution do not change.

图中显示：因为 y≤4,是冗余的，所以可行区域和最优解没有变化

The following programs (using the dedicated solver for nonnegative quadratic programs) will therefore again output

status: OPTIMAL

objective value: 8/1

variable values:

0: 2/1

1: 3/1

QP_solver/first_nonnegative_qp.cpp

QP_solver/first_nonnegative_qp_from_mps.cpp

QP_solver/first_nonnegative_qp_from_iterators.cpp

4.3 The Nonnegative Linear Programming Solver

Finally, a dedicated model and function is available for nonnnegative linear programs as well. Let's take our linear program from above and remove the constraint y≤4 to obtain a nonnegative linear program. At the same time we remove the constant objective function term to get a "minimal" input and a "shortest" program; the optimal value is −32(11/3)=−352/3.

最后，是专门为非负线性规划提供的模型和函数。本例始于上面的例子并去掉了约束y≤4得到了一个非负线性规划。同时，我们去除目录函数的常数项得到一个“最小”输入和一个“最短”规划，其最优值是−32(11/3)=−352/3

This can be solved by any of the following three programs

QP_solver/first_nonnegative_lp.cpp

QP_solver/first_nonnegative_lp_from_mps.cpp

QP_solver/first_nonnegative_lp_from_iterators.cpp

The output will always be

status: OPTIMAL

objective value: -352/3

variable values:

0: 10/3

1: 11/3

5 Working from Iterators

Here we present a somewhat more advanced example that emphasizes the usefulness of solving linear and quadratic programs from iterators. Let's look at a situation in which a linear program is given implicitly, and access to it is gained through properly constructed iterators.

The problem we are going to solve is the following: given points p₁,…p_n in d-dimensional space and another point p: is p in the convex hull of {p1,…,pn}? In formulas, this is the case if and only if there are real coefficients λ1,…,λn such that p is a convex combination of p₁,…,p_n:

The problem of testing the existence of such λ_j can be expressed as a linear program. It becomes particularly easy when we use the homogeneous representations of the points: if q₁,…,q_n,q∈R^d+1 are homogeneous coordinates for p₁,…,p_n,p with positive homogenizing coordinates h₁,…,h_n,h, we have

Now, nonnegative λ₁,…,λ_n are suitable coefficients for a convex combination if and only if

equivalently, if there are μ₁,…,μ_n (with μ_j=λ_j⋅h/h_j for all j) such that

The linear program now tests for the existence of nonnegative μj that satisfy the latter equation. Below is the code; it defines a function that solves the linear program, given p and p1,…,pn (through an iterator range). The only (mild) trickery involved is the construction of the nested iterator through a fixed column of the constraint matrix A. We get this from transforming the iterator through the points using a functor that maps a point to an iterator through its homogeneous coordinates.

这里我们给出一个更加高级的例子，它强调在解决线性或二次规划时迭代器的作用。当给出兵线性规划是隐含的，并且通过创建合适的迭代器来接入问题。

我们的问题如下：给定d维空间的点集p₁,…p_n和另一个点p，确定p点是否在{p₁,…p_n}点集形成的凸包（convex hull）中？ 由定理得知，当且仅当存在一个实数系统集λ1,…,λn 使得p是p₁,…p_n的凸组合

时本问题成立。

而测试是否存在的问题可以表达为一个线性规划问题。如果q₁,…,q_n,q∈R^d+1是具有正齐次化坐标h₁,…,h_n,h的点p₁,…,p_n,p的齐次坐标，我们使用齐次坐标表示点就十分方便（It becomes particularly easy when we use the homogeneous representations of the points: if q₁,…,q_n,q∈R^d+1 are homogeneous coordinates for p₁,…,p_n,p with positive homogenizing coordinates h₁,…,h_n,h）。我们有：

如果存在一组系数 λ₁,…,λ_n ，当且仅当

时，则 λ₁,…,λ_n 是适合的非负的凸组合系数。

同样地，如果有一组数 μ₁,…,μ_n (对于 所有的j有 μ_j=λ_j⋅h/h_j ) 满足

现在，线性规划测试满足最后等式的非负的 μ_j 是否存在，下面是其代码；它定义了一个给定p和p₁,…p_n(通过iterator给出)解线性规划的函数。这里有点技巧的是通过一个固定列的约束矩阵A创建内嵌的（nested iterator）迭代器。我们通过（We get this from transforming the iterator through the points using a functor that maps a point to an iterator through its homogeneous coordinates）

File QP_solver/solve_convex_hull_containment_lp.h

To see this in action, let us call it with p1=(0,0),p2=(10,0),p3=(0,10) fixed (they define a triangle) and all integral points p in [0,10]2. We know that p is in the convex hull of {p1,p2,p3} if and only if its two coordinates sum up to 10 at most. As the exact type, we use MP_Float or Gmpzf (which is faster and preferable if GMP is installed).

File QP_solver/convex_hull_containment.cpp

5.1 Using Makers

You already noticed in the previous example that the actual template arguments for Nonnegative_linear_program_from_iterators<A_it, B_it, R_it, C_it> can be quite elaborate, and this only gets worse if you plug more iterators into each other. In general, you want to construct a program from given expressions for the iterators, but the types of these expressions are probably very complicated and difficult to look up.

You can avoid the explicit construction of the type Nonnegative_linear_program_from_iterators<A_it, B_it, R_it, C_it> if you only need an expression of it, e.g. to pass it directly as an argument to the solving function. Here is an alternative version of QP_solver/solve_convex_hull_containment_lp.h that shows how this works. In effect, you get shorter and more readable code.

前一个例子中，你已经注意到Nonnegative_linear_program_from_iterators<A_it, B_it, R_it, C_it>的实际模板参数相当复杂。如果你在它们之间插入更多的iterator后。通常，我们想通过一个给定的关于iterator的表达式创建一个规划，但是表达式的类型可能会非常复杂并难以查找。

如果你只需要一个它的表达式，可以通过将它作为参数传递到解析函数就可以避免显式地创建Nonnegative_linear_program_from_iterators<A_it, B_it, R_it, C_it>类型。这里有一个QP_solver/solve_convex_hull_containment_lp.h的替代版本就是展示其方法的。通过这个方法可以得到更加可读的代码。
File QP_solver/solve_convex_hull_containment_lp2.h

6 Important Variables and Constraints

If you have a solution x∗ of a linear or quadratic program, the "important" variables are typically the ones that are not on their bounds. In case of a nonnegative program, these are the nonzero variables. Going back to the example of the previous Section Working from Iterators, we can easily interpret their importance: the nonzero variables correspond to points pj that actually contribute to the convex combination that yields p.

The following example shows how we can access the important variables, using the iterators basic_variable_indices_begin() and basic_variable_indices_end().

We generate a set of points that form a 4-gon in [0,4]2, and then find the ones that contribute to the convex combinations of all 25 lattice points in [0,4]2. If the lattice point in question is not in the 4-gon, we simply output this fact.

如果你有线性或二次规划的一个解x^∗，“重要的”变量是哪些典型地在边界之外的变量。在非负规划的情况，有非零变量。在Working from Iterators节的例子中，我们可以容易解释它们的重要性：非零变量与生成关于p的凸组合做出实际贡献的p_j对应。（Going back to the example of the previous Section Working from Iterators, we can easily interpret their importance: the nonzero variables correspond to points pj that actually contribute to the convex combination that yields p）

File QP_solver/important_variables.cpp

It turns out that exactly three of the four points contribute to any convex combination, even through there are lattice points that lie in the convex hull of less than three of the points. This shows that the set of basic variables that we access in the example does not necessarily coincide with the set of important variables as defined above. In fact, it is only guaranteed that a non-basic variable attains one of its bounds, but there might be basic variables that also have this property. In linear and quadratic programming terms, such a situation is called a degeneracy.

There is also the concept of an important constraint: this is typically a constraint in the system Ax⋛b that is satisfied with equality at x∗. Program QP_solver/first_qp_basic_constraints.cpp shows how these can be accessed, using the iterators basic_constraint_indices_begin() and basic_constraint_indices_end().

Again, we have a disagreement between "basic" and "important": it is guaranteed that all basic constraints are satisfied with equality at x∗, but there might be non-basic constraints that are satisfied with equality as well.

结果是四个点中的三个对任意 凸组合有贡献，即使存在格点处于少于3个点的凸包内（even through there are lattice points that lie in the convex hull of less than three of the points）。

7 Solution Certificates

Suppose the solver tells you that the problem you have entered is infeasible. Why should you believe this? Similarly, you can quite easily verify that a claimed optimal solution is feasible, but why is there no better one?

Certificates are proofs that the solver can give you in order to convince you that what it claims is indeed true. The archetype of such a proof is Farkas Lemma [1].

Farkas Lemma: Either the inequality system

has a solution x∗, or there exists a vector y such that

but not both.

Thus, if someone wants to convince you that the first system in the Farkas Lemma is infeasible, that person can simply give you a vector y that solves the second system. Since you can easily verify yourself that the y you got satisfies this second system, you now have a certificate for the infeasibility of the first system, assuming that you believe in Farkas Lemma.

Here we show how the solver can convince you. We first set up an infeasible linear program with constraints of the type Ax≤b,x≥0; then we solve it and ask for a certificate. Finally, we verify the certificate by simply checking the inequalities of the second system in Farkas Lemma.

File QP_solver/infeasibility_certificate.cpp

There are similar certificates for optimality and unboundedness that you can see in action in the programs QP_solver/optimality_certificate.cpp and QP_solver/unboundedness_certificate.cpp. The underlying variants of Farkas Lemma are somewhat more complicated, due to the mixed relations in ⋛ and the general bounds. The certificate section of Quadratic_program_solution<ET> gives the full picture and mathematically proves the correctness of the certificates.

8 Customizing the Solver

Sometimes it is necessary to alter the default behavior of the solver. This can be done by passing a suitably prepared object of the class Quadratic_program_options to the solution functions. Most options concern "soft" issues like verbosity, but there are two notable case where it is of critical importance to be able to change the defaults.

8.1 Exponent Overflow in Double Using Floating-Point Filters

The filtered version of the solver that is used for some problems by default on input type double internally constructs double-approximations of exact multiprecision values. If these exact values are extremely large, this may lead to infinite double values and incorrect results. In debug mode, the solver will notice this through a certificate cross-check in the end (or even earlier). In this case, it is advisable to explicitly switch to a non-filtered pricing strategy, see Quadratic_program_pricing_strategy.

Hint: If you have a program where the number of variables n and the number of constraints m have the same order of magnitude, the filtering will usually have no dramatic effect on the performance, so in that case you might as well switch to QP_PARTIAL_DANTZIG to be safe from the issue described here (see QP_solver/cycling.cpp for an example that shows how to change the pricing strategy).

8.2 The Solver Internally Cycles

Consider the following program. It reads a nonnegative linear program from the file cycling.mps (which is in the example directory as well), and then solves it in verbose mode, using Bland's rule, see Quadratic_program_pricing_strategy. 
File QP_solver/cycling.cpp

If you comment the line

options.set_pricing_strategy(CGAL::QP_BLAND); // Bland's rule

you will see that the solver cycles: the verbose mode outputs the same sequence of six iterations over and over again. By switching to QP_BLAND, the solution process typically slows down a bit (it may also speed up in some cases), but now it is guaranteed that no cycling occurs.

In general, the verbose mode can be of use when you are not sure whether the solver "has died", or whether it simply takes very long to solve your problem. We refer to the class Quadratic_program_options for further details.

9 Some Benchmarks for Convex Hull Containment

Here we want to show what you can expect from the solver's performance in a specific application; we don't know whether this application is typical in your case, and we make no claims whatsoever about the performance in other applications.

Still, the example shows that the performance can be dramatically affected by switching between pricing strategies, and we give some hints on how to achieve good performance in general.

The application is the one already discussed in Section Working from Iterators above: testing whether a point is in the convex hull of other points. To be able to switch between pricing strategies, we add another parameter of type Quadratic_program_options to the function solve_convex_hull_containment_lp that we pass on to the solution function:

File QP_solver/solve_convex_hull_containment_lp3.h

Now let us test containment of the origin in the convex hull of n random points in [0,1]d (it will most likely not be contained, and it turns out that this is the most expensive case). In the program below, we use d=10 and n=100,000, and we comment on some other combinations of n and d below (feel free to experiment with still other values).

File QP_solver/convex_hull_containment_benchmarks.cpp

If you compile with the macros NDEBUG or CGAL_QP_NO_ASSERTIONS set (this is essential for good performance!!), you will see runtimes that qualitatively look as follows (on your machine, the actual runtimes will roughly be some fixed multiples of the numbers in the table below, and they might vary with the random choices). The default choice of the pricing strategy in that case is QP_PARTIAL_FILTERED_DANTZIG.

Strategy	Runtime in seconds
QP_CHOOSE_DEFAULT	0.32
QP_DANTZIG	10.7
QP_PARTIAL_DANTZIG	3.72
QP_BLAND	3.65
QP_FILTERED_DANTZIG	0.43
QP_PARTIAL_FILTERED_DANTZIG	0.32

We clearly see the effect of filtering: we gain a factor of ten, roughly, compared to the next best non-filtered variant.

9.1 d=3, n=1,000,000

The filtering effect is amplified if the points/dimension ratio becomes larger. This is what you might see in dimension three, with one million points.

Strategy	Runtime in seconds
QP_CHOOSE_DEFAULT	1.34
QP_DANTZIG	47.6
QP_PARTIAL_DANTZIG	15.6
QP_BLAND	16.02
QP_FILTERED_DANTZIG	1.89
QP_PARTIAL_FILTERED_DANTZIG	1.34

In general, if your problem has a high variable/constraint or constraint/variable ratio, then filtering will typically pay off. In such cases, it might be beneficial to encode your problem using input type double in order to profit from the filtering (but see the issue discussed in Section Exponent Overflow in Double Using Floating-Point Filters).

9.2 d=100, n=100,000

Conversely, the filtering effect deteriorates if the points/dimension ratio becomes smaller.

Strategy	Runtime in seconds
QP_CHOOSE_DEFAULT	3.05
QP_DANTZIG	78.4
QP_PARTIAL_DANTZIG	45.9
QP_BLAND	33.2
QP_FILTERED_DANTZIG	3.36
QP_PARTIAL_FILTERED_DANTZIG	3.06

9.3 d=500, n=1,000

If the points/dimension ratio tends to a constant, filtering is no longer a clear winner. The reason is that in this case, the necessary exact calculations with multiprecision numbers dominate the overall runtime.

Strategy	Runtime in seconds
QP_CHOOSE_DEFAULT	2.65
QP_DANTZIG	5.55
QP_PARTIAL_DANTZIG	5.6
QP_BLAND	4.46
QP_FILTERED_DANTZIG	2.65
QP_PARTIAL_FILTERED_DANTZIG	2.61

In general, if you have a program where the number of variables and the number of constraints have the same order of magnitude, then the saving gained from using the filtered approach is typically small. In such a situation, you should consider switching to a non-filtered variant in order to avoid the rare issue discussed in Section Exponent Overflow in Double Using Floating-Point Filters altogether.