Не могу ничего толкового найти в инетах и книжках.
Дело в следующем:
Есть некие реальные измерения, представленные набором данных 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 , графички там в конце похожи на то, что мне нужно.
Ну или хотя бы подскажите что надо искать.