# How do I fit a constrained regression in R so all coefficients are positive and above 0 [duplicate]

I am trying to understand how to solve the below quadratic program:

This is my model Y=π1X1+π2X2+π3X3+ε,

My constraints are: - all weights πk>0 ; so each weight need to be above 0

My goal is minimize ∑(Yi−(π1Xi1+π2Xi2+π3Xi3))2

The code below show the first weight to be equal to 0 for meq=1, I am looking for a solution where all weights are above 0. I am not sure if meq should be 1 or 0. And why with meq =0, the unconstrained result match the constrained one?

require("quadprog")
X <- matrix(runif(300), ncol=3)
Y <- X %*% c(0.2,0.3,0.5) + rnorm(100, sd=0.2)
Dmat = t(X) %*% X
Amat = t(diag(3))
bvec = c(0,0,0)
dvec = t(Y) %*% X
solve.QP(Dmat = Dmat, dvec = dvec, Amat = Amat, bvec = bvec, meq = 0, factorized = F)


Results with meq=1

$solution [1] 0.0000000 0.4365420 0.5267267$unconstrained.solution
[1] 0.2052446 0.3196065 0.4612527


Results with meq=0

$solution [1] 0.1450482 0.2105224 0.5632482$unconstrained.solution
[1] 0.1450482 0.2105224 0.5632482

• What is meq? Is this some R-code specific variable? – Carl Nov 2 '16 at 17:34
• meq is part of the solve.QP function – user4797853 Nov 2 '16 at 17:36
• You might look at how good the fit is in each case. The best fit is likely the better answer, at least for the way you are making assumptions. However, your would have to compare the goodness of fit adjusted for the number of parameters used to obtain them. – Carl Nov 2 '16 at 17:41
• Thanks Carl, the thing is that I am looking for a solution where all betas are above 0, meq=0 seems to provide it but also seem to equal the unconstrained solution, which I do not understand why it's happening. – user4797853 Nov 2 '16 at 17:45
• It is possible that you are finding local minima. You might be better off using more global methods. One way is to use many different starting values of your parameters. Another is to use methods that find global solutions, for example, random search, which is rather slowly converging, and to search the entire parameter space. – Carl Nov 2 '16 at 17:55

First of all you should add set.seed(1) (or any other seed) at the start of your code, so that it's reproducible. Coming to your question:

From ?solve.QP:

meq     the first meq constraints are treated as equality constraints, all further as inequality constraints (defaults to 0).


Thus, when you set meq=1, you're telling solve.QP that the first constraint must not be a inequality, but an equality: in other words, you're saying that $\pi_1$ must be 0. And as a matter of fact, solve.QP returns you a solution where $\pi_1=0$.

solve.QP(Dmat = Dmat, dvec = dvec, Amat = Amat, bvec = bvec, meq=1, factorized = F)
$solution [1] 0.0000000 0.3582641 0.6358091$value
[1] -12.71193

$unconstrained.solution [1] 0.2024494 0.2646992 0.5357331  Anyway, since you only require$\pi_i\ge 0 \quad i=1,\dots,3$, you must set meq=0, i.e., you need to tell solve.QP that all constraints are inequality constraints. You could even avoid passing a value for meq=0 in the call to solve.QP, anyway, since meq=0 is the default. Now, since the unconstrained solution satisfies the constraints, then it's also the constrained solution, and that's why the two solutions are equal: solve.QP(Dmat = Dmat, dvec = dvec, Amat = Amat, bvec = bvec, meq=0, factorized = F)$solution
[1] 0.2024494 0.2646992 0.5357331

$value [1] -12.91278$unconstrained.solution
[1] 0.2024494 0.2646992 0.5357331


EDIT concerning the possibility of other solutions, and thus the need for global methods. This specific optimization problem belongs to the family of convex optimization problems. For this kind of problems, a local minimum is also a global minimum. Since you already found a local minimum, you don't need to do anything else. Note that for another, non-convex, optimization problem, this wouldn't be true anymore.

Two side notes:

1. you cannot impose strict inequality constraints with solve.QP, thus the right way to write your constraints is $\pi_i\ge 0 \quad i=1,\dots,3$, not $\pi_i\gt 0 \quad i=1,\dots,3$
2. I find that solve.QP from quadprog has a counterintuitive interface, and it can be less robust than alternatives. I suggest you have a look at lsei from package limSolve. You can find a worked example here. Since the interface is different, be sure to read the vignette carefully if you decide to use lsei.

EDIT: the OP seems to be expecially interested in the case of strict inequality constraints, thus I add a couple of well-known tricks to deal with them. First of all, solve the problem with not strict inequality constraints: if the constraints are not active in the solution of the constrained problem, i.e., if the solution satisfies the strict inequality constraints, then you're all set. For example, in the second case where meq=0, solve.QP, which uses non-strict inequality constraints, found a solution where all coefficients are strictly $>0$, i.e., the constraints are not active in the constrained solution.

If one or more constraints are active in the solution (for example, the solution found by solve.QP when meq=1), then a trick is to add a small positive $\epsilon$ to the constraint. For example, instead of the constraints $\pi_i\gt 0 \quad i=1,\dots,3$, you use the constraints $\pi_i\ge \epsilon \quad i=1,\dots,3$. The solution found in this case will solve the original strict inequality constraints, since $\epsilon>0$. The issue with this approach is that the solution will depend on $\epsilon$. However:

1. one could always solve a sequence of problems with decreasing $\epsilon$ and gain insight on the behavior of the solution with respect to $\epsilon$.
2. if the optimization problem models a real physical problem, then strict inequality constraints on physical variables usually require some "tolerance". For example, imagine you want to optimize the efficiency of some machine, whose efficiency is linked, among other things, to a clearance between a rotating shaft and its static casing (for example, this would be the case of rotor-stator seal clearances in a gas turbine). The distance $\delta$ between the rotating and the static part does need to be strictly positive (if rotor and stator touch, a disaster is guaranteed!). However, to be sure that rotor and stator don't touch, the clearance needs to be above some assembly/manufacturing tolerance: you cannot accept a solution with, say, $\delta=0.001$ microns, because you'll never be able to manufacture a design with such a small gap. You need to solve a problem with a constraint $\delta>\epsilon$ where, say, $\epsilon=100$ microns. Thus, when real physical systems are involved and it seems that strict inequality constraints are needed, usually what we really need are not strict inequality constraints with a "physically meaningful" tolerance.
• Thank you so much, why can't I impose πi > 0 ? – user4797853 Nov 2 '16 at 19:35
• Can I impose impose πi > 0 with lsei? – user4797853 Nov 2 '16 at 19:42
• no, you can't. Most (all those available in R?) linear and quadratic programming solvers don't allow strict constraints. The reason is that allowing for them often leads to ill-posed problems. For example, minimize $x^2$ subject to $x>0$: the minimizer ($x=0$) is on the boundary of the feasible set, but it's not in the feasible set. Another example: maximize $x$ subject to $x<3$. See here for more details. – DeltaIV Nov 2 '16 at 22:43
• I've modified my answer adding the usual strategies to deal with strict constraints. – DeltaIV Nov 2 '16 at 22:46
• This is extremely helpful. If I understand it correctly, I need to change bvec = c(0,0,0) to bvec = c(epsilon,epsilon,epsilon) and adjust meq=0. I have just one follow up if you don't mind. The previous commentator, Carl, mentioned that setting priors might help as well but I couldn't figure out how to adjust the code for prior. Do you know by any chance? In any case, thank you very much for all your help. – user4797853 Nov 3 '16 at 0:41