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

Question

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

Answer 1

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.