Нахождение коэффициентов аппроксимирующего полинома решением нормальной системы уравнений C#
Задача: имеется набор данных х[х1, х2, ..., xn] и y[y1, y2, ..., yn] (если конкретно -- в моей задаче это зависимость потоков от длины волны, т.е. спектр звезды). Также имеется набор других данных х[х1, х2, ..., xn] и r[r1, r2,..., rn]. Необходимо нормировать набор данных y[y1, y2, ..., yn] на набор данных r[r1, r2,..., rn].
Реализация: для этого найдем коэффициенты интерполирующего полинома и поделим набор y[yn] на него.
Пусть m() - столбец коэффициентов. Его размерность равна степени интерполирующего полинома Q {upd:} плюс 1, т.к. у еще необходим коэффициент для свободного члена, итого Q+1;
G() - матрица Вандермонда размерности n x (Q+1), составленная из набора данных x[xn];
r() - матрица размерности n x n, на главной диагонали которой располагается набор r[rn], все остальные элементы матрицы равны нулю;
y() - столбец данных y[yn].
Введем обозначения: G*r = A, (A)^T - транспонированная матрица
Чтобы найти коэффициенты полинома, необходимо решить нормальную систему уравнений вида
(A)^T * A * m() = (A)^T * y()
Проблема: матрица (A)^T * A является переопределенной - число столбцов в ней (Q+1, число неизвестных коэффициентов полинома) много меньше числа строк (n). Чтобы решить такую систему, например, методом оптимального исключения, необходимо сделать матрицу квадратной.
Вопрос: как привести матрицу к квадратному виду? Можно ли это сделать? Или же можно решать систему для каждого набора Q+1 уравнений отдельно?
Если эти операции сделать нельзя, то как решить такую систему?
Ответы (1 шт):
Если вам нужно найти min[ Summ{i=1,n} {P(x) * r(x) - y(x)}^2], то
Если
rиyзаданы в разных точках, тоrможно интерполировать кусочно линейными или кубическими сплайнами. Еслиrварьируется в умеренных пределах, то даже линейной приближение вполне подойдёт для практических задач.При аппроксимации полиномом задача наименьших квадратов сводится к системе линейных уравнений вот каким образом.
Полином P(x) степени m - это линейная форма от (m+1) коэффициента. Следовательно, Summ{i=1,n} {P(x) * r(x) - y(x)}^2 - это квадратичная форма от коэффициентов. В глобальном экстремуме частные производные по коэффициентам P должны быть равны нулю, но эти частные прозводные являются линейными формами. То есть приравняв нулю m+1 частную производную вы получите систему из m+1 линейного уравнения.
Интерполяция по n точкам - это полином степени n-1, то есть Q == n-1 и ваши матрицы G и A - квадратные. Ну а с квадратной матрицей всё понятно, метод Гаусса наше всё :)
Если в вашей задаче предполагается, что Q < n, то это не задача интерполяции, а задача аппроксимации, и она решается совсем иначе (например, методом наименьших квадратов).
Из общих соображений, использовать интерполяцию - плохая идея. Во-первых, задача интерполяции, как правило, вычислительно неустойчива. То есть при небольших изменениях входных данных интерполяционный полином может радикально измениться. Во-вторых, при интерполяции полиномами даже сравнительно небольших степеней могут возникать непредсказуемые по размеру осцилляции, многократно превышающие вариацию в исходных данных (пример - явление Рунге)
Не случайно в реальной жизни используют кусочную интерполяцию максимум кубическими сплайнами.