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

Question

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

(Причина этого в том, что на самом деле мне нужно подогнать тысячи кривых одновременно, исходя из гигабайта данных, которые я читаю с диска, и, следовательно, это sloooooow).

Обратите внимание, что количество полиномиальных коэффициентов будет ограничено (возможно, 5-10), поэтому точное совпадение будет крайне маловероятным, но это нормально, поскольку я пытаюсь найти базовый шаблон в данных с много случайного шума. Я понимаю, как можно использовать генетический алгоритм для подгонки кривой к набору данных, но для этого требуется много проходов через выборочные данные, и, следовательно, это не практично для моего приложения.

Есть ли способ согласовать кривую с одним проходом данных, где состояние, которое должно поддерживаться от образца к образцу, минимально?

Я должен добавить, что характер данных заключается в том, что точки могут лежать в любом месте на оси X между 0,0 и 1,0, но значения Y всегда будут либо 1,0, либо 0,0.

Итак, в Java я ищу класс со следующим интерфейсом:

public interface CurveFit {
   public void addData(double x, double y);
   public List<Double> getBestFit(); // Returns the polynomial coefficients
}

Класс, реализующий это, не должен хранить большой объем данных в своих полях экземпляра, не более килобайта даже для миллионов точек данных. Это означает, что вы не можете просто хранить данные, поскольку вы заставляете их делать несколько проходов позже.

edit: Некоторые полагают, что найти оптимальную кривую за один проход может быть невозможно, однако оптимальная подгонка не требуется, так же близко, как мы можем получить ее за один проход.

Голыми костями подхода может быть, если у нас есть способ начать с кривой, а затем способ изменить ее, чтобы она немного приблизилась к новым точкам данных по мере их появления - фактически в форме градиентного спуска. Есть надежда, что с достаточным количеством данных (и данных будет много), мы получим довольно хорошую кривую. Возможно, это вдохновляет кого-то на решение.

Answer 1

Если вы не возражаете, что у вас получится прямолинейная «кривая», тогда вам понадобится всего шесть переменных для любого количества данных. Вот исходный код, который входит в мою будущую книгу; Я уверен, что вы можете понять, как работает класс DataPoint:

Interpolation.h:

#ifndef __INTERPOLATION_H
#define __INTERPOLATION_H

#include "DataPoint.h"

class Interpolation
{
private:
  int m_count;
  double m_sumX;
  double m_sumXX;  /* sum of X*X */
  double m_sumXY;  /* sum of X*Y */
  double m_sumY;
  double m_sumYY;  /* sum of Y*Y */

public:
  Interpolation();

  void addData(const DataPoint& dp);

  double slope() const;
  double intercept() const;

  double interpolate(double x) const;
  double correlate() const;
};

#endif // __INTERPOLATION_H

Interpolation.cpp:

#include <cmath>

#include "Interpolation.h"

Interpolation::Interpolation()
{
  m_count = 0;
  m_sumX = 0.0;
  m_sumXX = 0.0;
  m_sumXY = 0.0;
  m_sumY = 0.0;
  m_sumYY = 0.0;
}

void Interpolation::addData(const DataPoint& dp)
{
  m_count++;
  m_sumX += dp.getX();
  m_sumXX += dp.getX() * dp.getX();
  m_sumXY += dp.getX() * dp.getY();
  m_sumY += dp.getY();
  m_sumYY += dp.getY() * dp.getY();
}

double Interpolation::slope() const
{
  return (m_sumXY - (m_sumX * m_sumY / m_count)) /
    (m_sumXX - (m_sumX * m_sumX / m_count));
}

double Interpolation::intercept() const
{
  return (m_sumY / m_count) - slope() * (m_sumX / m_count);
}


double Interpolation::interpolate(double X) const
{
  return intercept() + slope() * X;
}


double Interpolation::correlate() const
{
  return m_sumXY / sqrt(m_sumXX * m_sumYY);
}

Answer 2

Да, это проекция. Для

y = X beta + error        

где члены в нижнем регистре - векторы, а X - матрица, у вас есть вектор решения

\hat{beta} = inverse(X'X) X' y

согласно странице OLS . Вы почти никогда не хотите вычислять это напрямую, а используйте декомпозиции LR, QR или SVD. Ссылки на литературу многочисленны.

Если ваша задача имеет только один параметр (и, следовательно, x также является вектором), это сводится к суммированию перекрестных произведений между y и x.

Answer 3

Мне кажется, я нашел ответ на свой вопрос на основе модифицированной версии этого кода. Для тех, кто заинтересован, мой код Java здесь .

Answer 4

Предполагая, что вы не знаете, какая точка должна принадлежать какой кривой, что-то вроде преобразования Хафа может обеспечить то, что вам нужно.

Преобразование Хафа - это метод, позволяющий идентифицировать структуру в наборе данных. Одно из них предназначено для компьютерного зрения, где оно позволяет легко определять линии и границы в поле зрения.

Преимущества для этой ситуации:

• Каждую точку нужно рассматривать только один раз
• Вам не нужно хранить структуру данных для каждой строки-кандидата, только одну (сложную, многомерную) структуру
• Обработка каждой строки проста
• Вы можете остановиться в любой точке и вывести набор удачных совпадений
• Вы никогда не удаляете какие-либо данные, поэтому они не зависят от случайного расположения ссылок
• Вы можете выбирать между точностью и требованиями к памяти
• Не ограничивается точными совпадениями, но также выделяет частичные совпадения.

Подход

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

Answer 5

Вам нужно решение переопределенной линейной системы. Популярными методами являются нормальные уравнения (обычно не рекомендуется), QR-факторизация и разложение по сингулярным числам (SVD). В Википедии есть достойные объяснения, Трефетен и Бау - это очень хорошо. Ваши варианты:

1. Реализация вне ядра с помощью нормальных уравнений. Для этого требуется произведение A'A, где A имеет намного больше строк, чем столбцов (поэтому результат очень мал). Матрица A полностью определяется местами выборки, поэтому вам не нужно ее хранить, поэтому вычисление A'A является относительно дешевым (очень дешевым, если вам не нужно использовать память для расположения узлов). Как только A'A вычислено, вы получаете решение за один проход через ваши входные данные, но метод может быть нестабильным.

2. Внедрите факторизацию QR вне ядра. Классический Грам-Шмидт будет самым быстрым, но вы должны быть осторожны со стабильностью.

3. Делайте это в ядре с распределенной памятью (если у вас есть доступное оборудование). Библиотеки, такие как PLAPACK и SCALAPACK, могут это сделать, производительность должна быть намного выше 1. Параллельная масштабируемость не фантастическая, но она подойдет, если это проблема размера, о которой вы даже подумаете в последовательном порядке.

4. Используйте итерационные методы для вычисления SVD. В зависимости от спектральных свойств вашей системы (возможно, после предварительной обработки) это может очень быстро сходиться и не требует хранения для матрицы (которая в вашем случае имеет 5-10 столбцов, каждый из которых является размером ваших входных данных. Хорошая библиотека это SLEPc , вам нужно только найти произведение матрицы Вандермонда с вектором (так что вам нужно только хранить местоположения выборки). Это очень масштабируемо параллельно.

Answer 6

Учитывая характер ваших данных:

точки могут находиться в любом месте на оси X между 0,0 и 1,0, но значения Y всегда будут либо 1,0, либо 0,0.

Тогда вам не нужен даже один проход, так как эти две линии будут проходить точно через каждую точку:

• X = [0,0 ... 1,0], Y = 0,0
• X = [0,0 ... 1,0], Y = 1,0

Два коротких отрезка, длина единицы и каждая точка попадают на одну линию или другую.

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

Answer 7

Почему бы не использовать кольцевой буфер некоторого фиксированного размера (скажем, последние 1000 точек) и сделать так, чтобы стандартные наименьшие квадраты на основе QR-разложения соответствовали буферизованным данным? После заполнения буфера, каждый раз, когда вы получаете новую точку, вы заменяете самую старую и заново подбираете. Таким образом, у вас есть ограниченный рабочий набор, который все еще имеет некоторую локальность данных, без каких-либо проблем обработки живого потока (без памяти).

Answer 8

Ограничиваете ли вы количество полиномиальных коэффициентов (то есть соответствует ли максимальная степень х в вашем полиноме)?

Если нет, то вам не нужен алгоритм «наилучшего соответствия» - вы всегда можете точно сопоставить N точек данных ТОЧНО для многочлена из N коэффициентов.

Просто используйте матрицы для решения N одновременных уравнений для N неизвестных (N коэффициентов многочлена).

Если вы ограничиваете максимальное количество коэффициентов, каково ваше максимальное значение?

После ваших комментариев и редактирования:

Вам нужен фильтр нижних частот, чтобы отфильтровывать шум, а не подгонять полином к шуму.