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