PIQP solves quadratic programs of the form
\[ \begin{aligned} \min_{x} \quad & \frac{1}{2} x^\top P x + c^\top x \\ \text {s.t.}\quad & Ax=b, \\ & Gx \leq h, \\ & x_{lb} \leq x \leq x_{ub} \end{aligned} \]
with primal decision variables \(x \in \mathbb{R}^n\), matrices \(P\in \mathbb{S}_+^n\), \(A \in \mathbb{R}^{p \times n}\), \(G \in \mathbb{R}^{m \times n}\), and vectors \(c \in \mathbb{R}^n\), \(b \in \mathbb{R}^p\), \(h \in \mathbb{R}^m\), \(x_{lb} \in \mathbb{R}^n\), and \(x_{ub} \in \mathbb{R}^n\).
Consider:
\[ \begin{aligned} \min_{x} \quad & \frac{1}{2} x^\top \begin{bmatrix} 6 & 0 \\ 0 & 4 \end{bmatrix} x + \begin{bmatrix} -1 \\ -4 \end{bmatrix}^\top x \\ \text {s.t.}\quad & \begin{bmatrix} 1 & -2 \end{bmatrix} x = 1, \\ & \begin{bmatrix} 1 & -1 \\ 2 & 0 \end{bmatrix} x \leq \begin{bmatrix} 0.2 \\ -1 \end{bmatrix}, \\ & -1 \leq x_1 \leq 1. \end{aligned} \]
The data for this problem can be specified as below.
P <- matrix(c(6, 0, 0, 4), nrow = 2)
c <- c(-1, -4)
A <- matrix(c(1, -2), nrow = 1)
b <- 1
G <- matrix(c(1, 2, -1, 0), nrow = 2)
h <- c(0.2, -1)
x_lb <- c(-1, -Inf) ## 2 variables
x_ub <- c(1, Inf) ## 2 variablesThe problem can now be solved via a call to
solve_piqp().
sol <- solve_piqp(P, c, A, b, G, h, x_lb = x_lb, x_ub = x_ub, backend = "auto")
cat(sprintf("(Solution status, description): = (%d, %s)\n",
sol$status, sol$info$status_desc))
#> (Solution status, description): = (1, solved)
cat(sprintf("Objective: %f, solution: (x1, x2) = (%f, %f)\n", sol$info$primal_obj, sol$x[1], sol$x[2]))
#> Objective: 6.160000, solution: (x1, x2) = (-0.600000, -0.800000)sol contains many components as str(sol)
will display but the most important ones are:
status : 1 if all goes well (more below),x : solution vectory : dual solution for the equality constraintsz : dual solution for the inequality constraintsz_lb : dual solution of lower bound box
constraintsz_ub : dual solution of upper bound box
constraintsinfo$status_desc: a descriptive string of the
statusinfo$primal_pobj : primal objective valueinfo$run_time : total runtime, if asked for in settings
(see below).One can always construct the descriptive string for the status using:
Note that PIQP can handle infinite box constraints well, i.e. when
elements of \(x_{lb}\) or \(x_{ub}\) are \(-\infty\) or \(\infty\), respectively. On the contrary,
infinite values in the general inequalities \(Gx \leq h = \pm \infty\) can cause
problems, i.e., they are converted internally to -1e30 and
1e30, respectively.
Users who wish to solve QP problems will mostly use the
solve_piqp() function. Behind the scenes,
solve_piqp() creates a solver object and calls methods on
the object to obtain the solution. The solver object can be created
explitly using piqp() and provides more elaborate
facilities for updating problem data. This can be very efficient when
one is solving the same kind of problem over and over.
The above problem could be solved using the solver model object thus:
model <- piqp(P, c, A, b, G, h, x_lb = x_lb, x_ub = x_ub)
sol2 <- model$solve()
identical(sol, sol2)
#> [1] TRUEIndeed, this is exactly what solve_piqp() does. But this
interface allows us to update the settings, the bounds, etc. and resolve
the same problem more efficiently.
model$update(x_lb = c(0, 0))
sol3 <- model$solve()
cat(sprintf("(Solution status, description): = (%d, %s)\n",
sol3$status, status_description(sol3$status)))
#> (Solution status, description): = (-2, The problem is primal infeasible.)Setting the lower bounds made the problem infeasible.
## Try to give an inappropriate b value
model$update(b = c(5, 2))
#> Error in model$update(b = c(5, 2)): Update parameters not match original problem dimensionsThe error message correctly indicates that b has wrong
dimensions.
The methods exposed by the model object can be seen in the
documentation for the object piqp_model.
For example, we could query the problem dimensions.
PIQP supports dense and sparse problem formulations. For small and dense problems the dense interface is preferred since vectorized instructions and cache locality can be exploited more efficiently, but for sparse problems the sparse interface and result in significant speedups.
Either interface can be requested explicitly via the
backend parameter which can take on any value among
"dense", "sparse", or "auto", the
default. The last value will automatically switch to a sparse interface
if any of the supplied inputs (\(A\),
\(P\), or \(G\)) is a sparse matrix; otherwise it uses
the dense interface.
sparse_sol <- solve_piqp(P, c, A, b, G, h, x_lb, x_ub, backend = "sparse")
str(sparse_sol)
#> List of 15
#> $ status: int 1
#> $ x : num [1:2] -0.6 -0.8
#> $ y : num -11.8
#> $ z : num [1:2] 1.64e+01 1.01e-08
#> $ z_lb : num [1:2] 5.61e-10 0.00
#> $ z_ub : num [1:2] 3.2e-10 0.0
#> $ s : num [1:2] 5.19e-11 2.00e-01
#> $ s_lb : num [1:2] 0.4 Inf
#> $ s_ub : num [1:2] 1.6 Inf
#> $ zeta : num [1:2] -0.6 -0.8
#> $ lambda: num -11.8
#> $ nu : num [1:2] 1.64e+01 1.01e-08
#> $ nu_lb : num [1:2] 5.61e-10 0.00
#> $ nu_ub : num [1:2] 3.2e-10 0.0
#> $ info :List of 24
#> ..$ status_desc : chr "solved"
#> ..$ iter : num 6
#> ..$ rho : num 1e-10
#> ..$ delta : num 1e-10
#> ..$ mu : num 9.01e-10
#> ..$ sigma : num 1e-06
#> ..$ primal_step : num 0.99
#> ..$ dual_step : num 0.99
#> ..$ primal_inf : num 4.44e-11
#> ..$ primal_rel_inf : num 1.6
#> ..$ dual_inf : num 8.11e-10
#> ..$ dual_rel_inf : num 7.2
#> ..$ primal_obj : num 6.16
#> ..$ dual_obj : num 6.16
#> ..$ duality_gap : num 3.88e-09
#> ..$ duality_gap_rel : num 11.8
#> ..$ factor_retires : num 0
#> ..$ reg_limit : num 1e-10
#> ..$ no_primal_update: num 0
#> ..$ no_dual_update : num 0
#> ..$ setup_time : num 0
#> ..$ update_time : num 0
#> ..$ solve_time : num 0
#> ..$ run_time : num 0Suppose that we want to solve the following 2-dimensional quadratic programming problem:
\[ \begin{array}{ll} \text{minimize} & 3x_1^2 + 2x_2^2 - x_1 - 4x_2\\ \text{subject to} & -1 \leq x \leq 1, ~ x_1 = 2x_2 \end{array} \]
Since the solver expects the objective in the form \(\frac{1}{2}x^\top P x + c^\top x\), we define
\[ P = 2 \cdot \begin{bmatrix} 3 & 0 \\ 0 & 2\end{bmatrix} \mbox{ and } q = \begin{bmatrix} -1 \\ -4\end{bmatrix}. \]
We have one equality constraint and box constraints. This leads to the following straight-forward formulation.
P <- matrix(2 * c(3, 0, 0, 2), nrow = 2, ncol = 2)
c <- c(-1, -4)
A <- matrix(c(1, -2), ncol = 2)
b <- 0
x_lb <- rep(-1.0, 2)
x_ub <- rep(1.0, 2)
sol <- solve_piqp(P = P, c = c, A = A, b = b, x_lb = x_lb, x_ub = x_ub)
cat(sprintf("(Solution status, description): = (%d, %s)\n",
sol$status, sol$info$status_desc))
#> (Solution status, description): = (1, solved)
cat(sprintf("Objective: %f, solution: (x1, x2) = (%f, %f)\n", sol$info$primal_obj, sol$x[1], sol$x[2]))
#> Objective: -0.642857, solution: (x1, x2) = (0.428571, 0.214286)But we can also choose to move the upper box constraints into the inequalities.
G <- diag(2)
h <- c(1, 1)
sol <- solve_piqp(P = P, c = c, A = A, b = b, G = G, h = h,
x_lb = c(-1, -1), x_ub = c(Inf, Inf))
cat(sprintf("(Solution status, description): = (%d, %s)\n",
sol$status, sol$info$status_desc))
#> (Solution status, description): = (1, solved)
cat(sprintf("Objective: %f, solution: (x1, x2) = (%f, %f)\n", sol$info$primal_obj, sol$x[1], sol$x[2]))
#> Objective: -0.642857, solution: (x1, x2) = (0.428571, 0.214286)Or we can move both of them into the inequalities.
G <- Matrix::Matrix(c(1, 0, -1, 0, 0, 1, 0, -1), byrow = TRUE,
nrow = 4, sparse = TRUE)
h <- rep(1, 4)
sol <- solve_piqp(A = A, b = b, c = c, P = P, G = G, h = h)
cat(sprintf("(Solution status, description): = (%d, %s)\n",
sol$status, status_description(sol$status)))
#> (Solution status, description): = (1, Solver solved problem up to given tolerance.)
cat(sprintf("Objective: %f, solution: (x1, x2) = (%f, %f)\n", sol$info$primal_obj, sol$x[1], sol$x[2]))
#> Objective: -0.642857, solution: (x1, x2) = (0.428571, 0.214286)All of them will yield the same result.
PIQP has a number of parameters that control its behavior, including
verbosity, tolerances, etc.; see help on piqp_settings().
As an example, in the last problem, we can reduce the number of
iterations.
s <- solve_piqp(P = P, c = c, A = A, b = b, G = G, h = h,
settings = list(max_iter = 3)) ## Reduced number of iterations
cat(sprintf("(Solution status, description): = (%d, %s)\n",
s$status, s$info$status_desc))
#> (Solution status, description): = (-1, max iterations reached)
cat(sprintf("Objective: %f, solution: (x1, x2) = (%f, %f)\n", s$info$primal_obj, s$x[1], s$x[2]))
#> Objective: -0.642857, solution: (x1, x2) = (0.428294, 0.214147)Note the different status, which should always be checked in code.