分段最小二乘法

Question

给出一个算法，它将平面（x_1，y_1），（x_2，y_2），...，（x_n，y_n）中的点序列和整数k作为输入，并返回最佳分段线性函数f由至多k个最小化平方误差的块组成。你可以假设你有机会获得一种算法，它通过Θ（n）时间内的一组n个点来计算一个段的平方和误差。解决方案应该使用O（n^2k）时间和O（nk）空间。

任何人都可以帮我解决这个问题吗？非常感谢！

Answer 1

（这是为时已晚，你的功课，但希望它有助于反正。）
首先是在Python/numpy的动态编程仅k = 4， 帮助您了解动态规划工作;一旦你明白了，为任何k写一个循环应该很容易。
另外，Cost[]是2d矩阵，空间O（n^2）; 看到音符在结束用于获得下至空间为O（n k）的

#!/usr/bin/env python 
""" split4.py: min-cost split into 4 pieces, dynamic programming k=4 """ 

from __future__ import division 
import numpy as np 

__version__ = "2014-03-09 mar denis" 

#............................................................................... 
def split4(Cost, verbose=1): 
    """ split4.py: min-cost split into 4 pieces, dynamic programming k=4 
     min Cost[0:a] + Cost[a:b] + Cost[b:c] + Cost[c:n] 
     Cost[a,b] = error in least-squares line fit to xy[a] .. xy[b] *including b* 
     or error in lsq horizontal lines, sum (y_j - av y) ^2 for each piece -- 

      o-- 
     o- 
      o--- 
       o---- 
     | | | | 
     0 2 5 9 

     (Why 4 ? to walk through step by step, then put in a loop) 
    """ 
     # speedup: maxlen 2 n/k or so 

    Cost = np.asanyarray(Cost) 
    n = Cost.shape[1] 
     # C2 C3 ... costs, J2 J3 ... indices of best splits 
    J2 = - np.ones(n, dtype=int) # -1, NaN mark undefined/bug 
    C2 = np.ones(n) * np.NaN 
    J3 = - np.ones(n, dtype=int) 
    C3 = np.ones(n) * np.NaN 

     # best 2-splits of the left 2 3 4 ... 
    for nleft in range(1, n): 
     J2[nleft] = j = np.argmin([ Cost[0,j-1] + Cost[j,nleft] for j in range(1, nleft+1)]) + 1 
     C2[nleft] =     Cost[0,j-1] + Cost[j,nleft] 
      # an idiom for argmin j, min value c together 

     # best 3-splits of the left 3 4 5 ... 
    for nleft in range(2, n): 
     J3[nleft] = j = np.argmin([ C2[j-1] + Cost[j,nleft] for j in range(2, nleft+1)]) + 2 
     C3[nleft] =     C2[j-1] + Cost[j,nleft] 

     # best 4-split of all n -- 
    j4 = np.argmin([ C3[j-1] + Cost[j,n-1] for j in range(3, n)]) + 3 
    c4 =    C3[j4-1] + Cost[j4,n-1] 
    j3 = J3[j4] 
    j2 = J2[j3] 
    jsplit = np.array([ 0, j2, j3, j4, n ]) 
    if verbose: 
     print "split4: len %s pos %s cost %.3g" % (np.diff(jsplit), jsplit, c4) 
     print "split4: J2 %s C2 %s" %(J2, C2) 
     print "split4: J3 %s C3 %s" %(J3, C3) 
    return jsplit 

#............................................................................... 
if __name__ == "__main__": 
    import random 
    import sys 
    import spread 

    n = 10 
    ncycle = 2 
    plot = 0 
    seed = 0 

     # run this.py a=1 b=None c=[3] 'd = expr' ... in sh or ipython 
    for arg in sys.argv[1:]: 
     exec(arg) 
    np.set_printoptions(1, threshold=100, edgeitems=10, linewidth=100, suppress=True) 
    np.random.seed(seed) 
    random.seed(seed) 

    print "\n", 80 * "-" 
    title = "Dynamic programming least-square horizontal lines %s n %d seed %d" % (
      __file__, n, seed) 
    print title 

    x = np.arange(n + 0.) 
    y = np.sin(2*np.pi * x * ncycle/n) 
     # synthetic time series ? 
    print "y: %s av %.3g variance %.3g" % (y, y.mean(), np.var(y)) 

    print "Cost[j,k] = sum (y - av y)^2 --" # len * var y[j:k+1] 
    Cost = spread.spreads_allij(y) 
    print Cost # .round().astype(int) 

    jsplit = split4(Cost) 
     # split4: len [3 2 3 2] pos [ 0 3 5 8 10] 

    if plot: 
     import matplotlib.pyplot as pl 
     title += "\n lengths: %s" % np.diff(jsplit) 
     pl.title(title) 
     pl.plot(y) 
     for js, js1 in zip(jsplit[:-1], jsplit[1:]): 
      if js1 <= js: continue 
      yav = y[js:js1].mean() * np.ones(js1 - js + 1) 
      pl.plot(np.arange(js, js1 + 1), yav) 
     # pl.legend() 
     pl.show() 

然后，下面的代码做Cost[]对于水平线只，斜率0; 将它延伸到任何斜率的线段，在时间O（n）中，作为练习留下。

""" spreads(all y[:j]) in time O(n) 

define spread(y[]) = sum (y - average y)^2 
    e.g. spread of 24 hourly temperatures y[0:24] i.e. y[0] .. y[23] 
    around a horizontal line at the average temperature 
     (spread = 0 for constant temperature, 
     24 c^2 for constant + [c -c c -c ...], 
     24 * variance(y)) 

    How fast can one compute all 24 spreads 
    1 hour (midnight to 1 am), 2 hours ... all 24 ? 

    A simpler problem: compute all 24 averages in time O(n): 
     N = np.arange(1, len(y)+1) 
     allav = np.cumsum(y)/N 
      = [ y0, (y0 + y1)/2, (y0 + y1 + y2)/3 ...] 
    An identity: 
     spread(y) = sum(y^2) - n * (av y)^2 
    Voila: the code below, all spreads() in time O(n). 

    Exercise: extend this to spreads around least-squares lines 
    fit to [ y0, [y0 y1], [y0 y1 y2] ... ], not just horizontal lines. 
""" 

from __future__ import division 
import sys 
import numpy as np 


#............................................................................... 
def spreads(y): 
    """ [ spread y[:1], spread y[:2] ... spread y ] in time O(n) 
     where spread(y[]) = sum (y - average y)^2 
      = n * variance(y) 
    """ 
    N = np.arange(1, len(y)+1) 
    return np.cumsum(y**2) - np.cumsum(y)**2/N 

def spreads_allij(y): 
    """ -> A[i,j] = sum (y - av y)^2, spread of y around its average 
     for all y[i:j+1] 
     time, space O(n^2) 
    """ 
    y = np.asanyarray(y, dtype=float) 
    n = len(y) 
    A = np.zeros((n,n)) 
    for i in range(n): 
     A[i,i:] = spreads(y[i:]) 
    return A 

到目前为止，我们有一个n×n的成本矩阵，空间O（n^2）。 要踏踏实实地空间O（NK）， 在Cost[i,j]格局仔细一看访问的DYN-PROG代码：

for nleft .. to n: 
    Cost_nleft = Cost[j,nleft ] -- time nleft or nleft^2 
    for k in 3 4 5 ...: 
     min [ C[k-1, j-1] + Cost_nleft[j] for j .. to nleft ] 

这里Cost_nleft是完整的N×N的成本矩阵的一行，〜N段，根据需要生成。 这可以在时间O（n）线段完成。 但是，如果“一个段通过一组n点需要O（n）时间”的错误， 似乎我们到时间O（n^3）。评论任何人？

Answer 2

如果您可以对n^2中的某个细分区域做最小平方，可以通过动态编程轻松完成n^2 k^2中的任何操作。您可能只能将其优化为单个k。