Главная Методы условной оптимизации [0] [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [ 29 ] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73] [74] [75] [76] [77] [78] [79] [80] [81] [82] [83] задачи выбора параметров теперь будет выглядеть так; t, найти m\n\(if(x, t)-S(t))dt, где f (t)- «целевая траектория». Заменив интеграл суммой с помощью какой-нибудь подходящей квадратурной формулы, мы снова приходим к задаче о наименьших квадратах; найти т!п5!(ф(х, t,)~7 (t,)), где через ф и обозначены произведения ф и на веса принятой квадратурной схемы. В данном случае оптимальное значение критерия будет мальш, если не требовать от объекта «невозможного». Хотя для минимизации функций (4.58) можно использовать универсальные методы, как правило, этого не делают, а обращаются к специальньш алгоритмам, разработаннььм именно для задач о наименьших квадратах. Последние достаточно специфичны, и естественно попытаться учесть это. В частности имеется в виду особая структура градиента функции (4.58) и ее матрицы Гессе. Обозначим через У(л:)mxn-матрицу Якоби для f{x), и пусть С((х) - матрица Гессе для /,(х). Тогда выражения для градиента g{x) и матршщ Гессе G(x) функции (4.58) будут выглядеть так: g(x) = J(xYnx)\ (4.59а) G{x) = J(xYJ(x) + Q{x). (4.59b) Здесь через Q{x) обозначена сумма 11!1(х)С,{х). Из (4.59Ь) видно, что элементы G(x) образованы характерными комбинациями первых и вторых производных функций fj. Большинство алгоритмов решения задач о наименьших квадратах опирается на предположение о том, что слагаемое У (xY J (х) в (4.59Ь) рано или поздно станет доминирующим. Это предположение не соблюдается, если минимальные невязки велики, т. е., грубо говоря, если норма (/(х*)! сравнима с максимальным собственным значением матрицы J {xTY J {х). В данном случае преимущества специальных методов над универсальными теряются. Однако чаще встречаются ситуации, когда не вязки в решении достаточно малы и применение специальной техники оправданно. 4.7.2. МЕТОД ГАУССА - НЬЮТОНА Обозначим через х текущую оценку решения задачи минимизации функции (4.58) и условимся, что нижний индекс k у любой величины будет означать ее принадлежность Хд. Тогда ньютонов- ская система (4.17) в силу (4.59) примет вид (4.60) Ее решением будет некоторый вектор рд. (ньютоновское направление). Если прн стремлении Од к оптимуму норма Ц/дЦ приближается к нулю, то матрица Q,, тоже будет приближаться к нулевой. Значит, ньютоновское направление можно атроксимироеагпь решением системы JlhPH - Jlh- (4-61) Заметим, что эта система заведомо разрешима н включает только первые производные от f. Решение (4.61) представляет собой оптимальный вектор в линейной задаче о наименьших квадратах вида (4.62) найти min-i-yjp + ft. Если матрица имеет полный столбцовый ранг, этот вектор определен однозначно. Мы будем обозначать его через рд.. Метод, в котором Род, используется в качестве направления поиска, известен под названием метод Гаусса-Ньютона, а сам вектор Pojv принято называть направлением Гаусса-Ньютона. Когда норма IQ (х близка к ну.лю и матрица J. имеет полный столбцовый ранг, направление Гаусса -Ньютона мало отличается от ньютоновского. Точнее говори, если при достаточно малом положительном е выполнено равенство Q (х = е, то тем самым обеспечено соотношение Ру° =0(£). Следовательно, если норма f (х) равна нулю и столбцы матрицы / (х*) линейно независимы, метод Гаусса-Ньютона может достигать квадратичной скорости сходимости, хотя при расчете PoN учитываются только первые производные. В ранних реализациях метода Гаусса-Ньютона для решения системы (4.61), как правило, использовали какую-нибудь универсальную процедуру, предварительно формируя матрицу JIJ,, явным образом. Этот подход обладает одним серьезным недостатком: поскольку число обусловленности произведения JlJ равно квадрату числа обусловленности матрицы /д, оно может оказаться очень большим даже при неплохо обусловленной Уд, а это означает, что расчет направления поиска будет сопровождаться большими ошибками, которых можно было бы избежать. Например, работая на машине с двенадцатиразрядной точностью при числе обусловленности матрицы Уд, равном 10«, мы вправе рассчитывать на шестиразрядную точность вычисления рдд,, однако при решении системы (4.61) обычным способом в рассматриваемой ситуации можно получить ответ, не содержащий ни одной правильной цифры. Угроза неоправданной потери точности при «лобовой» реализации метода Гаусса - Ньютона вполне реальна, так как для прикладных нелинейных задач о наименьших квадратах плохая обусловленность является скорее правилом, чем исключением. Это относится прежде всего к задачам, возникающим из потребностей оценивания параметров математических моделей, и связано с тем, что такие модели часто бывают некорректно определенными. Чтобы не утяжелять последствий плохой обусловленности, вектор PoN надо искать как решение задачи (4.62) с использованием полной ортогональной факторизации (см. разд. 2.2.5.3) или разложения по особым числам (см. разд. 2.2.5.5). Особенно тщательно метод Гаусса - Ньютона должен быть реализован в случаях, когда нет гарантии линейной независимости столбцов Jft. Сразу отметим, что если /ь имеет дефект ранга, то у задачи (4.62) будет целое многообразие решений. Следовательно, требуется правило выбора конкретного вектора из этого многообразия. Лучше всего брать в качестве р решение задачи (4.62) с минимальной евклидовой нормой. Используя полную ортогональную факторизацию или разложение по особым числам, построить его нетрудно. Каждая реализация метода Гаусса - Ньютона, в которой при линейной зависимости столбцов направление Pow определяется как решение (4.62) минимальной нормы, долиша включать некую стратегию оценивания ранга J. Это означает, что в процедуре факторизации, используемой для расчета рдг, надо будет предусмотреть какие-то дополнительные действия, которых нет в стандартных процедурах, предполагающих, что ранг г факторизуемой матриць! известен заранее и ее первые г столбцов линейно неззвисимьт. К сожалению, определение понятия «ранг» в контексте приближенной арифметики с плавающей запятой теряет универсальность и становится проблемно-зависимым. Вопрос о ранге матрицы в конкретном случае остается открытым до тех пор, пока не введено явное масштабирование, т. е. не установлено, какие величины следует считать «пренебрежимо малыми». Пример 4.10. Суть дела легко понять на примере матрицы где е - ненулевое число, малое по модулю относительно единицы. С формальной точки зрения столбцы матрицы J линейно независимы, и ее ранг равен двум. Однако в действительности второй столбец может быть машинной версией первого, и тогда их следует признать «численно эквивалентными»; в этом случае <фанг» должен считаться равньш единице. Значит, вывод о значении ранга матрицы / будет определяться тем, следует ли пренебречь числом е в задаче, откуда она возникла. Если бы не погрешности машинной арифметики, выявить линейную зависимость столбцов матрицы во время построения ее полного ортогонапьного разложения было бы просто; когда факторизация осуществляется последовательными преобразованиями Хаусхолдера (см. разд. 2.2.5.3), зависимость очередного, (А-Ы)-го столбца от k предыдущих проявилась бы в равенстве нулю всех его преобразованных компонент с номерами от k-\-l до т. Таким образом, при точных вычислениях необходимым условием точной линейной зависимости является обращение в нуль некоторых величин. Естественно возникает вопрос; если вычтюлення осуществляются приближенно, го не будет ли «близосты> этих величин к нулю необходимым условием «приближенной» линейной зависимости? Оказывается, не будет, и это легко доказать от противного. Иначе была бы справедлива теорема о тем, что плохо обусловленная треугольная матрица должна иметь хотя бы один малый диагональный элемент, а известно, что можно построить треугольную матрицу, обусловленную сколь угодно плохо и не имеющую «малых» элементов на диагонали. Итак, исходя при оценивании ранга матрицы во время ее (ЭЛ-факторизации из «малости» элементов необнуленпого субдиагонального блока, мы рискуем дать завышенную оценку. Несколько проще оценивать ранг иа основе разложения по особым числам: близость вырождения обязательно проявится в малости этих чисел. Однако и тут могут во:тннкнуть определенные трудности. Возьмем, к примеру, разложение с особыми числами 1, 10-, 10-..... 10-»°, 10-"..... 10-"" Для него выбор порога «малости» существенно повлияет на оценку ранга, и, поскольку этот выбор в значигельнон степени произволен, здесь легко ошибиться. Проблема построения оценок ранга для метода Гаусса - Ньютона важна потому, что от них очень сильно зависит его функционирование. Это видно из следующего примера. Пример 4.11. Пусть / - двумерная функция, компоненты которой fl и fl являются в некоторой точке величинами порядка еди-ннць1, а ее матрица Якоби в этой точке равна где е - число, малое относительно единицы. Тогда, считая, что ранг J равен двум, в качестве направления Гаусса - Ньютона мы получим вектор Если же пренебречь величиной е и принять единичную оценку ранга J, направлением поиска будет (4.64) Угол между векторами (4.63) и (4.64) близок к прямому, причем первый из них почти ортогонален градиенту. Таким образом, изменение оценки ранга на единицу повлекло за собой разворот направления поиска почти на девяносто градусов. К сожалению, приходится признать, что такой стратегии оценивания ранга в методе Гаусса - Ньютона, которую можно было бы порекомендовать на все случаи жизни, не существует. Например, когда F похожа на плохо обусловленную квадратичную форму, лучше всего оценивать ранг по максимуму Это следует из соображений близости к методу Ньютона. Однако в других случаях данная стратегия может оказаться неподходящей, так как есть доводы и в пользу заниженных оценок. В частности, при почти вырожденной матрице Ju завышенная оценка ранга может привести к очень большому по норме решению (4.62), в то время как меньшая оценка часто дает вектор более разумной длины, практически не изменив минимальной невязки в (4.62). Кроме того, анализ возмущений для задачи (4.62) показывает, что относительное изменение ее точного решения пропорционально относительному изменению /д с коэффициентом (cond (Уд)) 3. Соответственно, занижая оценку ранга и уменьшая тем самым величину cond (Уд), мы можем избежать необходимости решения при расчетеплохо обусловленной задачи. 4.7.3. МЕТОД ЛЕВЕНБЕРГА - МАРКАРДТА Распространенной альтернативой методу Гаусса - Ньютона является метод Левенберга - Маркардта. В нем направление поиска определяется как решение системы уравнений вида (llJh+Kl)Pk=-Jlh, (4.65) где Хд-некоторое неотрицательное число. В этом методе шаг вдоль рд всегда полагается единичным, т. е. очередной точкой 4+1 будет д;д-Ь Рд. Можно показать, что рд представляет собой решение задачи на условный минимум вида найти тгп-1Удр-[-/д при ограничении lps<A, где Д-параметр, связанный с \. Таким образом, метод Левенберга-Маркардта относится к классу методов доверительной окрестности, обсуждавшихся в замечаниях к разд. 4.4. Монотонное убывание минимизируемой функции достигается в нем за счет подбора «хороших» значений Ад. При Х», равном нулю, рд будет направлением Гаусса-Ньютона, когда Хд стремится в бесконечность, норма рд стремится к нулю и вектор Рд в пределе становится параллельным антиградиенту. Следовательно, неравенство (Хд-Ьрд)<д всегда можно обеспечить, выбрав Хд достаточно большим. Обозначим через рьм (К) решение системы (4.65) при каких-то Хд и положительном Хд. Оказывается, что если матрица Уд имеет дефект ранга, то в общбм случае независимо от величин \\Q (Хд)11 и Хд будет справедливо соотношение 1£л1£оиМ1=0(1). Ip/vII *4.7.4. КВАЗИНЬЮТОНОВСКИЕ МЕТОДЫ Методы Гаусса-Ньютона и Левенберга-Маркардта опираются на предположение о доминирующей роли слагаемого УдУд в правой части (4.Б9Ь); исходя из него, матрицей Сд предлагается пренебречь. Это предположение ие вполне оправданно для так называемых задач с большой невязкой, в которых норму 7(х*) нельзя считать «малой», но которые еще не окончательно утрачивают специфики задач о наименьших квадратах. Если нормы 0;(л:) для 1=1, т являются величинами порядка единицы, то Q (jf) будет существенной составляющей матрицы Гессе функции (4.58), когда [/(х*)! превышает линилшльное собственное значение У [xY J (х*). При этом максимальное собственное число последней может существенно превосходить 1/()11. и именно такие задачи мы назовем задачами с большой невязкой. Для них и предназначены описанные ниже специальные методы. Что же касается задач с очень большими значениями 7(.«*)li. то никакой выгоды нз особой структуры минимизируемых функций для них извлечь не удается, и здесь следует обращаться к универсальным методам. Отказавшись от вычисления вторых производных, метод решения задач с большой невязкой можно строить иа основе квазиньютоновских приближений Мд неизвестной составляющей Q (х) матрицы Гессе. Надо отметить, что использование точного значения другой составляющей несколько усложняет формулы по сравнению со случаем, когда квазиньюгоновская аппроксимация применяется для оценки матрицы Гессе в целом. Система уравнений для расчета направления поиска при рассматриваемой организации вычислений выглядит так: (УгУд-ЬМд)рд = -Уд7д, а аналогичным выражению (4.30) условием, которому должно подчиняться очередное квазиньютоновское приближение Л1д4.„ [0] [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [ 29 ] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73] [74] [75] [76] [77] [78] [79] [80] [81] [82] [83] 0.0011 |