10
$\begingroup$

I am using the nnls() function from the nnls package in R to do a linear regression for regressors $x_i$ and observations $y$. The function delivers beta coefficients $\beta_i\geq{0}, \forall i$. However, is it possible to apply the constraints only to some regressors so that

$$\beta_k \geq 0 \quad k \in \{1...,10\}, k\neq i \\ \beta_i \in \mathbb{R} \quad i \in \{1...,10\}$$

given that I have 10 regressor variables?

The nnls package also offers the possibility to enforce some coefficients to be negative and others to be positive via the function nnnpls(). However, I only want the positive constraint for some parameters, the other ones can be either positive or negative (i.e., unconstrained).

$\endgroup$
3
  • 2
    $\begingroup$ Investigate the optim() function in R. This function will allow you to maximize a likelihood function given some constraints. $\endgroup$ Feb 6, 2015 at 12:36
  • $\begingroup$ Added other option via reformulating the objective as a quadratic programming problem - this should be faster than the other options mentioned here... $\endgroup$ Jun 28, 2019 at 10:00
  • $\begingroup$ I think the option I posted as a solution should be checked as the correct answer, as it's the fastest & doesn't require starting values to be given... $\endgroup$ May 10, 2023 at 23:34

3 Answers 3

6
$\begingroup$

Regardless of what your computing platform may be, you can trick any linearly constrained linear model solver into doing what you want without having to modify any code at all. Observe that the model

$$\mathbb{E}(Y) = X \beta$$

is equivalent to the model

$$\mathbb{E}(Y + X\gamma) = X(\beta + \gamma)$$

for any fixed coefficients $\gamma$, because the constant vector $X\gamma$ has been added to both sides. Moreover, fitting the second model will produce the same covariance matrix of estimates, etc., as the first because no change has been made in the errors whatsoever. Consequently, after estimating $\beta+\gamma$ as $\widehat{\beta+\gamma}$ using any linear procedure (not just nnls), the estimate of $\beta$ is obtained as

$$\hat\beta = \widehat{\beta+\gamma} -\gamma.$$

Exploit this flexibility to make sure that some of the coefficients of $\widehat{\beta+\gamma}$ will naturally be positive (and therefore not constrained by nnls). This can be done by estimating what those coefficients might be, performing the adjusted procedure, and checking that it has not constrained the corresponding estimates. If any have been constrained, increase the amount of adjustment and repeat until success is achieved.

I propose using an initial ordinary least squares fit to start the procedure. To be conservative, change the estimates by some small multiple $\rho$ of their standard errors. If iteration is needed, keep increasing $\rho$. The following code doubles $\rho$ at each iteration. The entire algorithm is contained within the short repeat block in the middle of the code example below.

As an example, I generated a problem with $200$ observations of seven variables (and included a constant term). The following tableau summarizes the results:

            AIntercept  AX1  AX2   AX3  AX4  AX5 AX6 AX7
True              -3.5 -2.5 -1.5 -0.50 0.50 1.50 2.5 3.5
OLS               -3.4 -2.5 -1.6 -0.52 0.44 1.49 2.4 3.5
NNLS               0.0  0.0  0.0  0.00 0.63 0.87 2.3 4.0
NNLS.0            -3.5 -2.7  0.0  0.00 0.75 1.56 2.3 3.6
Constraints        0.0  0.0  1.0  1.00 1.00 1.00 1.0 1.0
Bound              0.0  0.0  1.0  1.00 0.00 0.00 0.0 0.0

The true values of $\beta$ are listed first. The first half are negative, the second half positive. These are followed by their OLS estimates, their NNLS estimates, and the modified NNLS estimates wherein the first two coefficients (AIntercept and AX1) were not constrained to be positive. The last two lines summarize the constraints, printing "1.0" where positivity constraints were applied. The constraints actually imposed in the solution, using the same indicator method, appear last. In this case the procedure worked well.

#
# Describe a problem.
#
p <- 7                              # Number of variables
n <- 200                            # Number of observations
beta <- 0:p - p/2                   # True coefficients
constraints <- c(0, 0, rep(1, p-1)) # Positivity constraint indicator
#
# Generate data.
#
set.seed(17)
A <- cbind(Intercept=rep(1, n), 
           matrix(rnorm(p*n), n, dimnames=list(c(), paste0("X", 1:p))))
b <- A %*% beta + rnorm(n)
#
# OLS for reference.
#
fit.lm <- lm(b ~ A - 1)
#
# NNLS.
#
library(nnls)
fit.nnls <- nnls(A, b)
#
# NNLS with selective constraints.
#
rho <- 2 
coefficients <- coef(fit.lm)
se <- coef(summary(fit.lm))[, "Std. Error"]
repeat {
  beta.0 <- (coefficients - rho*se) * (1 - constraints)
  b.0 <- b - A %*% beta.0
  fit.nnls.0 <- nnls(A, b.0)
  if (all(constraints[fit.nnls.0$bound] == 1)) break #$
  if (rho > 1000) stop("Unable to find a solution.")
  rho <- rho*2
}
fit.nnls.0$x <- fit.nnls.0$x + beta.0
#
# Compare.
#
bound <- rep(0, p+1); bound[fit.nnls.0$bound] <- 1
print(rbind(True=beta, OLS=coef(fit.lm), NNLS=coef(fit.nnls),
      NNLS.0=coef(fit.nnls.0), 
      Constraints=constraints, Bound = bound), digits=2)
$\endgroup$
1
  • 3
    $\begingroup$ I know this comment goes against the guidelines, but this is a wonderful explanation. $\endgroup$
    – Mihai
    Apr 25, 2021 at 10:03
5
$\begingroup$

The fastest option (faster than the other solutions posted here) is to observe that nnls can be reformulated as a quadratic programming problem, and that when written as a quadratic programming problem, one can apply individual constraints. Such problems in R can be fit using the quadprog or the osqp package. As the osqp package is the fastest, and allows the covariate matrix to be either dense or sparse, I'll use that one here.

Using that package, here is a function that does a constrained least squares fit using vectors of lower and upper bounds lower and upper on the fitted coefficients and allowing X to be either dense or sparse:

constrainedLS_osqp <- function(y, X,
                               lower=rep(0, ncol(X)), upper=rep(Inf, ncol(X)),
                               x.start = NULL, y.start = NULL) {
  require(osqp)
  require(Matrix)
  XtX = crossprod(X,X)
  Xty = crossprod(X,y)
  
  settings = osqpSettings(verbose = FALSE, eps_abs = 1e-8, eps_rel = 1e-8, linsys_solver = 0L,
                          warm_start = FALSE)
  pff = .sparseDiagonal(ncol(X))
  model <- osqp(XtX, -Xty, pff, l=lower, u=upper, pars=settings)
  
  if (!is.null(x.start)) model$WarmStart(x=x.start, y=y.start)

  coefs = model$Solve()$x # fitted coefficients

  coefs = pmax(lower, pmin(coefs, upper) ) # fitted coefficients sometimes go very slightly outside constraint zone due to numerical inaccuracies in solver - this is fixed here via clipping

  return(coefs)
}

Alternatively, another option is to use glmnet's bigGlm function and set options lower.limits & upper.limits to specify the desired box constraints. That uses coordinate descent. It also allows your covariate matrix to be either dense or sparse.

$\endgroup$
14
  • 1
    $\begingroup$ PS if your covariate matrix would happen to be sparse it is more efficient to use the osqp quadratic programming solver in the R package of the same name... $\endgroup$ Apr 17, 2021 at 20:02
  • 1
    $\begingroup$ While I find @whuber's answer good for building intuition, I find your solution extremely helpful and efficient. $\endgroup$
    – Mihai
    Apr 25, 2021 at 15:00
  • 1
    $\begingroup$ @Mihai Have a try with the function constrainedLS_osqp that I included in my edited post above. Would that work for you? $\endgroup$ Apr 26, 2021 at 8:13
  • 1
    $\begingroup$ @Mihai One other function that might work for you is rdrr.io/cran/zetadiv/man/glm.cons.html $\endgroup$ Apr 26, 2021 at 8:38
  • 1
    $\begingroup$ Ha no, sorry that was a mistake - I had wanted to keep it sparse for efficiency reasons, made a mistake by quickly copying and pasting some bits together... If you put pff = t(.sparseDiagonal(ncol(X))) as an argument, does that work? X in any case can be either sparse or dense... $\endgroup$ Apr 28, 2021 at 21:05
2
$\begingroup$

The solution is constrOptim(), which allows for individual constraining of coefficients.

For example:

If I want to regress variable $y$ with the two regressors $x_1, x_2$, I can define the linear regression minimization problem (minimizing squared residuals) myself via

min.RSS <- function(data, par){
  with(data, sum((par[1]*x1 + par[2]*x2 - y)^2))
}

All data required is stored in a data frame, column one and two contains regressors $x_1, x_2$, and column three contains the observations $y$:

dat = data.frame(x1=c(1,2), x2=c(2,3), y=c(5,6))

Then the constrained optimization can be conducted via

myObj        = constrOptim(theta, f, grad, data, ui, ci)
coefficients = myObj$par

where theta contains the starting values for par[1] and par[2], f is the function which is to be minimized, grad the gradient of f or NULL, data contains the data frame dat, ui is the constraint matrix, and ci the constraint vector.

$\endgroup$
1
  • $\begingroup$ Disadvantage of this method is that it requires starting values; the method I posted below uses the quadratic programming solver osqp and doesn't need that & is computationally also much faster (it also supports sparse covariate matrices). $\endgroup$ Aug 3, 2022 at 17:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.