Здравствуйте.
У меня научная задача по байесовой подгонке данных под модель.
Читала книги, статьи, написала программу, не пашет. Консультировалась с единственным человеком который этим занимался, оказалось, что минимизацией функции с плохо обусловленным гессианом - не занимался. Никто на кафедре и близлежащих не смог помочь. Я реализовала программно алгоритм, пытаюсь инкорпорировать алгоритмы из свободного доступа (но они во-первых часто написаны на старом фортране, который ужасно разбирать, во-вторых, часто не в том виде, который нужен - потому тупо пишу алгоритм из книги и корректирую, если разбираюсь что они сделали), проверяла и так и этак, нарисовала срезы функции графически, зависимость функции по направлению итерации (и удостоверилась, что ньютоновский шаг хотя бы правильный и сама функция в этом направлении убывает), проверяла степень вычислительной стабильности обращения матрицы (см.ниже), но что-то явно упускаю из виду.
Итоговая проблема сводится к следующему:
Есть функция наименьших квадратов (плохо обусловленная, т.е. гессиан - симметричный, положительно определенный, но с разбросом собственных значений на порядки) плюс регуляризующая параболическая функция (просто
,
- вектор начального приближения, A - диагональная матрица, пока пропорциональна единичной, т.е. просто число)
Нужно искать минимум получившейся функции для ряда значений A и искать оптимальное A. Условие оптимальности в вопрос к форуму не входит. В вопрос входит - как получить минимум плохо обусловленной функции при малом (минимально достаточном) регуляризаторе.
Авторы работы "Bayes-Turchin approach to x-ray absorption fine structure data analysis, H. J. Krappe and H. H. Rossner,PHYSICAL REVIEW B 66, 184303 (2002)" и еще ряда работ на ту же тему 1999 и 2004 годов, с которой я, собственно, списывала алгоритм, ухитрились получить адекватный минимум при A=1-100 и максимальном собственном значении гессиана 1e5.
Для этого они использовали метод Ньютона с линейным поиском длины шага (damped Newton, не знаю, корректно ли перевожу). Оный метод сходится у них все медленнее при уменьшении регуляризатора, но до A=1 они добили как-то. Метод есть в ряде разных форм, но сходимость доказана в случае выбора одной из форм "аффинной инвариантности", т.е. нечувствительности метода к замене переменных либо в пространстве переменных, либо в пространстве градиента функции F(x) (правой части системы нелинейных уравнений, эквивалентных минимизации).
Вариант, который они использовали, следующий:
берем "функцию уровня" - квадратную норму
упрощенного ньютоновского шага, но с гессианом в точке
:
где x берется вдоль ньютоновского шага
и начиная с длины шага 1 (т.е. равной ньютоновскому шагу) уменьшаем оную длину p, пока не окажется, что в новой точке оная функция уровня меньше, чем в начале шага - монотонное убывание. Ну там еще множитель, чтобы убывание было строгим.
Суть в том, что этот метод должен сходиться и для плохо обусловленных матриц, и доказано, что направление шага Ньютона является для этой функции просто направлением наискорейшего спуска (и независимо от обусловленности):
Здесь обратный гессиан берется в точке с которой стартует итерация
, т.е. он постоянен по x и по нему не дифференцируем.
Доказательство из книги Deuflhard-а, на которую авторы той статьи ссылаются. Newton methods for nonlinear problems
Вроде все банально и работать должно?
Ан нет, я получаю, что при хорошей обусловленности (A на порядок меньше, чем максимальное с.з. гессиана) норма шага и правда убывает, а как только обусловленность достигает двух-трех порядков, норма шага внезапно растет в этом направлении в прямое нарушение теоремы.(выделила вопрос, чтобы не потерялся)
Хотя вроде бы доказано, что от обусловленности зависит, если мы функцию уровня возьмем с единицей вместо обратной матрицы (второй подход), а вот здесь же специально сделано так что именно в точке
производная по направлению будет отрицательной! Как обусловленность влияет на прямое перемножение обратной матрицы с ней самой?
Я думала, что надо что-то перемасштабировать во избежание ошибок округления и соответственно чувствительности задачи к малым ошибкам, раз она плохо обусловлена. Но простой тест, когда я добавляю к гессиану и градиенту к каждой компоненте (поочередно все компоненты проверяю) малый сдвиг (1e-8, при двойной точности), показал, что вычисленный шаг и его норма не меняются, т.е. все стабильно. Хотя поскольку я физик-расчетчик, а не программист и не математик-теоретик, возможно, надо более подробно этот вопрос проработать, но я не знаю пока, как это лучше сделать.
Плохую обусловленность саму по себе убрать нельзя - в том и вся загвоздка, это фишка метода, что он позволяет минимально регуляризовать задачу, так, чтобы сделать ее сходящейся, но не повлиять на те компоненты решения, которые хорошо обусловлены. Т.е. минимальный\оптимальный регуляризатор по Тихонову.
Может ли кто-нибудь помочь понять, как проработать влияние ошибок, и как проверить вот это доказательство, которое я тут привела, впрямую? Занимается ли тут кто-нибудь конкретно этим методом?
Код на фортране и очень длинный, если включать саму модель. Могу скинуть кусок кода, связанный собственно с итерацией и поиском шага, и полученные матрицы-вектора, если надо. Мне кажется, лучше в алгоритмах говорить.
Если кто найдется, буду очень благодарна, потому что хода назад нет, задачу мне надо решить, а я уже пару лет над ней бьюсь, так как в институте чуть ли никто этим не занимается в таком виде. Хоть меняй направленность, но уже заявлено все... А оптимизацию нам не читали, приходится самой книги читать, равно как и численные методы на компьютере были, но не углубленно. Я могу проработать простое доказательство в матрицах, могу написать программу, но "плаваю" в операторном анализе безнадежно (ну, могу освоить курс, но времени мало, надо скипнуть то, что не нужно на данном этапе. Но вроде нужно только вышенаписанное доказательство, а с ним же все просто), да и с контролем точности в программе проблема - я знаю только, что все желательно отнормировать и что надо как-то проверять чувствительность к начальным данным, и что, как и в задаче в целом, в мелочах неточность зависит от степени обусловленности матриц как раз. Но не знаю, как это дотошно проверить.
т.е. либо как лучше проверить оное доказательство численно, либо как проверить стабильность по отношению к малым ошибкам в данном случае, если простая проверка ничего не дала, либо, может, этот метод кто-то знает и скажет что он не сходится в этом случае или его надо модифицировать.
Есть еще вариант - взять метод другой совсем, найти с его помощью минимум функции, а на последнем этапе включить Ньютон (локально-то можно обойтись без выбора длины шага), т.к. итоговый результат в окрестности решения все равно требует вычисления гессиана и т.п. Вопрос: какой из методов оптимизации не убьется о плохо обусловленную матрицу и не уйдет в таких условиях на край области определения? Симплекс в оригинальном программном обеспечении, которым мы раньше пользовались, радостно убивается. Собственно потому я за это и вынуждена была взяться. Но я просто не знаю всего спектра методов, то, что копала, не подходит как раз по сходимости на плохо обусловленных матрицах. Но ведь если в статье четко метод Ньютона с выбором шага, и там написано какая обусловленность, то и у меня должен работать, зачем далеко ходить.
Проверяю метод на тестовом случае, т.е. там стопроцентно есть минимум и он в центре области определения. Кажется, если гессиан везде положительно определенный, то там вообще не должно скатываться вбок, т.е. ситуация, когда на краю области определения ньютонов шаг указывает за ее пределы - невозможна, там же один минимум. Или возможна? Это я тоже не могу понять как доказать. Что минимум один. Но скатывается он из-за плохого выбора шага, я так понимаю. Или есть еще проблема?
маленькое, то мы имеем плохую обусловленность промежуточных задач, но зато большую длину шага и лучшую скорость сходимости. И, соответственно, наооборот. Всё нормально. Так и должно быть.