coneproj (version 1.18)

qprog: 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)

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\).

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

Author

Mary C. Meyer and Xiyue Liao

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

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