LINUX.ORG.RU

аппроксимация полиномом 3-й степени, добавить монотонность

 , ,


0

4

Не могу ничего толкового найти в инетах и книжках.

Дело в следующем:

Есть некие реальные измерения, представленные набором данных x0,x1,x2...xN и y0,y1,y2...yN. Ближе к x=0 зависимость почти линейная, потом y начинает насыщаться.

Это дело аппроксимируется полиномом 3-й степени (y = a0 + a1*x + a2*x*x + a3*x*x*x) методом наименьших квадратов посредством QR-разложения (LSE polynomial fit using QR decomposition). В результате находятся коэффициэнты полинома a1, a2, a3. Т.к. полином должен проходить через начало координат, то a0 = 0.

Всё отлично работает и в общем-то даже понятно как и почему это работает.

Но есть проблема - для данной задачи полином всегда должен быть монотонно возрастающим (без перегибов) на участке от 0 до xN. Вот тут-то внезапно обнаружились грабли. Чтение литературы разъяснило, что в LSE-решалку надо всего-навсего добавить constraint в котором указать что производная полинома для используемого диапазона x должна быть больше нуля. Но, сцуко, я нигде не нашёл _как_конкретно_ это сделать, в частности, для решалки LSE через QR (да пофиг, хоть SVD, хоть LU).

В принципе, я догадываюсь, что в матрицу, которая, собственно, разлагается этим самым QR-разложением, надо добавить дополнительных строчек и в вектор y добавить соответствующих им значений. Но вот что это должны быть за строчки - я в понятном виде нигде не нашёл. А неравенство (например, «для каждого x из набора данных производная полинома должна быть > 0») я не представляю как и куда можно добавить в эту straight-forward решалку.

Есть кто тут сведущий в матане, способный ткнуть меня в литературу, где описание решения этой задачи не заканчивается словами «just constraint solver to f'(x) >= 0», а есть описание того как именно это cделать. И без ссылок на всякие сраные Матлабы (use MATLAB CONSTR function with constraints defined earlier) и пр. Типичный пример - http://ws680.nist.gov/publication/get_pdf.cfm?pub_id=17206 , графички там в конце похожи на то, что мне нужно.

Ну или хотя бы подскажите что надо искать.

★★★★★

Может тебе совсем не полином нужен? Попробуй логистической функцией например.

dikiy ★★☆☆☆ ()

just constraint solver to f'(x) >= 0

ну а что тебе еще надо-то? Добавь условия соответствующие. Мне кажется, что как раз ман к нужной функции матлаба это то, что доктор прописал.

dikiy ★★☆☆☆ ()
Ответ на: комментарий от dikiy

Может тебе совсем не полином нужен? Попробуй логистической функцией например.

С превеликим удовольствием. S-like curve изумительно подойдёт. Более того, с точки зрения физики, там вполне может быть что-то типа (1-exp(-x)/1+exp(-x))

Вот только никак не могу найти S-like функцию у которой наклон линейной части и «острота» перехода к предельному значению могут регулироваться разными параметрами независимо.

Stanson ★★★★★ ()
Ответ на: комментарий от Stanson

Вот только никак не могу найти S-like функцию у которой наклон линейной части и «острота» перехода к предельному значению могут регулироваться разными параметрами независимо.

острота перехода, вероятно, регулируется коэффициентом у экспоненты. Параметр a в exp(-ax).

Но это по нюху. Он может меня и обманывать.

dikiy ★★☆☆☆ ()
Ответ на: комментарий от knovich

А задать условия на сами коэффициенты нельзя? Типа f(a1,a2,a3)==0

f(a1,a2,a3,x)==0 не поможет. Тебе неравенство надо.

но вообще это изврат какой-то. Возьми S-образную функцию какую-нить и подкрути.

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

Потом получившуюся функцию проинтегрируй со сглаживающим ядром, и получишь нужную сглаженность отсроты.

dikiy ★★☆☆☆ ()

по-моему

благородному Дону нужно искать по словами constrained optimization. Ну и еще мне кажется что LSE using QR decomposition - это как раз про то «как добавить» дополнительных строчек в вектор (сама декомпозиция это и делает). А вот дополнительное ограничение - это как раз поломало бы исходный алгоритм.

sshestov ()

меня IDL (это что-то типа матлаба в некотором виде)

ссылает на нижепреведенную литературу, а метод называет Generalized Reduced Gradient Method

1. Lasdon, L.S., Waren, A.D., Jain, A., and Ratner, M., “Design and Testing of a Generalized Reduced Gradient Code for Nonlinear Programming”, ACM Transactions on Mathematical Software, Vol. 4, No. 1, March 1978, pp. 34-50.

2. Lasdon, L.S. and Waren, A.D., “Generalized Reduced Gradient Software for Linearly and Nonlinearly Constrained Problems”, in “Design and Implementation of Optimization Software”, H. Greenberg, ed., Sijthoff and Noordhoff, pubs, 1979.

3. Abadie, J. and Carpentier, J. “Generalization of the Wolfe Reduced Gradient Method to the Case of Nonlinear Constraints”, in Optimization, R. Fletcher (ed.), Academic Press London; 1969, pp. 37-47.

4. Murtagh, B.A. and Saunders, M.A. “Large-scale Linearly Constrained Optimization”, Mathematical Programming, Vol. 14, No. 1, January 1978, pp. 41-72.

sshestov ()
Ответ на: комментарий от dikiy

ну а что тебе еще надо-то? Добавь условия соответствующие.

Куда добавить-то???

Вот псевдокод:

есть набор данных:
x[] = { x0, x2, x3, ... xN };
y[] = { y0, y2, y3, ... yN };

Заполняем конкретными величинами матрицу Vandermonde (единичного столбца нету, потому что a0 всегда 0, график должен через начало координат проходить):

A[][] = {
{ x1, x1*x1, x1*x1*x1 },
{ x2, x2*x2, x2*x2*x2 },
{ x3, x3*x3, x3*x3*x3 },
...
{ xN, xN*xN, xN*xN*xN }
};

Натравливаем на A Householder'а:

QR_Householder( A, Q, R );

Транспонируем Q

Transpose( Q, QT );

Перемножаем QT и y, получаем вектор R * a

Ra = QT*y;

Обратной подстановкой находим вектор a из Ra и треугольной R - R*a = Ra

Получаем готовый к использованию полиноминал вида

f(x) = a[0]*x + a[1]*x*x + a[2]*x*x*x;

Куда здесь это неравенство надо засунуть?

Stanson ★★★★★ ()
Ответ на: комментарий от Stanson

Куда здесь это неравенство надо засунуть?

туда ты его не засунешь есессно. Тебе надо составить функцию Лагранжа для начала. Но я еще раз настоятельно предлагаю сменить подход :)

dikiy ★★☆☆☆ ()
Последнее исправление: dikiy (всего исправлений: 1)
Ответ на: комментарий от dikiy

острота перехода, вероятно, регулируется коэффициентом у экспоненты. Параметр a в exp(-ax).

Не, это четверть наклона, если про сигмоид речь. Наклон линейной части, в общем.

Stanson ★★★★★ ()

Ближе к x=0 зависимость почти линейная, потом y начинает насыщаться.
Но есть проблема - для данной задачи полином всегда должен быть монотонно возрастающим (без перегибов) на участке от 0 до xN.
just constraint solver to f'(x) >= 0

Не забудь к f'(x) >= 0 добавить ещё f"(x) <= 0. Это чтобы перегибов не было

Crocodoom ★★★★ ()
Ответ на: комментарий от Stanson

короче, берешь кусок линии, берешь две горизонтали. Определяешь кусочно гладкую функцию f из этих вот кусков. Берешь потом сглаживающее ядро phi(x,t)=Aexp(-B/(1-(x-t)^2)). Определяешь функцию

F(x)=int f(t)phi(x,t)dt.

F — будет той функцией, что ты ищешь.

Подробнее можно прочитать например у Смирнова «Курс высшей математики», том V, стр.218 и далее.

Интеграл возьмешь ручками. Коэффициенты A и B в принципе то, что тебя потом интересует при аппроксимации получившейся функцией F.

dikiy ★★☆☆☆ ()
Последнее исправление: dikiy (всего исправлений: 1)
Ответ на: комментарий от dikiy

туда ты его не засунешь есессно.

Не вопрос.

Тебе надо составить функцию Лагранжа для начала.

Вот в этом и заключается глобальный конфликт между теоретиками и практиками. :)

Можно как-то более приземлённо про составление из набора точек и конкретного требования про монотонность, и главное - что потом с этой функцией делать, чтобы в итоге получить искомые a1, a2 и a3?

Stanson ★★★★★ ()
Ответ на: комментарий от Stanson

Можно как-то более приземлённо про составление из набора точек и конкретного требования про монотонность, и главное - что потом с этой функцией делать, чтобы в итоге получить искомые a1, a2 и a3?

В треде тебе не ответить на этот вопрос. Твои варианты:

1)Идешь и читаешь книжку по оптимизации.

2)Идешь и делаешь в матлабе

3)Внимаешь советам про S-образную функцию

dikiy ★★☆☆☆ ()
Ответ на: комментарий от Crocodoom

Не забудь к f'(x) >= 0 добавить ещё f"(x) <= 0. Это чтобы перегибов не было

Формально да, но в данном конкретном случае мне более чем достаточно чтобы для каждого x из набора конкретнейших x0...xN f'(x) >= 0. Так что не в этом проблема.

Мне не понятно, каким образом и куда это условие запихивать в имеющуюся решалку.

В одной статейке попалось итерационное решение типа если прямое решение дало немонотонный на нужном участке полином, то добавить точку за пределами участка в которой производная форсирована в нуль, и типа продолжать этим заниматься, пока не получится полином монотонный, но что-то как-то это всё на костыль похоже.

Stanson ★★★★★ ()
Ответ на: комментарий от Zubok

Я кучу подобного уже прочитал, но так и не понял, как выведенные _неравенства_ (1) натянуть на прямое вычисление коэффициэнтов полинома.

Во всех подобных статьях всё начинается с «we try to fit polynomial f(x) = a0 + a1*x + a2*x^2 + a3*x^3 to given points (x,y)n», что мне понятно и я знаю как это сделать, потом каким-либо образом дело доходит до неравенств для производной, а потом, без объявления войны и без объяснения как именно эти неравенства прикрутить в решению изначальной, известной мне задачи, выдаются красивые монотонные графички.

В данной статье всё то же самое. Как конкретно получить из заданных точек синенькую кривую - я прекрасно знаю. А как конкретно они получили красную - мне совершенно не понятно. Что нужно изменить в процессе получения синенькой кривой, чтобы получить красненькую?

Я как бы не математик, и то, что у них там само собой разумеется, для меня вовсе не очевидно.

Stanson ★★★★★ ()
Ответ на: комментарий от dikiy

Идешь и читаешь книжку по оптимизации.

Какую именно? Желательно не для теоретиков, а для практиков. Что-то типа «численных методов».

Идешь и делаешь в матлабе

Мне это надо делать сотни тысяч раз, причём, желательно, как можно быстрее. Матлаб вообще ни разу не подойдёт.

Внимаешь советам про S-образную функцию

Вот над этим я подумаю. Но опять же, как рассчитывать численно коэффициэнты S-образной функции, особенно в случаях, когда данные есть только на промежутке, когда до насыщения, возможно, ещё очень далеко.

Stanson ★★★★★ ()
Ответ на: комментарий от Stanson

Ну, сначала ты спрашивал, как задать условие, что производная должна быть >= нуля. С этойточки сдвинулись. Приходим к ситуации, когда надо решить ту же оптимизационную задачу, но с тремя линейными ограничениями. Теперь надо смотреть Quadratic programming algorithms.

https://en.wikipedia.org/wiki/Quadratic_programming

Они как раз и решают эту конкретную задачу. В математических пакетах должны быть функции.

Zubok ★★★★★ ()
Последнее исправление: Zubok (всего исправлений: 1)
Ответ на: комментарий от Zubok

Они как раз и решают эту конкретную задачу. В математических пакетах должны быть функции.

С вики как раз в конце есть список кое-какой.

Zubok ★★★★★ ()
Ответ на: комментарий от Stanson

Но опять же, как рассчитывать численно коэффициэнты S-образной функции, особенно в случаях, когда данные есть только на промежутке, когда до насыщения, возможно, ещё очень далеко.

Примерно так же, как и для полиномов. У тебя всего два коэффициента будет, я думаю.

Какую именно? Желательно не для теоретиков, а для практиков. Что-то типа «численных методов».

без теории далеко не уедешь.

dikiy ★★☆☆☆ ()
Ответ на: комментарий от Zubok

Ну, сначала ты спрашивал, как задать условие, что производная должна быть >= нуля.

Ну уж производную полиноминала взять - не rocket science, это и так понятно.

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

Ну вот и спрашиваю, куда эти неравенства в готовую решалку засовывать.

Теперь надо смотреть Quadratic programming algorithms.

Там неравенство присутствует только в разделе Problem formulation. А я ищу Solution method.

В математических пакетах должны быть функции.

Видимо, таки придётся в итоге сырцы читать, а не умные книжки и статьи. :)

Stanson ★★★★★ ()
Ответ на: комментарий от Stanson

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

Книжка: Jorge Nocedal Stephen J. Wright. Numerical Optimization.

Со стр. 463 и до конца раздела рассматриваются методы с ограничениями в виде равенств и неравенств. Есть описание алгоритма (например, на стр. 472 для active-set method). Ну и там дальше есть - разделы 17, 18 и 19. Там те методы, которые также в википедии перечислены.

Но вдруг кто-то тебе поможет с более простым вариантом. Ты зато знаешь, что тебе надо искать. :)

Zubok ★★★★★ ()
Последнее исправление: Zubok (всего исправлений: 3)
Ответ на: комментарий от Zubok

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

Вот поэтому и не хочу в сырцы лезть.

Книжка: Jorge Nocedal Stephen J. Wright. Numerical Optimization.

Знатная книжка, хоть поподробнее чем обычно, спасибо.

Но вдруг кто-то тебе поможет с более простым вариантом.

Видимо простого нету.

Просто надо поискать.

Да обыскался уже. Английских математических терминов теперь больше чем русских знаю. Странно, но в рунете как-то тухло с чем-то мало-мальски выходящим за рамки обычных ВУЗовских курсов.

Stanson ★★★★★ ()

Мне кажется, ты не до конца понимаешь что надо и как делать. Мне честно ковырять это лень да и некогда, но можешь начать с забивания в гугл LSE QR. У меня одним из результатов является книжка (не весть текст) Linear Regression Analysis: Theory and Computing от Xin Yan. Может, если её почитать, станет немного легче, но я не читал её, так что прямо точно не подскажу.

peregrine ★★★★★ ()
Ответ на: комментарий от peregrine

но можешь начать с забивания в гугл LSE QR

C LSE QR давно всё в порядочке. Тут задача несколько шире.

Stanson ★★★★★ ()

Напиши функционал типа минимизация невязки по данным с условием на коэффициенты, которые означают монотонную производную. Потом делаешь метод мн. Лагранжа.

Но лучше взять и использовать монотонные сплайны того же третьего порядка. В Гугле ищется. Ключевое слово pchip, описание и алгоритм есть на википедии. Делается очень легко.

anonymous ()
Ответ на: комментарий от anonymous

Сплайны не очень годятся. Нужна максимально простая функция, да ещё и легко дифференцируемая (Jacobian и всё такое), ибо она и её дифференциал будут использоваться в аццких миллионах итерационных циклов, и любой просёр в производительности в ней приведёт к весьма ощутимым тормозам.

Начальный расчёт коэффициэнтов функции тоже не должен быть охочим до процессора, потому что при начальном расчёте из исходных данных нужно построить десятки тыщ разных кривых, так что весьма желателен быстрый и надёжный неитерационный способ нахождения решения типа однопроходного LSE QR.

Из чего-то относительно разжёванного я ещё нашёл NNLS (Non-Negative Least Square) который, в принципе, даст годную апроксимацию, но, сцуко, он опять же итерационный и прожорливый в смысле флопсов, и даже в виде FNNLS выглядит как-то печально. Частный случай Quadratic programming в общем-то.

Stanson ★★★★★ ()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.