Подгонка данных к полиному 3-й степени

Question

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

Часть проблемы заключается в том, что я не могу использовать различные числовые пакеты, такие как 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?Я предполагаю, что я пытаюсь переписать вышеупомянутый алгоритм в более простой случай, который работает для подгонки к кубическому полиному.Я не ищу решение кода, но если кто-то может указать мне правильное направление, я был бы признателен.

Answer 1

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

Если вы не можете работать с матрицами 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) th запись (я начинаю с (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 для вашей кубики.

Удачи!

Answer 2

Да, мы можем ограничить проблему вычислениями с «матрицей 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-го столбец А, должно быть ясно, что:

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

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

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

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

, где k = 0, .., 6. Матрица 4x4 C имеет «полосатую» структуру, что означает, что отображаются только эти семь значений. Зацикливая список 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).

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

Answer 3

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