python - 分段线性回归的优化

Question

我正在尝试创建分段线性回归以最小化 MSE(最小平方误差)，然后直接使用线性回归。该方法应该使用动态规划来计算不同的分段大小和组的组合，以实现整体 MSE。我认为算法运行时间是 O(n²)，我想知道是否有办法将其优化为 O(nLogN)？

import numpy as np
from sklearn.metrics import mean_squared_error
from sklearn import linear_model
import pandas as pd
import matplotlib.pyplot as plt

x = [3.4, 1.8, 4.6, 2.3, 3.1, 5.5, 0.7, 3.0, 2.6, 4.3, 2.1, 1.1, 6.1, 4.8,3.8]
y = [26.2, 17.8, 31.3, 23.1, 27.5, 36.5, 14.1, 22.3, 19.6, 31.3, 24.0, 17.3, 43.2, 36.4, 26.1]
dataset = np.dstack((x,y))
dataset = dataset[0]
d_arg = np.argsort(dataset[:,0])
dataset = dataset[d_arg]

def calc_error(dataset):
    lr_model = linear_model.LinearRegression()
    x = pd.DataFrame(dataset[:,0])
    y = pd.DataFrame(dataset[:,1])
    lr_model.fit(x,y)
    predictions = lr_model.predict(x)
    mse = mean_squared_error(y, predictions)
    return mse

#n is the number of points , m is the number of groups, k is the minimum number of points in a group
#（15，5，3）returns 【【3，3，3，3，3】】
#（15，5，2） returns [[2,4,3,3,3],[3,2,4,2,4],[4,2,3,3,3]....]
def all_combination(n,m,k):
    result = []
    if n < k*m:
        print('There are not enough elements to split.')
        return
    combination_bottom = [k for q in range(m)] 
    #add greedy algorithm here?
    if n == k*m:
        result.append(combination_bottom.copy())
    else:
        combination_now = [combination_bottom.copy()]
        j = k*m+1
        while j < n+1:
            combination_last = combination_now.copy()
            combination_now = []
            for x in combination_last:
                for i in range (0, m):
                    combination_new = x.copy()
                    combination_new[i] = combination_new[i]+1
                    combination_now.append(combination_new.copy())
            j += 1
        else:
            for x in combination_last:
                for i in range (0, m):
                    combination_new = x.copy()
                    combination_new[i] = combination_new[i]+1
                    if combination_new not in result:
                        result.append(combination_new.copy())
            
    return result #2-d list

def calc_sum_error(dataset,cb):#cb = combination
    mse_sum = 0    
    for n in range(0,len(cb)):
        if n == 0:
            low = 0
            high = cb[0]
        else:
            low = 0
            for i in range(0,n):
                low += cb[i]
            high = low + cb[n]
        mse_sum += calc_error(dataset[low:high])
    return mse_sum
            
    
    
#k is the number of points as a group
def best_piecewise(dataset,k):
    lenth = len(dataset)
    max_split = lenth // k
    min_mse = calc_error(dataset)
    split_cb = []
    all_cb = []
    for i in range(2, max_split+1):
        split_result = all_combination(lenth, i, k)
        all_cb += split_result
        for cb in split_result:
            tmp_mse = calc_sum_error(dataset,cb)
            if tmp_mse < min_mse:
                min_mse = tmp_mse
                split_cb = cb
    return min_mse, split_cb, all_cb

min_mse, split_cb, all_cb = best_piecewise(dataset, 2)

print('The best split of the data is '+str(split_cb))
print('The minimum MSE value is '+str(min_mse))

x = np.array(dataset[:,0])
y = np.array(dataset[:,1])

plt.plot(x,y,"o")
for n in range(0,len(split_cb)):
    if n == 0:
        low = 0 
        high = split_cb[n]
    else:
        low = 0
        for i in range(0,n):
            low += split_cb[i]
        high = low + split_cb[n]
    x_tmp = pd.DataFrame(dataset[low:high,0])
    y_tmp = pd.DataFrame(dataset[low:high,1])
    lr_model = linear_model.LinearRegression()
    lr_model.fit(x_tmp,y_tmp)
    y_predict = lr_model.predict(x_tmp)
    plt.plot(x_tmp, y_predict, 'g-')

plt.show()

如果我在任何部分没有说清楚，请告诉我。

Answer 1

我花了一些时间才意识到 您描述的问题正是决策树回归器试图解决的问题。
不幸的是，最佳决策树的构建是 NP 难的，这意味着即使使用动态编程，您也无法将运行时间降低到 O(NlogN) 之类的东西。
好消息是您可以直接使用任何维护良好的决策树实现，DecisionTreeRegressor的 sklearn.tree例如，可以确定在 O(NlogN) 时间复杂度中获得最佳性能。要强制每组的最小点数，请使用 min_samples_leaf范围。您还可以控制其他几个属性，例如最大数量。与 max_leaf_nodes 的群组, 优化 w.r.t 使用 criterion 的不同损失函数等等
如果您好奇 Scikit-learn 的决策树与您的算法学习到的决策树相比如何(即代码中的 split_cb):

X = np.array(x).reshape(-1,1)
dt = DecisionTreeRegressor(min_samples_leaf=MIN_SIZE).fit(X,y)
split_cb = np.unique(dt.apply(X),return_counts=True)[1]

然后使用您使用的相同绘图代码。请注意，由于您的时间复杂度远高于 O(NlogN)*，因此您的实现通常会比 scikit-learn 的贪婪算法找到更好的分割。

[1] Hyafil, L. 和 Rivest, R. L. (1976)。构建最优二叉决策树是 np 完全的。信息处理快报，5(1)，15-17
*虽然我不确定您的实现的确切时间复杂度，但肯定比 O(N^2)、 all_combination(21,4,2) 更糟花了超过5分钟。