对“有谁知道这是为什么？ ”的相当反气候的答案是，根本没有人足够关心实施非负岭回归例程。主要原因之一是人们已经开始实施 非负弹性网络例程（例如这里和这里）。弹性网络包括岭回归作为一种特殊情况（基本上将 LASSO 部分设置为具有零权重）。这些作品相对较新，因此尚未纳入 scikit-learn 或类似的通用包中。您可能想向这些论文的作者询问代码。

编辑：

正如@amoeba 和我在评论中讨论的那样，实际实现相对简单。假设有以下回归问题：

$y = 2 x_1 - x_2 + \epsilon, \qquad \epsilon \sim N(0,0.2^2)$

其中和都是标准法线，例如：。请注意，我使用标准化的预测变量，因此我不必在之后进行标准化。为简单起见，我也不包括截距。我们可以立即使用标准线性回归解决这个回归问题。所以在 R 中它应该是这样的：$x_1$$x_2$$x_p \sim N(0,1)$

rm(list = ls()); 
library(MASS); 
set.seed(123);
N = 1e6;
x1 = rnorm(N)
x2 = rnorm(N)
y = 2 * x1 - 1 * x2 + rnorm(N,sd = 0.2)

simpleLR = lm(y ~ -1 + x1 + x2 )
matrixX = model.matrix(simpleLR); # This is close to standardised
vectorY = y
all.equal(coef(simpleLR), qr.solve(matrixX, vectorY), tolerance = 1e-7)  # TRUE

注意最后一行。几乎所有线性回归例程都使用 QR 分解来估计。我们想对岭回归问题使用相同的方法。此时阅读@whuber 的这篇文章；我们将完全执行此程序。简而言之，我们将用对角矩阵和我们的响应向量用个零这样我们就可以将原来的岭回归问题重新表达为 其中$\beta$$X$$\sqrt{\lambda}I_p$$y$$p$$(X^TX + \lambda I)^{-1} X^Ty$$(\bar{X}^T\bar{X})^{-1} \bar{X}^T\bar{y}$$\bar{}$象征着增强版。检查这些笔记中的幻灯片 18-19是否完整，我发现它们非常简单。因此，在 R 中，我们会喜欢以下内容：

myLambda = 100;  
simpleRR = lm.ridge(y ~ -1 + x1 + x2, lambda = myLambda)
newVecY = c(vectorY, rep(0, 2))
newMatX = rbind(matrixX, sqrt(myLambda) * diag(2))
all.equal(coef(simpleRR), qr.solve(newMatX, newVecY), tolerance = 1e-7)  # TRUE

它有效。好的，所以我们得到了岭回归部分。不过，我们可以用另一种方式解决，我们可以将其表述为一个优化问题，其中残差平方和是成本函数，然后针对它进行优化，即。。果然我们可以做到：$\displaystyle \min_{\beta} || \bar{y} - \bar{X}\beta||_2^2$

myRSS <- function(X,y,b){ return( sum( (y - X%*%b)^2 ) ) }
bfgsOptim = optim(myRSS, par = c(1,1), X = newMatX, y= newVecY, 
                  method = 'L-BFGS-B')
all.equal(coef(simpleRR), bfgsOptim$par, check.attributes = FALSE, 
          tolerance = 1e-7) # TRUE

正如预期的那样再次起作用。所以现在我们只需要：其中。这只是相同的优化问题，但受到约束，因此解决方案是非负的。$\displaystyle \min_{\beta} || \bar{y} - \bar{X}\beta||_2^2$$\beta \geq 0$

bfgsOptimConst = optim(myRSS, par = c(1,1), X=newMatX, y= newVecY, 
                       method = 'L-BFGS-B', lower = c(0,0))
all(bfgsOptimConst$par >=0)  # TRUE
(bfgsOptimConst$par) # 2.000504 0.000000

这表明原始的非负岭回归任务可以通过将其重新表述为一个简单的约束优化问题来解决。一些警告：

1. 我使用（实际上）归一化的预测变量。您需要自己考虑标准化。
2. 截距的非标准化也是如此。
3. 我使用了optim's L-BFGS-B参数。它是接受边界的最普通的 R 求解器。我相信你会找到几十个更好的求解器。
4. 通常，约束线性最小二乘问题被提出为二次优化任务。这对这篇文章来说有点过头了，但请记住，如果需要，您可以获得更快的速度。
5. 如评论中所述，您可以跳过岭回归作为增强线性回归部分，直接将岭成本函数编码为优化问题。这会简单得多，而且这篇文章要小得多。为了论证，我也附上了第二个解决方案。
6. 我在 Python 中不是完全会话式的，但基本上你可以通过使用 NumPy 的linalg.solve和 SciPy 的优化函数来复制这项工作。
7. 要选择超参数等，您只需执行您在任何情况下都会执行的常规 CV 步骤；没有什么变化。$\lambda$

第 5 点的代码：

myRidgeRSS <- function(X,y,b, lambda){ 
                return( sum( (y - X%*%b)^2 ) + lambda * sum(b^2) ) 
              }
bfgsOptimConst2 = optim(myRidgeRSS, par = c(1,1), X = matrixX, y = vectorY,
                        method = 'L-BFGS-B', lower = c(0,0), lambda = myLambda)
all(bfgsOptimConst2$par >0) # TRUE
(bfgsOptimConst2$par) # 2.000504 0.000000

R 包 glmnet 实现了弹性网络，因此 lasso 和 ridge 允许这样做。使用参数lower.limits和upper.limits，您可以分别为每个权重设置最小值或最大值，因此如果将下限设置为 0，它将执行非负弹性网（套索/脊）。

还有一个 python 包装器https://pypi.python.org/pypi/glmnet/2.0.0

回想一下我们正在尝试解决的问题：

$$\text{minimize}_{x}\,\,\,\,\left\Vert Ax-y\right\Vert _{2}^{2}+ \lambda \| x \|^2_2 \,\,\,\,\text{s.t. }x>0$$

相当于：

$$\text{minimize}_{x}\,\,\,\,\left\Vert Ax-y\right\Vert _{2}^{2}+ \lambda x^\top I x \,\,\,\,\text{s.t. }x>0$$

再加上一些代数：

$$\text{minimize}_{x}\,\,\,\,x^{T}\left(A^{T}A+ \lambda I \right)x+\left(-2A^{T}y\right)^{T}x\,\,\,\,\text{s.t. }x>0$$

伪python中的解决方案很简单：

Q = A'A + lambda*I
c = - A'y
x,_ = scipy.optimize.nnls(Q,c)

---

请参阅：正则化器进行稀疏非负最小二乘？$K$$x^\top R_k x$

以获得更一般的答案。

我知道这是一个老歌，但我在这里举了一个例子，说明如何在 Python 中通过使用ElasticNet（近似）或通过nnls().

基本上，对于ElasticNet，您可以使用：

from sklearn import linear_model as lm

eln = lm.ElasticNet(l1_ratio=0, fit_intercept=False)
act_alphas, coefs, dual_gaps = eln.path(X, y, alphas=alphas, positive=True)

但ElasticNet请注意，对于较大的 alpha 值会发出很多警告，因为它不会收敛。

通过nnls()：

from scipy.optimize import nnls

def nnls_ridge(X, y, alpha):
    """return non-negative Ridge coefficients"""
    p = X.shape[1]
    Xext = np.vstack((X, alpha * np.eye(p)))
    yext = np.hstack((y, np.zeros(p)))
    coefs, _ = nnls(Xext, yext)
    return coefs

笔记

（遵循这篇优秀文章的注释）。

虽然在不受约束的环境中，Ridge 的 OLS 公式是正确的：

$$(X_{*}^\prime X_{*})\beta = X_{*}^\prime y_{*}.$$

相当于：

$$(X^\prime X + \lambda I)\beta = X^\prime y.$$

在受约束的设置（例如bounds = (0, np.inf)或使用nnls()）中，两者不等价，因为已经通过（无约束的）LS 投影$(X^\prime X + \lambda I)$$p \times p$

特别是，我最初的尝试nnls_ridge()是错误的：

def nnls_ridge_wrong(X, y, alpha):
    p = X.shape[1]
    Xmod = X.T @ X + alpha * np.eye(p)
    ymod = X.T @ y
    coefs, _ = nnls(Xmod, ymod)
    return coefs

示例（使用本文开头 GitHub 链接中描述的功能）：

X = np.arange(6).reshape((-1, 2))
y = np.array([8,1,1])
a = 0.001

>>> ols_ridge(X, y, a)
array([-8.56345838,  6.81837427])

>>> nn_ridge_path_via_elasticnet(X, y, [a])[1].T
array([[0.        , 0.45708041]])

>>> nnls_ridge(X, y, a)
array([0.        , 0.45714284])

>>> nnls_ridge_wrong(X, y, a)
array([0.        , 0.37663842])

对于和，其中所有分量均为正 ，我们看到（左下图对）存在一些巨大差异：$X_{12,10}$$y_{12}$$N(0,1)$$R^2$nnls_ridge_wrong