Sub-quadrati c алгоритм подбора кривой с двумя линиями

Question

Проблема состоит в том, чтобы найти наилучшее соответствие действительной двумерной кривой (заданной набором точек) с ломаной линией, состоящей из двух линий.

Подход грубой силы заключался бы в нахождении "левая" и "правая" линейные аппроксимации для каждой точки кривой и выбор пары с минимальной ошибкой. Я могу вычислить две линейные аппроксимации постепенно, перебирая точки кривой, но я не могу найти способ постепенного вычисления ошибки. Таким образом, этот подход приводит к квадратичной c сложности.

Вопрос в том, существует ли алгоритм, который обеспечит подквадратичную c сложность?

Второй вопрос: есть ли удобная библиотека C ++ для таких алгоритмов?

---

EDIT Для подгонки одной строкой есть формулы:

m = (Σxiyi - ΣxiΣyi/N) / (Σxi2 - (Σxi)2/N)
b = Σyi/N - m * Σxi/N

где m - наклон, а b - смещение линии. Наличие такой формулы для ошибки подбора решило бы проблему наилучшим образом.

Answer 1

Заявление об ограничении ответственности: мне не хочется разбираться в том, как это сделать на C ++, поэтому я буду использовать нотацию Python (numpy). Концепции полностью переносимы, поэтому у вас не должно возникнуть проблем с переводом обратно на язык по вашему выбору.

Допустим, у вас есть пара массивов x и y, содержащих точки данных, и что x монотонно увеличивается. Предположим также, что вы всегда будете выбирать точку разделения, которая оставляет не менее двух элементов в каждом разделе, поэтому уравнения разрешимы.

Теперь вы можете вычислить некоторые соответствующие величины:

N = len(x)

sum_x_left = x[0]
sum_x2_left = x[0] * x[0]
sum_y_left = y[0]
sum_y2_left = y[0] * y[0]
sum_xy_left = x[0] * y[0]

sum_x_right = x[1:].sum()
sum_x2_right = (x[1:] * x[1:]).sum()
sum_y_right = y[1:].sum()
sum_y2_right = (y[1:] * y[1:]).sum()
sum_xy_right = (x[1:] * y[1:]).sum()

Причина, по которой нам нужны эти величины (а это O(N) для инициализации), заключается в том, что вы можете использовать их напрямую для вычисления некоторых хорошо известных формул для параметров линейной регрессии. Например, оптимальные m и b для y = m * x + b задаются как

μ<sub>x</sub> = Σx<sub>i</sub>/N
μ<sub>y</sub> = Σy<sub>i</sub>/N
m = Σ(x<sub>i</sub> - μ<sub>x</sub>)(y<sub>i</sub> - μ<sub>y</sub>) / Σ(x<sub>i</sub> - μ<sub>x</sub>)<sup>2</sup>
b = μ<sub>y</sub> - m * μ<sub>x</sub>

Сумма квадратов ошибок определяется как

e = Σ(y<sub>i</sub> - m * x<sub>i</sub> - b)<sup>2</sup>

Их можно расширить с помощью простую алгебру в следующее:

m = (Σx<sub>i</sub>y<sub>i</sub> - Σx<sub>i</sub>Σy<sub>i</sub>/N) / (Σx<sub>i</sub><sup>2</sup> - (Σx<sub>i</sub>)<sup>2</sup>/N)
b = Σy<sub>i</sub>/N - m * Σx<sub>i</sub>/N
e = Σy<sub>i</sub><sup>2</sup> + m<sup>2</sup> * Σx<sub>i</sub><sup>2</sup> + N * b<sup>2</sup> - m * Σx<sub>i</sub>y<sub>i</sub> - b * Σy<sub>i</sub> + m * b * Σx<sub>i</sub>

Таким образом, вы можете l oop по всем возможностям и записать минимальное e:

for p in range(1, N - 3):
    # shift sums: O(1)
    sum_x_left += x[p]
    sum_x2_left += x[p] * x[p]
    sum_y_left += y[p]
    sum_y2_left += y[p] * y[p]
    sum_xy_left += x[p] * y[p]

    sum_x_right -= x[p]
    sum_x2_right -= x[p] * x[p]
    sum_y_right -= y[p]
    sum_y2_right -= y[p] * y[p]
    sum_xy_right -= x[p] * y[p]

    # compute err: O(1)
    n_left = p + 1
    slope_left = (sum_xy_left - sum_x_left * sum_y_left * n_left) / (sum_x2_left - sum_x_left * sum_x_left / n_left)
    intercept_left = sum_y_left / n_left - slope_left * sum_x_left / n_left
    err_left = sum_y2_left + slope_left * slope_left * sum_x2_left + n_left * intercept_left * intercept_left - slope_left * sum_xy_left - intercept_left * sum_y_left + slope_left * intercept_left * sum_x_left

    n_right = N - n_left
    slope_right = (sum_xy_right - sum_x_right * sum_y_right * n_right) / (sum_x2_right - sum_x_right * sum_x_right / n_right)
    intercept_right = sum_y_right / n_right - slope_right * sum_x_right / n_right
    err_right = sum_y2_right + slope_right * slope_right * sum_x2_right + n_right * intercept_right * intercept_right - slope_right * sum_xy_right - intercept_right * sum_y_right + slope_right * intercept_right * sum_x_right

    err = err_left + err_right
    if p == 1 || err < err_min
        err_min = err
        n_min_left = n_left
        n_min_right = n_right
        slope_min_left = slope_left
        slope_min_right = slope_right
        intercept_min_left = intercept_left
        intercept_min_right = intercept_right

Вероятно, есть другие упрощения, которые вы можете make, но этого достаточно, чтобы иметь алгоритм O(n).

Answer 2

На случай, если это поможет, вот какой-то код C, который я использовал для такого рода вещей. Это мало что добавляет к тому, что сказал Безумный физик.

Во-первых, формула. Если вы поместите линию y ^: x-> a * x + b через некоторые точки, тогда ошибка будет выражена следующим образом:

E = Sum{ sqr(y[i]-y^(x[i])) }/ N = Vy - Cxy*Cxy/Vx
where 
Vx is the variance of the xs
Vy that of the ys 
Cxy the covariance of the as and the ys

В приведенном ниже коде используется структура, которая содержит средние значения, дисперсии, ковариация и счет.

Функция moms_acc_pt () обновляет их, когда вы добавляете новую точку. Функция moms_line () возвращает a и b для строки и ошибку, как указано выше. Fmax (0,) в возврате используется в случае почти идеального совпадения, когда ошибка округления может отправить (математически неотрицательный) результат отрицательный.

Хотя можно было бы иметь функцию, которая удаляет точку из моментовT, мне легче решать, к каким моментамT добавить точку, делая копии, накапливая точки в копиях, получая линии, а затем сохраняя копию для той стороны, в которую точка лучше всего подходит, а оригинал для другой

typedef struct
{   int n;      // number points
    double  xbar,ybar;  // means of x,y
    double  Vx, Vy;     // variances of x,y
    double  Cxy;        // covariance of x,y
}   momentsT;

// update the moments to include the point x,y
void    moms_acc_pt( momentsT* M, double x, double y)
{   M->n += 1;
double  f = 1.0/M->n;
double  dx = x-M->xbar;
double  dy = y-M->ybar;
    M->xbar += f*dx;
    M->ybar += f*dy;
double  g = 1.0 - f;
    M->Vx   = g*(M->Vx  + f*dx*dx);
    M->Cxy  = g*(M->Cxy + f*dx*dy);
    M->Vy   = g*(M->Vy  + f*dy*dy);
}

// return the moments for the combination of A and B (assumed disjoint)
momentsT    moms_combine( const momentsT* A, const momentsT* B)
{
momentsT    C;
    C.n = A->n + B->n;
double  alpha = (double)A->n/(double)C.n;
double  beta = (double)B->n/(double)C.n;
    C.xbar = alpha*A->xbar + beta*B->xbar;
    C.ybar = alpha*A->ybar + beta*B->ybar;
double  dx = A->xbar - B->xbar;
double  dy = A->ybar - B->ybar;
    C.Vx = alpha*A->Vx + beta*B->Vx + alpha*beta*dx*dx;
    C.Cxy= alpha*A->Cxy+ beta*B->Cxy+ alpha*beta*dx*dy;
    C.Vy = alpha*A->Vy + beta*B->Vy + alpha*beta*dy*dy;
    return C;
}

// line is y^ : x -> a*x + b; return Sum{ sqr( y[i] - y^(x[i])) }/N
double  moms_line( momentsT* M, double* a, double *b)
{   *a = M->Cxy/M->Vx;
    *b = M->ybar - *a*M->xbar;
    return fmax( 0.0, M->Vy - *a*M->Cxy);
}