Присоединение данных к полиному третьей степени

В настоящее время я пишу программу на С++, где у меня есть векторы независимых и зависимых данных, которые я хотел бы поместить в кубическую функцию. Однако у меня возникают проблемы с созданием полинома, который может соответствовать моим данным.

Часть проблемы заключается в том, что я не могу использовать различные числовые пакеты, такие как GSL (длинный рассказ); возможно, это может быть даже излишним для моего дела. Мне не нужно очень обобщенное решение для финализации наименьших квадратов. Я специально хочу подгонять мои данные к кубической функции. У меня есть доступ к библиотеке векторов Sony, которая поддерживает матрицы 4x4 и может, среди прочего, вычислять их обратные.

Во время прототипирования в Scilab я использовал такую функцию, как:

function p = polyfit(x, y, n)
    m = length(x);
    aa = zeros(m, n+1)
    aa(:,1) = ones(m,1)
    for k = 2:n+1
        aa(:,k) = x.^(k-1)
    end
    p = aa\y
endfunction

К сожалению, это не соответствует моей текущей среде. Приведенный выше пример должен поддерживать матрицу размеров M x N + 1. В моем случае это M x 4, где M зависит от того, сколько у меня данных выборки. Там также проблема левого деления. Мне понадобится библиотека матриц, поддерживающая обратную матрицу произвольных размеров.

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

Ответ 1

Да, мы можем ограничить проблему вычислением "матрицей 4x4". Наименьшая квадратичная подгонка кубической, даже для M данных точек, требует только решения четырех линейных уравнений с четырьмя неизвестными. Предполагая, что все x-координаты различны, матрица коэффициентов обратима, поэтому в принципе систему можно решить, инвертируя матрицу коэффициентов. Предположим, что M больше 4, как это обычно бывает для приемов наименьших квадратов.

Здесь написана запись для Maple, Установка кубика в данные, которая почти полностью скрывает детали того, что решается. Минимальные критерии первого порядка (первые производные по коэффициентам как параметры ошибки квадратов квадратов) дают нам четыре линейных уравнения, часто называемые нормальными уравнениями .

Вы можете "собрать" эти четыре уравнения в коде, а затем применить свою обратную матрицу или более сложную стратегию решения. Очевидно, что вам нужно хранить данные в определенной форме. Одна из возможностей состоит в двух линейных массивах: одна для x-координат и одна для y-координат, а длина M - количество точек данных.

NB: Я собираюсь обсудить эту матричную сборку в терминах индексов массива на основе 1. Полиномиальные коэффициенты на самом деле являются одним из приложений, в которых индексы основанных на 0 оснований делают вещи более чистыми и более простыми, но переписывание их на C или любой другой язык, который поддерживает индексы на основе 0, остается в качестве упражнения для читателя.

Линейная система нормальных уравнений наиболее легко выражается в матричной форме, обращаясь к массиву Mx4 A, чьи записи являются полномочиями x-координатных данных:

A (i, j) = x-координата i-й точки данных, поднятой до мощности j-1

Пусть A 'обозначает транспонирование A, так что A'A является матрицей 4x4.

Если d - столбец длины M, содержащий y-координаты точек данных (в данном порядке), то система нормальных уравнений такова:

A'A u = A 'd

где u = [p0, p1, p2, p3] '- столбец коэффициентов для кубического многочлена с минимальными квадратами:

P (x) = p0 + p1 * x + p2 * x ^ 2 + p3 * x ^ 3

Ваши возражения, похоже, сосредоточены на трудности хранения и/или манипулирования массивом Mx4 A или его транспонированием. Поэтому мой ответ будет сосредоточен на том, как собрать матрицу A'A и столбец A'd без явного хранения всего A за один раз. Другими словами, мы будем делать указанные матрично-матричные и матрично-векторные умножения неявно, чтобы получить систему 4x4, которую вы можете решить:

C u = f

Если вы думаете о том, что запись C (i, j) является произведением i-й строки A 'с j-м столбцом A, плюс тот факт, что i-я строка A' на самом деле является просто транспонированием i-го столбца A, должно быть ясно, что:

C (i, j) = SUM x ^ (i + j-2) по всем точкам данных

Это, безусловно, одно место, где изложение будет упрощено с помощью индексов на основе 0!

Возможно, имеет смысл накапливать записи для матрицы C, которые зависят только от значения я + j, т.е. так называемой ганкелевой матрицы, в линейный массив длины 7 такой, что:

W (k) = SUM x ^ k по всем точкам данных

где k = 0,..., 6. Матрица C 4x4 имеет "полосатую" структуру, которая означает, что появляются только эти семь значений. Перейдя по списку x-координат точек данных, вы можете накапливать соответствующие вклады каждой мощности каждой точки данных в соответствующей записи W.

Аналогичную стратегию можно использовать для сборки столбца f = A 'd, а именно для петли над точками данных и накопления следующих четырех сумм:

f (k) = SUM (x ^ k) * y во всех точках данных

где k = 0,1,2,3. [Конечно, в приведенных выше суммах значения x, y являются координатами для общей точки данных.]

Предостережения: Это удовлетворяет цели работы только с матрицей 4x4. Однако обычно пытаются избежать явного формирования матрицы коэффициентов для нормальных уравнений, потому что эти матрицы часто используются в численном анализе как плохо обусловленные. В частности, случаи, когда x-координаты расположены близко друг от друга, могут вызвать трудности при попытке решить систему путем инвертирования матрицы коэффициентов.

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

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

Ответ 2

Здесь - это страница, с которой я работаю, хотя сама эта страница напрямую не касается вашего вопроса. Резюме моего ответа будет:

Если вы не можете напрямую работать с матрицами Nx4, тогда эти матрицы вычислений "вручную", пока у вас не будет проблемы с чем-то, у которого есть только 4x4 или меньше матриц. В этом ответе я расскажу, как выполнять конкретные вычисления в матрице, которые вам нужны "вручную".

Предположим, что у вас есть куча точек данных (x1,y1)...(xn,yn), и вы ищете кубическое уравнение y = ax^3 + bx^2 + cx + d, которое лучше всего подходит для этих точек.

Затем, следуя ссылке выше, вы должны написать это уравнение:

Я напишу A, x и B для этих матриц. Затем, следуя приведенной выше ссылке, вы должны умножить на транспонирование A, которая даст вам матрицу 4x4 A ^T * A, которую вы можете инвертировать. В уравнениях приведен следующий план:

A * x = B.................... [с чего мы начали с]

(A ^T * A) * x = A ^T * B..... [умножить на A ^T]

x = (A ^T * A) ^-1 * A ^T * B... [умножить на обратное к A ^T * A]

Вы сказали, что довольны инвертированием матриц 4x4, поэтому, если мы сможем закодировать способ получить эти матрицы без фактического использования объектов матрицы, мы должны быть в порядке.

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

Напишите явное уравнение для каждого из 16 элементов матрицы 4x4. Запись (i,j) -го (я начинаю с (0,0)) дается выражением x ₁ⁱ * x ₁^j + x ₂ⁱ * x ₂^j +... + x _Nⁱ * x _{Nсуб > ^J.}
Инвертируйте эту матрицу 4x4, используя вашу матричную библиотеку. То есть (A ^T * A) ^-1.
Теперь нам нужно A ^T * B, что является матрицей 4x1. Его i-я запись дается выражением x ₁ⁱ * y ₁ + x ₂ⁱ * y ₂ +... + x _Nⁱ * y _N.
Умножьте нашу созданную вручную матрицу 4x4 (A ^T * A) ^-1 нашей созданной вручную матрицей 4x1 A ^T * B, чтобы получить матрицу 4x1 коэффициентов наименьших квадратов для вашей кубики.

Удачи!

Ответ 3

Вам не следует делать полную инверсию матрицы по причинам стабильности. Сделайте LU-разложение и обратную замену. В противном случае другие решения будут найдены.