Подгонка полинома с использованием L1-нормы
У меня есть n точек (x0, y0), (x1, y1) … (xn, yn). n мало (10-20). Я хочу подогнать эти точки полиномом низкого порядка (3-4): P (x) = a0 + a1 * x + a2 * x ^ 2 + a3 * x ^ 3.
Я выполнил это, используя метод наименьших квадратов в качестве метрики ошибки, то есть минимизировать f = (p0-y0) ^ 2 + (p1-y1) ^ 2 + … + (pn-yn) ^ 2. Мое решение использует разложение по сингулярным числам (SVD).
Теперь я хочу использовать норму L1 (абсолютное значение расстояния) в качестве метрики ошибки, то есть минимизировать f = | p0-y0 | + | p1-y1 | + … + | pn-yn |.
Существуют ли какие-либо библиотеки (желательно с открытым исходным кодом), которые могут это делать и которые можно вызывать из C ++? Есть ли доступный исходный код, который можно быстро изменить в соответствии с моими потребностями?
Решение
Регрессия L_1 на самом деле довольно просто формулируется как линейная программа. Вы хотите
minimize error
subject to x_1^4 * a_4 + x_1^3 * a_3 + x_1^2 * a_2 + x_1 * a_1 + a_0 + e_1 >= y_1
x_1^4 * a_4 + x_1^3 * a_3 + x_1^2 * a_2 + x_1 * a_1 + a_0 - e_1 <= y_1
.
.
.
x_n^4 * a_4 + x_n^3 * a_3 + x_n^2 * a_2 + x_n * a_1 + a_0 + e_n >= y_n
x_n^4 * a_4 + x_n^3 * a_3 + x_n^2 * a_2 + x_n * a_1 + a_0 - e_n <= y_n
error - e_1 - e_2 - ... - e_n = 0.
Ваши переменные a_0, a_1, a_2, a_3, a_4, errorи все e переменные. x а также y данные вашей проблемы, так что это не проблема, что x Появляется второй, третьей и четвертой степени.
Вы можете решить задачи линейного программирования с GLPK (GPL) или lp_solve (LGPL) или любое количество коммерческих пакетов. Мне нравится GLPK, и я рекомендую использовать его, если его лицензия не является проблемой.
Другие решения
Да, это должно быть выполнимо. Стандартный способ постановки задач подбора полиномов в виде множественной линейной регрессии состоит в определении переменных x1, x2 и т. Д., Где xn определяется как x. ^ N (поэлементное возведение в степень в нотации Matlab). Затем вы можете объединить все эти векторы, включая пересечение, в матрицу проектирования X:
X = [1 x1 x2 x3]
Тогда ваша проблема подбора полиномов — это проблема регрессии:
argmin_a (| y — X * a |)
где | | нотация — это желаемая функция стоимости (для вашего случая, норма L1), а a — вектор весов (извините, SO, насколько я могу судить, не имеет хороших математических разметок). Подобные регрессии известны как «устойчивые регрессии», и у Numeric Recipes есть процедура для их вычисления: http://www.aip.de/groups/soe/local/numres/bookfpdf/f15-7.pdf
Надеюсь это поможет!
Проблема с нормой L1 состоит в том, что она не дифференцируема, поэтому любые минимизаторы, которые зависят от производных, могут потерпеть неудачу. Когда я пытался свести к минимуму такие функции, используя, например, Минимизация сопряженного градиента, я обнаружил, что ответ застревает на излом, то есть х = 0 в функции у = | х |.
Я часто решаю эти математические задачи из первых принципов. Одна идея, которая может сработать, заключается в том, что целевая функция будет кусочно-линейной по коэффициентам вашего полинома низкого порядка. Так что это можно решить, начав с полинома, который получается из наименьших квадратов, и затем улучшив решение, решив ряд линейных задач, но каждый раз лишь переходя от вашего текущего лучшего решения к ближайшему кинку.