无论您的计算平台是什么，您都可以欺骗任何线性约束的线性模型求解器去做您想做的事情，而无需修改任何代码。观察模型

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

相当于模型

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

对于任何固定系数$\gamma$, 因为常数向量$X\gamma$已添加到双方。此外，拟合第二个模型将产生与第一个模型相同的估计协方差矩阵等，因为误差没有发生任何变化。因此，估计后$\beta+\gamma$作为$\widehat{\beta+\gamma}$使用任何线性程序（不仅仅是nnls），估计$\beta$获得为

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

利用这种灵活性来确保某些系数$\widehat{\beta+\gamma}$自然会是正数（因此不受 约束nnls）。这可以通过估计这些系数可能是什么、执行调整过程并检查它是否没有限制相应的估计来完成。如果有任何限制，请增加调整量并重复直到成功。

我建议使用初始普通最小二乘拟合来启动该过程。为了保守起见，将估计值更改一些小的倍数$\rho$他们的标准误。如果需要迭代，继续增加$\rho$. 以下代码加倍$\rho$在每次迭代中。整个算法包含在repeat下面代码示例中间的短块中。

例如，我产生了一个问题$200$七个变量的观察结果（包括一个常数项）。下表总结了结果：

            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

的真实价值观$\beta$首先列出。前半段为负，后半段为正。随后是他们的 OLS 估计、他们的 NNLS 估计和修改后的 NNLS 估计，其中前两个系数 (AIntercept和AX1) 不受约束为正。最后两行总结了约束，在应用了积极约束的地方打印“1.0”。解决方案中实际施加的约束，使用相同的指标方法，出现在最后。在这种情况下，程序运行良好。

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

解决方案是constrOptim()，它允许单独约束系数。

例如：

如果我想回归变量$y$与两个回归器$x_1, x_2$，我可以通过自己定义线性回归最小化问题（最小化平方残差）

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

所需的所有数据都存储在数据框中，第一列和第二列包含回归量$x_1, x_2$，第三列包含观察结果$y$：

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

然后可以通过以下方式进行约束优化

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

其中theta包含 和 的起始值par[1]，是要最小化的函数，par[2]或NULL的梯度，包含数据帧 dat，是约束矩阵和约束向量。fgradfdatauici

最快的选择（比此处发布的其他解决方案更快）是观察nnls可以重新表述为二次规划问题，并且当写为二次规划问题时，可以应用单独的约束。R 中的此类问题可以使用 thequadprog或osqp包来解决。由于osqp包是最快的，并且允许协变量矩阵要么密集要么稀疏，我将在这里使用那个。

lower使用该包，这是一个函数，它使用下限和上限向量以及upper拟合系数进行约束最小二乘拟合，并允许X密集或稀疏：

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)
  return(model$Solve()$x) # fitted coefficients
}