qprog function

Quadratic Programming

Given a positive definite $n$ by $n$ matrix $Q$ and a constant vector $c$ in $R^n$, the object is to find $\theta$ in $R^n$ to minimize $\theta'Q\theta - 2c'\theta$ subject to $A\theta \ge b$, for an irreducible constraint matrix $A$. This routine transforms into a cone projection problem for the constrained solution.

qprog(q, c, amat, b, face = NULL, msg = TRUE)

Arguments

• q: A $n$ by $n$ positive definite matrix.
• c: A vector of length $n$.
• amat: A $m$ by $n$ constraint matrix. The rows of amat must be irreducible.
• b: A vector of length $m$. Its default value is 0.
• face: A vector of the positions of edges, which define the initial face for the cone projection. For example, when there are $m$ cone edges, then face is a subset of $1,\ldots,m$. The default is face = NULL.
• msg: A logical flag. If msg is TRUE, then a warning message will be printed when there is a non-convergence problem; otherwise no warning message will be printed. The default is msg = TRUE

Details

To get the constrained solution to $\theta'Q\theta - 2c'\theta$ subject to $A\theta \ge b$, this routine makes the Cholesky decomposition of $Q$. Let $U'U = Q$, and define $\phi = U\theta$ and $z = U^{-1}c$, where $U^{-1}$ is the inverse of $U$. Then we minimize $||z - \phi||^2$, subject to $B\phi \ge 0$, where $B = AU^{-1}$. It is now a cone projection problem with the constraint cone $C$ of the form $\{\phi: B\phi \ge 0 \}$. This routine gives the estimation of $\theta$, which is $U^{-1}$ times the estimation of $\phi$.

The routine qprog dynamically loads a C++ subroutine "qprogCpp".

Returns

• df: The dimension of the face of the constraint cone on which the projection lands.

• thetahat: A vector minimizing $\theta'Q\theta - 2c'\theta$.

• steps: The number of iterations before the algorithm converges.

• xmat: The rows of the matrix are the edges of the face of the polar cone on which the residual of the projection onto the constraint cone lands.

• face: A vector of the positions of edges, which define the face on which the final projection lands on. For example, when there are $m$ cone edges, then face is a subset of $1,\ldots,m$.

References

Goldfarb, D. and A. Idnani (1983) A numerically stable dual method for solving strictly convex quadratic programs. Mathematical Programming 27, 1--33.

Fraser, D. A. S. and H. Massam (1989) A mixed primal-dual bases algorithm for regression under inequality constraints application to concave regression. Scandinavian Journal of Statistics 16, 65--74.

Fang,S.-C. and S. Puthenpura (1993) Linear Optimization and Extensions. Englewood Cliffs, New Jersey: Prentice Hall.

Silvapulle, M. J. and P. Sen (2005) Constrained Statistical Inference. John Wiley and Sons.

Meyer, M. C. (2013b) A simple new algorithm for quadratic programming with applications in statistics. Communications in Statistics 42(5), 1126--1139.

Liao, X. and M. C. Meyer (2014) coneproj: An R package for the primal or dual cone projections with routines for constrained regression. Journal of Statistical Software 61(12), 1--22.

Author(s)

Mary C. Meyer and Xiyue Liao

See Also

coneA

Examples

# load the cubic data set
    data(cubic)

# extract x
    x <- cubic$x

# extract y
    y <- cubic$y

# make the design matrix
    xmat <- cbind(1, x, x^2, x^3)

# make the q matrix
    q <- crossprod(xmat)

# make the c vector
    c <- crossprod(xmat, y)

# make the constraint matrix to constrain the regression to be increasing, nonnegative and convex
    amat <- matrix(0, 4, 4)
    amat[1, 1] <- 1; amat[2, 2] <- 1
    amat[3, 3] <- 1; amat[4, 3] <- 1
    amat[4, 4] <- 6
    b <- rep(0, 4)

# call qprog 
    ans <- qprog(q, c, amat, b)

# get the constrained fit of y
    betahat <- fitted(ans)
    fitc <- crossprod(t(xmat), betahat)

# get the unconstrained fit of y
    fitu <- lm(y ~ x + I(x^2) + I(x^3))

# make a plot to compare fitc and fitu
    par(mar = c(4, 4, 1, 1))
    plot(x, y, cex = .7, xlab = "x", ylab = "y")
    lines(x, fitted(fitu))
    lines(x, fitc, col = 2, lty = 4)
    legend("topleft", bty = "n", c("constr.fit", "unconstr.fit"), lty = c(4, 1), col = c(2, 1))
    title("Qprog Example Plot")