qprog: Quadratic Programming

In coneproj: Primal or Dual Cone Projections with Routines for Constrained Regression

Quadratic Programming

Description

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.

Usage

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".

Value

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.

Author(s)

Mary C. Meyer and Xiyue Liao

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.

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")

coneproj documentation built on Oct. 15, 2023, 9:06 a.m.