用 scikit 学习线性模型约束系数的总和
Constraint the sum of coefficients with scikit learn linear model
我正在做一个有 1000 个系数的 LassoCV。 Statsmodels 似乎无法处理这么多系数。所以我正在使用 scikit 学习。 Statsmodel 允许 .fit_constrained("coef1 + coef2...=1")。这将系数的总和限制为 = 1。我需要在 Scikit 中执行此操作。我也将截距保持为零。
from sklearn.linear_model import LassoCV
LassoCVmodel = LassoCV(fit_intercept=False)
LassoCVmodel.fit(x,y)
如有任何帮助，我们将不胜感激。
如评论中所述：文档和来源并未表明 sklearn 支持此功能！
我刚刚尝试了使用现成的凸优化求解器的替代方法。这只是一种简单的类似原型的方法，可能不适合您的（未完全定义的）任务（样本大小？）。
一些评论：
- implementation/model-formulation很简单
- 这个问题比我想象的更难解决
- 求解器 ECOS 遇到一般问题
- 求解器 SCS 达到了很好的精度（与 sklearn 相比更差）
- 但是：调整迭代以提高精度会破坏求解器
- 问题将对 SCS 不可行！
- 基于 SCS + bigM 的公式（约束在 objective 中作为惩罚项发布）看起来可用；但可能需要调整
- 仅测试了开源求解器，商业求解器可能会好得多
进一步尝试：
- 解决巨大的问题（与稳健性和准确性相比，性能变得更加重要），（加速的）投影随机梯度方法看起来很有希望
代码
""" data """
from time import perf_counter as pc
import numpy as np
from sklearn import datasets
diabetes = datasets.load_diabetes()
A = diabetes.data
y = diabetes.target
alpha=0.1
print('Problem-size: ', A.shape)
def obj(x): # following sklearn's definition from user-guide!
return (1. / (2*A.shape[0])) * np.square(np.linalg.norm(A.dot(x) - y, 2)) + alpha * np.linalg.norm(x, 1)
""" sklearn """
print('\nsklearn classic l1')
from sklearn import linear_model
clf = linear_model.Lasso(alpha=alpha, fit_intercept=False)
t0 = pc()
clf.fit(A, y)
print('used (secs): ', pc() - t0)
print(obj(clf.coef_))
print('sum x: ', np.sum(clf.coef_))
""" cvxpy """
print('\ncvxpy + scs classic l1')
from cvxpy import *
x = Variable(A.shape[1])
objective = Minimize((1. / (2*A.shape[0])) * sum_squares(A*x - y) + alpha * norm(x, 1))
problem = Problem(objective, [])
t0 = pc()
problem.solve(solver=SCS, use_indirect=False, max_iters=10000, verbose=False)
print('used (secs): ', pc() - t0)
print(obj(x.value.flat))
print('sum x: ', np.sum(x.value.flat))
""" cvxpy -> sum x == 1 """
print('\ncvxpy + scs sum == 1 / 1st approach')
objective = Minimize((1. / (2*A.shape[0])) * sum_squares(A*x - y))
constraints = [sum(x) == 1]
problem = Problem(objective, constraints)
t0 = pc()
problem.solve(solver=SCS, use_indirect=False, max_iters=10000, verbose=False)
print('used (secs): ', pc() - t0)
print(obj(x.value.flat))
print('sum x: ', np.sum(x.value.flat))
""" cvxpy approach 2 -> sum x == 1 """
print('\ncvxpy + scs sum == 1 / 2nd approach')
M = 1e6
objective = Minimize((1. / (2*A.shape[0])) * sum_squares(A*x - y) + M*(sum(x) - 1))
constraints = [sum(x) == 1]
problem = Problem(objective, constraints)
t0 = pc()
problem.solve(solver=SCS, use_indirect=False, max_iters=10000, verbose=False)
print('used (secs): ', pc() - t0)
print(obj(x.value.flat))
print('sum x: ', np.sum(x.value.flat))
输出
Problem-size: (442, 10)
sklearn classic l1
used (secs): 0.001451024380348898
13201.3508496
sum x: 891.78869298
cvxpy + scs classic l1
used (secs): 0.011165673357417458
13203.6549995
sum x: 872.520510561
cvxpy + scs sum == 1 / 1st approach
used (secs): 0.15350853891775978
13400.1272148
sum x: -8.43795102327
cvxpy + scs sum == 1 / 2nd approach
used (secs): 0.012579569383536493
13397.2932976
sum x: 1.01207061047
编辑
只是为了好玩，我使用加速投影梯度（代码中的备注！）的方法实现了一个缓慢的非优化原型求解器。
尽管这里的行为很慢（因为没有优化），但这个方法应该可以更好地解决大问题（因为它是一阶方法）。应该很有潜力吧！
警告： 对某些人来说可能被视为高级数值优化:-)
编辑 2： 我忘了在投影上添加非负约束（sum(x) == 1 如果 x 可以是非负数！)。这使得求解变得更加困难（数值问题），而且很明显，应该使用那些快速的专用投影之一（我现在太懒了；我认为 n*log n algs 可用）。再次强调：这个 APG 求解器是一个原型，还没有为实际任务做好准备。
代码
""" accelerated pg -> sum x == 1 """
def solve_pg(A, b, momentum=0.9, maxiter=1000):
""" remarks:
algorithm: accelerated projected gradient
projection: proj on probability-simplex
-> naive and slow using cvxpy + ecos
line-search: armijo-rule along projection-arc (Bertsekas book)
-> suffers from slow projection
stopping-criterion: naive
gradient-calculation: precomputes AtA
-> not needed and not recommended for huge sparse data!
"""
M, N = A.shape
x = np.zeros(N)
AtA = A.T.dot(A)
Atb = A.T.dot(b)
stop_count = 0
# projection helper
x_ = Variable(N)
v_ = Parameter(N)
objective_ = Minimize(0.5 * square(norm(x_ - v_, 2)))
constraints_ = [sum(x_) == 1]
problem_ = Problem(objective_, constraints_)
def gradient(x):
return AtA.dot(x) - Atb
def obj(x):
return 0.5 * np.linalg.norm(A.dot(x) - b)**2
it = 0
while True:
grad = gradient(x)
# line search
alpha = 1
beta = 0.5
sigma=1e-2
old_obj = obj(x)
while True:
new_x = x - alpha * grad
new_obj = obj(new_x)
if old_obj - new_obj >= sigma * grad.dot(x - new_x):
break
else:
alpha *= beta
x_old = x[:]
x = x - alpha*grad
# projection
v_.value = x
problem_.solve()
x = np.array(x_.value.flat)
y = x + momentum * (x - x_old)
if np.abs(old_obj - obj(x)) < 1e-2:
stop_count += 1
else:
stop_count = 0
if stop_count == 3:
print('early-stopping @ it: ', it)
return x
it += 1
if it == maxiter:
return x
print('\n acc pg')
t0 = pc()
x = solve_pg(A, y)
print('used (secs): ', pc() - t0)
print(obj(x))
print('sum x: ', np.sum(x))
输出
acc pg
early-stopping @ it: 367
used (secs): 0.7714511330487027
13396.8642379
sum x: 1.00000000002
我很惊讶之前没有人在评论中说明这一点，但我认为你的问题陈述中存在概念上的误解。
让我们从 Lasso 估计器的 定义 开始，例如 Statistical Learning with Sparsity The Lasso and Generalizations 中给出的例子Hastie、Tibshirani 和 Wainwright：
Given a collection of N predictor-response pairs {(xi,yi)}, the
lasso finds the fit coefficients (β0,βi) to the least-square
optimization problem with the additional constraint that the L1-norm
of the vector of coefficients βi is less than or equal to t.
其中系数向量的L1范数是所有系数大小的总和。 如果你的系数都是正的，这正好解决了你的问题。
现在，这个 t 和 scikit-learn 中使用的 alpha 参数有什么关系？好吧，事实证明，根据拉格朗日对偶性，t 的每个值与 alpha.
的每个值之间存在一对一的对应关系
这意味着当您使用 LassoCV 时，由于您使用的是 alpha 的值范围，因此根据定义，您使用的是所有系数之和的允许值范围！
综上所述，所有系数之和等于 1 的条件等同于对 alpha.
的特定值使用 Lasso
我正在做一个有 1000 个系数的 LassoCV。 Statsmodels 似乎无法处理这么多系数。所以我正在使用 scikit 学习。 Statsmodel 允许 .fit_constrained("coef1 + coef2...=1")。这将系数的总和限制为 = 1。我需要在 Scikit 中执行此操作。我也将截距保持为零。
from sklearn.linear_model import LassoCV
LassoCVmodel = LassoCV(fit_intercept=False)
LassoCVmodel.fit(x,y)
如有任何帮助，我们将不胜感激。
如评论中所述：文档和来源并未表明 sklearn 支持此功能！
我刚刚尝试了使用现成的凸优化求解器的替代方法。这只是一种简单的类似原型的方法，可能不适合您的（未完全定义的）任务（样本大小？）。
一些评论：
- implementation/model-formulation很简单
- 这个问题比我想象的更难解决
- 求解器 ECOS 遇到一般问题
- 求解器 SCS 达到了很好的精度（与 sklearn 相比更差）
- 但是：调整迭代以提高精度会破坏求解器
- 问题将对 SCS 不可行！
- 基于 SCS + bigM 的公式（约束在 objective 中作为惩罚项发布）看起来可用；但可能需要调整
- 仅测试了开源求解器，商业求解器可能会好得多
进一步尝试：
- 解决巨大的问题（与稳健性和准确性相比，性能变得更加重要），（加速的）投影随机梯度方法看起来很有希望
代码
""" data """
from time import perf_counter as pc
import numpy as np
from sklearn import datasets
diabetes = datasets.load_diabetes()
A = diabetes.data
y = diabetes.target
alpha=0.1
print('Problem-size: ', A.shape)
def obj(x): # following sklearn's definition from user-guide!
return (1. / (2*A.shape[0])) * np.square(np.linalg.norm(A.dot(x) - y, 2)) + alpha * np.linalg.norm(x, 1)
""" sklearn """
print('\nsklearn classic l1')
from sklearn import linear_model
clf = linear_model.Lasso(alpha=alpha, fit_intercept=False)
t0 = pc()
clf.fit(A, y)
print('used (secs): ', pc() - t0)
print(obj(clf.coef_))
print('sum x: ', np.sum(clf.coef_))
""" cvxpy """
print('\ncvxpy + scs classic l1')
from cvxpy import *
x = Variable(A.shape[1])
objective = Minimize((1. / (2*A.shape[0])) * sum_squares(A*x - y) + alpha * norm(x, 1))
problem = Problem(objective, [])
t0 = pc()
problem.solve(solver=SCS, use_indirect=False, max_iters=10000, verbose=False)
print('used (secs): ', pc() - t0)
print(obj(x.value.flat))
print('sum x: ', np.sum(x.value.flat))
""" cvxpy -> sum x == 1 """
print('\ncvxpy + scs sum == 1 / 1st approach')
objective = Minimize((1. / (2*A.shape[0])) * sum_squares(A*x - y))
constraints = [sum(x) == 1]
problem = Problem(objective, constraints)
t0 = pc()
problem.solve(solver=SCS, use_indirect=False, max_iters=10000, verbose=False)
print('used (secs): ', pc() - t0)
print(obj(x.value.flat))
print('sum x: ', np.sum(x.value.flat))
""" cvxpy approach 2 -> sum x == 1 """
print('\ncvxpy + scs sum == 1 / 2nd approach')
M = 1e6
objective = Minimize((1. / (2*A.shape[0])) * sum_squares(A*x - y) + M*(sum(x) - 1))
constraints = [sum(x) == 1]
problem = Problem(objective, constraints)
t0 = pc()
problem.solve(solver=SCS, use_indirect=False, max_iters=10000, verbose=False)
print('used (secs): ', pc() - t0)
print(obj(x.value.flat))
print('sum x: ', np.sum(x.value.flat))
输出
Problem-size: (442, 10)
sklearn classic l1
used (secs): 0.001451024380348898
13201.3508496
sum x: 891.78869298
cvxpy + scs classic l1
used (secs): 0.011165673357417458
13203.6549995
sum x: 872.520510561
cvxpy + scs sum == 1 / 1st approach
used (secs): 0.15350853891775978
13400.1272148
sum x: -8.43795102327
cvxpy + scs sum == 1 / 2nd approach
used (secs): 0.012579569383536493
13397.2932976
sum x: 1.01207061047
编辑
只是为了好玩，我使用加速投影梯度（代码中的备注！）的方法实现了一个缓慢的非优化原型求解器。
尽管这里的行为很慢（因为没有优化），但这个方法应该可以更好地解决大问题（因为它是一阶方法）。应该很有潜力吧！
警告： 对某些人来说可能被视为高级数值优化:-)
编辑 2： 我忘了在投影上添加非负约束（sum(x) == 1 如果 x 可以是非负数！)。这使得求解变得更加困难（数值问题），而且很明显，应该使用那些快速的专用投影之一（我现在太懒了；我认为 n*log n algs 可用）。再次强调：这个 APG 求解器是一个原型，还没有为实际任务做好准备。
代码
""" accelerated pg -> sum x == 1 """
def solve_pg(A, b, momentum=0.9, maxiter=1000):
""" remarks:
algorithm: accelerated projected gradient
projection: proj on probability-simplex
-> naive and slow using cvxpy + ecos
line-search: armijo-rule along projection-arc (Bertsekas book)
-> suffers from slow projection
stopping-criterion: naive
gradient-calculation: precomputes AtA
-> not needed and not recommended for huge sparse data!
"""
M, N = A.shape
x = np.zeros(N)
AtA = A.T.dot(A)
Atb = A.T.dot(b)
stop_count = 0
# projection helper
x_ = Variable(N)
v_ = Parameter(N)
objective_ = Minimize(0.5 * square(norm(x_ - v_, 2)))
constraints_ = [sum(x_) == 1]
problem_ = Problem(objective_, constraints_)
def gradient(x):
return AtA.dot(x) - Atb
def obj(x):
return 0.5 * np.linalg.norm(A.dot(x) - b)**2
it = 0
while True:
grad = gradient(x)
# line search
alpha = 1
beta = 0.5
sigma=1e-2
old_obj = obj(x)
while True:
new_x = x - alpha * grad
new_obj = obj(new_x)
if old_obj - new_obj >= sigma * grad.dot(x - new_x):
break
else:
alpha *= beta
x_old = x[:]
x = x - alpha*grad
# projection
v_.value = x
problem_.solve()
x = np.array(x_.value.flat)
y = x + momentum * (x - x_old)
if np.abs(old_obj - obj(x)) < 1e-2:
stop_count += 1
else:
stop_count = 0
if stop_count == 3:
print('early-stopping @ it: ', it)
return x
it += 1
if it == maxiter:
return x
print('\n acc pg')
t0 = pc()
x = solve_pg(A, y)
print('used (secs): ', pc() - t0)
print(obj(x))
print('sum x: ', np.sum(x))
输出
acc pg
early-stopping @ it: 367
used (secs): 0.7714511330487027
13396.8642379
sum x: 1.00000000002
我很惊讶之前没有人在评论中说明这一点，但我认为你的问题陈述中存在概念上的误解。
让我们从 Lasso 估计器的 定义 开始，例如 Statistical Learning with Sparsity The Lasso and Generalizations 中给出的例子Hastie、Tibshirani 和 Wainwright：
Given a collection of N predictor-response pairs {(xi,yi)}, the lasso finds the fit coefficients (β0,βi) to the least-square optimization problem with the additional constraint that the L1-norm of the vector of coefficients βi is less than or equal to t.
其中系数向量的L1范数是所有系数大小的总和。 如果你的系数都是正的，这正好解决了你的问题。
现在，这个 t 和 scikit-learn 中使用的 alpha 参数有什么关系？好吧，事实证明，根据拉格朗日对偶性，t 的每个值与 alpha.
这意味着当您使用 LassoCV 时，由于您使用的是 alpha 的值范围，因此根据定义，您使用的是所有系数之和的允许值范围！
综上所述，所有系数之和等于 1 的条件等同于对 alpha.