用scikit学习线性模型约束系数之和

Question

我正在做一个有1000个系数的LassoCV.Statsmodels似乎无法处理这么多系数.所以我正在使用scikit学习.Statsmodel允许.fit_constrained("coef1 + coef2 ... = 1").这将coefs的总和约束为= 1.我需要在Scikit中执行此操作.我也将拦截保持为零.

from sklearn.linear_model import LassoCV

LassoCVmodel = LassoCV(fit_intercept=False)
LassoCVmodel.fit(x,y)

任何帮助,将不胜感激.

Answer 1

正如评论中所提到的:文档和来源并不表示sklearn支持这一点!

我刚刚尝试了使用离架凸优化求解器的替代方案.它只是一个简单的类似原型的方法,它可能不适合您的(未完全定义的)任务(样本大小？).

一些评论:

• 实施/模型制定很容易
• 问题比我想象的更难解决
  • 解决者ECOS有一般麻烦
  • 求解器SCS达到良好的准确性(与sklearn相比更差)
  • 但是:调整迭代以提高准确性会破坏求解器
    • 问题对于SCS来说是不可行的!
  • 基于SCS + bigM的公式(约束在目标范围内作为惩罚条款发布)看起来可用; 但可能需要调整
  • 只测试了开源求解器,商业解算器可能要好得多

还有一些事情要尝试:

• 解决巨大的问题(与稳健性和准确性相比,性能变得更加重要),(加速)预测随机梯度方法看起来很有希望

码

""" 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:我忘了在投影上添加非负约束(如果x可以是非负的,则sum(x)== 1没有多大意义!).这使得解决更加困难(数值问题)并且显而易见的是,应当使用其中一个快速专用投影(我现在太懒了;我认为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