Нелинейный метод наименьших квадратов



Скачать 177.81 Kb.
страница1/2
Дата02.06.2013
Размер177.81 Kb.
ТипДокументы
  1   2
Глава 4 Нелинейный метод наименьших квадратов
Постановка задачи. Что будет, если зависимость наблюдаемых значений yi от параметров нелинейная, т.е.

yi = h(xi ; ) + i = h i ( ) + i , i = 1, ..., n, (1)

где xi  значение (скалярного или векторного) аргумента в i-ом наблюдении, h  функция заданного вида, зависящая от параметров = ( 1 , ..., k)T, истинное значение неизвестно и подлежит оценке по yi; ошибки измерений i независимы и в совокупности для = (1 , ..., n )T выполняются условия E = 0, cov( , ) = 2In ? Свойства оптимальности оценки НК, доказанные в линейном случае, целиком основывались на линейности модели (5’) гл.1 по и на линейности НК-оценок

= Arg min || yh () || 2 , (2)

где h () = (h 1 ( ) , h 2 ( ) , ..., h n ( ) )T Rn . Ни то, ни другое не имеет места для h(x ; ) нелинейного вида. Оценка (2) обычно не может быть записана в явной форме, а получается численно в итерационной процедуре минимизации

Q ( ) = || yh () || 2. (3)

Поэтому в общем случае МНК не обоснован, и метод МП является единственным обоснованным методом :

МП = Arg max L ( ) , (4)

и в зависимости от ф.р. ошибок получим разные оценки (4) (см. пример 3 гл.1 для линейного случая).
В частности, если Nn (0 , 2 I n ), то

МП = Arg max [exp {1/(2 2 )  || yh () || 2 }] = Arg min Q ( ) = , (5)

т.е. получаем обоснование нелинейного МНК в случае нормальности ошибок.

Однако, в отличие от линейного случая, оценка (2) не будет нормально распределённой (вспомним, что результат теоремы 2.4 был основан на линейности = L y) . Поэтому вся теория нормальной регрессии тут не проходит. В частности, несмещённость и эффективность НК-оценки лишь асимптотическая, как и для любой МП-оценки.

Таким образом, рассмотрение нелинейного случая в основном касается вычислительной стороны  как искать (2) ? Однако прежде, чем искать, поставим вопрос о существовании и единственности НК-оценки.

Пример 1. Оценка (2) может не существовать. Пусть y1 = y2 = ... = yn = 0, R1 , h (x ; ) = , xi > 0, i = 1, ..., n . Тогда

Q ( ) = , inf Q ( ) = 0, но он не достигается ни в какой конечной точке , а лишь при    . 

Пример 2 . (решение (2) может быть не единственным). Пусть опять yi = 0, i = 1, ..., n , и пусть = (1 , 2 , 3 )T ,

h (x; ) = exp{ -1x } - exp{ -2x } - 3 (exp{ -x } - exp{ -10x }).

Нетрудно видеть, что опять inf Q ( ) = 0, но он достигается при

(а) 1 = 1, 2 = 10, 3 = 1; (б) 1 = 10, 2 = 1, 3 = -1; (в) 1 = a , 2 = -a, 3 = 0 ;   a    . 

Наконец, из (5) вытекает, что под НК-оценкой понимается , реализующая глобальный минимум Q( ) (3). Однако при нелинейной h(x; ) , как правило, Q( ) имеет несколько минимумов, а найденное * может давать не глобальный, а локальный минимум Q() (неподходящее решение). Задача глобальной минимизации очень трудна, и до настоящего времени в общем случае не решена. Поэтому, чтобы иметь хоть какую-то гарантию успеха, на практике пробуют несколько разных начальных приближений в одной из описанных ниже итерационных процедур, и убеждаются, что все эти начальные приближения приводят к одному и тому же решению (в противном случае за принимают то, которое даёт наименьшее значение Q() ).

Пример 3. Пусть n = 5, опять yi = 0; xi = i, i = 1, ..., 5 , и пусть R1 , h (x ; ) = 1 – cos(2 ) + 0.36 x 2. Ясно, что min Q ( ) = 0, и он достигается в точке = 0 (глобальный минимум). Но имеются ещё четыре локальных минимума вблизи точек -2, -1, 1, 2 (в этих точках функция f() = 1 – cos(2) принимает наименьшее значение, равное 0). Если использовать, как это обычно делается, итерационный алгоритм, и взять в качестве начального приближения для значение (0) , отделённое от нуля одним из локальных минимумов, то процедура локальной минимизации найдёт в качестве * локальный минимум.




Действительно, автор провёл вычислительный эксперимент, используя функцию FindMinimum из пакета Mathematica (version 2.2). Получились такие результаты:

(0) = 0.5 Q( (0) ) = 25.8455 * = 0.0251734 Q( *) = 8.70710 -8

(0) = 1.2 Q( (0) ) = 27.9141 * = 0.936574 Q( *) = 6.25756

(0) = 2.2 Q( (0) ) = 205.484 * = 1.86034 Q( *) = 99.5182.

Как же минимизировать целевую функцию (3) ? Для некоторых частных видов модельной функции h(x;) можно использовать подходящее преобразование переменных, делающее модельную функцию линейной. Например, если в примере 1 взять

g (x; ) = ln h (x; ) = - x , (*)

и если данные не искусственно сконструированные (y1 = y2 = ... = yn = 0), а более реальные (y1 , y2 , ... yn > 0) , можно рассмотреть zi = lnyi . Такова же ситуация в случае модельной функции – гиперболы

h(x;) = 1/ (1 + 2 x)

(g(x; ) = 1/ h(x;), zi = 1/ yi ). Однако при преобразованиях типа (*) теряется аддитивность ошибок:

zi = ln yi ln h(xi ; ) + i’ = g(xi ; ) + i’,

где i– ошибка для преобразованных данных. Поэтому минимизация целевой функции

(**)

    не то же самое, что минимизация (3), так что решение (**) не есть НК-оценка. Правда, если | i | << |h(x; )|, можно взять по формуле Тейлора линейную (по i ) часть разложения yi : zi = ln h(xi ; ) + ln (1+ )  ln h(xi ; ) + = g(xi ; ) + i’. Однако и в этом случае преобразованная ошибка i не имеет постоянной дисперсии, Var{i} = , так что нужно использовать оценку Эйткена (20) из гл.1, а не обычную НК-оценку (10) гл.1. Таким образом, преобразование переменных требует осторожности, и правильнее рассматривать решение линеаризованной задачи (**) в качестве метода выбора начального приближения для , а затем решать нелинейную задачу (3) итерационными методами.

Градиентные методы. Какими? Обычно модельная функция h(x;)  гладкая функция (в дальнейшем нам понадобится, чтобы h ( ), h ( ) и h ( ) были непрерывны), а Q() квадратичная по h (т.е. тоже гладкая) . Теория и практика оптимизации показывают, что для локальной минимизации гладкой функции (3) более эффективны градиентные методы, а не алгоритмы типа случайного поиска. При этом решение ищется итерационно :

(l+1) = (l) +  (l) , (6)

где l  номер шага, (l)  приближение для на l -ой итерации, (l+1)  на l+1-ой итерации,  (l)Rkитерационный шаг. Итерационный шаг определяется направлением l (единичный вектор, || l || = 1, лежащий на луче, проведённом из (l) ) , и величиной шага l > 0, так что

 (l) = l l . (7)
Хотелось бы, чтобы очередной шаг приближал нас к решению * (обозначаем его *, а не потому, что ищется не глобальный, а локальный минимум (3)). Однако мы не знаем * , поэтому считаем шаг удачным, если он уменьшил целевую функцию (3).

Определение 1. Шаг  (l) допустимый, если Q( (l) +  (l)) < Q( (l)) .

Определение 2. Градиентный метод называется допустимым , если он состоит из одних лишь допустимых шагов.

В дальнейшем мы будем рассматривать только допустимые градиентные методы. Общую итерационную процедуру такого метода можно описать так:

  1. Пусть начальное приближение для есть (1) . Положим l := 1.

  2. Находим l направление допустимого шага.

  3. Определяем величину шага l , вычисляем по (7)  (l) и (l+1) = (l) +  (l)

  4. Проверяем условия останова. Если они выполнены, полагаем *:= (l+1) и идём на шаг 5. Иначе присваиваем (l):= (l+1), l := l+1 и идём на шаг 2.

  5. Стоп.

    В книге Гилл Ф., Мюррей У., Райт М. “Практическая оптимизация”. М.: Мир, 1985 при изложении теории оптимизации приводятся условия локального минимума функции в точке *

    А) Необходимые:

    1. Вектор градиента g( *) = 0 ;

    2. матрица вторых производных (матрица Гессе (Hesse), гессиан, играющая в оптимизации ключевую роль) H 0 в точке *.

    Б) Достаточные:

1. g( *) = 0 ;

2. H > 0 в точке *.
Поэтому в качестве условий останова обычно фигурируют:

а) малость изменения (относительного) целевой функции  (Q( (l)) – Q( (l+1))) / Q( (l))  1 ;

б) малость относительного изменения параметров  ;

в) малость градиента Q (т.к.  Q / = 0 в стационарной точке)  || gl+1 || = 3 .

Ни одно отдельно взятое условие останова не гарантирует достижения минимума. Так, если R1 , для Q(), имеющей горизонтальную точку перегиба , в этой точке выполнены условия останова a) и в), а если эта Q() быстро убывает, то в точке , расположенной на таком участке на склоне, выполнено условие останова б), но ни та, ни другая точка не является точкой локального минимума. Надёжным условием останова является совместное выполнение всех трёх критериев (при этом нужно согласовывать допуски 1, 2, 3   ).

Вернёмся к рассмотрению допустимого шага. Для заданного направления рассмотрим f() = Q( + ) (индекс итерации для краткости опускаем). Пусть мы ещё не в точке минимума, так что градиент не равен нулю, || g || > 0. Имеем

= g T . (8)

При рассмотрении малой окрестности точки можно записать
f() = Q( + ) f(0) +     | = 0 = Q( ) +  g T .

Чтобы шаг был допустимым, требуется, чтобы Q( + ) < Q( ), т.е.  g T должно быть меньше 0 при достаточно малых . Так как > 0, отсюда следует условие

g T < 0 (9)

(в допустимом шаге направление шага составляет с градиентом тупой угол).

Теорема 4.1. Для того, чтобы шаг  был допустимым, необходимо и достаточно выполнение условия

= R g , (10) где R R k x k , R > 0 .

Доказательство (а) Достаточность . Покажем, что из (10) следует (9). Поскольку || g || > 0 , имеем

g T = g TR g < 0 .
(б) Необходимость. Пусть g T < 0 . Покажем, что

RI k (11)  требуемая матрица. Действительно, R = R T очевидно, и выполняется равенство R g = g + g+ = , т.е. (10). Покажем, что для всех zR k , || z || > 0 квадратичная форма z TR z > 0. Действительно,

z TR zzTz .

­ –––––––––– ––––––––– Первое подчёркнутое слагаемое   согласно неравенству Коши-Шварца : (|| z|| 2 || g || 2 - (z T g) 2 ) / || g || 2   , причём равенство реализуется при z = g. Но в этом случае второе слагаемое равно - 2( Tg) > 0 (согласно (9)). Напротив, если второе слагаемое z T равно 0, то z ортогонален , но при этом первое слагаемое больше 0, ибо оно равно 0 лишь при z = g , но тогда было бы g , что противоречит (9). 

Итак, различные допустимые методы отличаются видом матрицы R и способом выбора величины шага на итерациях. Таким образом,
(l+1) = (l) l R l g l , где l > 0, R l > 0 . (12)
1. Метод наискорейшего спуска. Согласно (10), = R g , с R > 0 , || || = 1. Простейший способ выбора R – положить её единичной матрицей, с точностью до константы:

R = = , (10’)

т.е. направление шага выбирается по антиградиенту. Это метод наискорейшего спуска, предложенный Коши 150 лет назад. Величина шага в нём произвольна:

= const > 0  (l+1) = (l) l g l . (12’)

Метод является допустимым, и значения целевой функции вдали от минимума уменьшаются довольно быстро. Но вблизи от точки * для функций с сильно вытянутыми линиями уровня, например, для функции (долины) Розенброка:

Q() = 100 (2 1 2 )2 + (1 - 1 )2

итерации наискорейшего спуска приводят к шагам от борта к борту долины, так что количество итераций до достижения * = Arg min Q() велико, сходимость крайне медленная. Поэтому метод наискорейшего спуска редко используют при решении практических задач.

2. Метод Ньютона. Недостаток метода наискорейшего спуска в том, что он “видит ближайшую перспективу”, даваемую градиентом. Между тем, как вы знаете из теории дифференцируемых функций одной или нескольких переменных, поведение графика (поверхности) функции, радиус её кривизны определяются вторыми производными. Рассматривается квадратичное приближение для Q( ) :
Q( )   Q( (l)) + ( (l)) T g l + 1/2( (l)) TH l ( (l)) , (13)

где H lR k x k  матрица вторых производных Q( ) (матрица Гессе ). В качестве (l+1) принимается стационарная точка , т.е. решение уравнения Эйлера

0 =  / gl + H l ( (l))  (l+1) = (l)H l -1 gl . (14)

Шаг (14) будет допустимым, если H l > 0 . Но достаточным условием минимума является H l > 0 в точке минимума *, а так как H ( ) непрерывна, то в близкой окрестности * обеспечено H l > 0 . Таким образом, метод Ньютона будет допустимым вблизи * с R = H l -1, l = 1.

Его достоинства:

1. Для квадратичной функции  он сходится за одну итерацию.

2. Для не квадратичных функций он имеет максимальную скорость сходимости (квадратичную) среди общих градиентных методов [2] :

||  (l+1)||  const ||  (l)|| 2.

Недостатки:

1. Не гарантировано выполнение условия H l > 0 вне близкой окрестности *.  Нужны хорошие начальные оценки.

  1. Вычисление в итерациях матрицы Hl часто более трудоёмко, чем её обращение (или решение системы (14)), т.к. ищем вторые производные h( )  всего k(k+1)/2 производных в точках x1 , ..., xn , n k .

Чтобы не вычислять вторые производные целевой функции, но в то же время использовать быструю скорость сходимости метода Ньютона, в общих методах оптимизации были разработаны квази-ньютоновские методы, первый из которых – метод Дэвидона-Флетчера-Пауэлла (или переменной метрики). В этих методах последовательно (на итерациях) строятся приближения E(l) к матрице Гессе H , причём переход к следующему приближению E(l+1) производится путём обновления (модификации) матрицами С (l) ранга 1 или ранга 2 по формуле

E(l+1) = E(l) + С (l) ,
где, например, для случая поправок ранга 1
С (l) = Const (b(l) - E(l) s(l)) (b(l) - E(l) s(l))T,

а

s(l) = (l+1) (l) =  (l), b(l) = g (l+1) - g (l)



(обновление E(l+1) делается уже после вычисления нового приближения (l+1)).

Таким образом, для вычисления очередного приближения E(l+1) требуются только первые производные целевой функции Q( ) для построения вектора b(l). Показано, что последовательность матриц {E(l)} сходится к H за конечное, довольно малое, число шагов. В качестве начального приближения можно брать единичную матрицу: E(1) = Ik . Из (14) ясно, что на каждом итерационном шаге метода Ньютона нужно вычислять обратную к матрице H матрицу, что требует O(n 3) флоп работы. При использовании вместо неё матрицы E(l+1) , например. в случае поправок ранга 1, можно использовать формулу Шермана-Моррисона (см., например, Голуб Г., Ван Лоун Ч. “Матричные вычисления” М.: Мир, 1999. Стр. 58): если матрица A R m x m , а u , v – векторы длины m, то
(A + u v T ) – 1 = A – 1 A – 1 u (1 + v T A – 1 u) – 1 v T A – 1.


Если матрица A – 1 (у нас E(l)) уже вычислена, то обращение матрицы A + u v T (у нас E(l+1)) по формуле Шермана-Моррисона потребует всего лишь O(n 2) флоп. Таким образом, при использовании квази-ньютоновских методов не только не нужно вычислять вторые производные функции регрессии, но и обращение матрицы E(l+1) можно делать экономично.

Более подробно материалы по теории и (что особенно ценно !) по практике оптимизации можно посмотреть в превосходной книге Гилл Ф., Мюррей У., Райт М. “Практическая оптимизация”. М.: Мир, 1985. При изложении метода МП в конфирматорном факторном анализе мы ещё вспомним квази-ньютоновские методы...

Название “метод Ньютона” связано с задачей поиска изолированного корня x0 функции f(x), т.е. в простейшем одномерном случае ищется решение уравнения

f(x) = 0. (*)

Ньютон использовал (если излагать метод в современных обозначениях) линейное разложение по формуле Тейлора в окрестности точки t

f(x) = f(t) +(x - t) = 0 ,

где точка x* расположена между x и t. Отсюда, если производная – гладкая функция, а точки x и t близки, можно взять производную в точке x = t , а не в x* и получить очередное приближение для корня x0 в виде:

x = t - . (**)

В многомерном случае x вектор, матрица, а операция деления заменяется в (**) на обращение этой матрицы и умножения результата обращения на x.

Как мы видели, среди и необходимых, и достаточных условий минимума (само по себе оно – признак стационарной точки) – 0 = g( *) = Q / . Полагая в формуле (*) f = g (так что вместо нужно подставить H = ), придём к итерационному шагу типа (**) метода Ньютона, совпадающему с (14).
  1   2

Похожие:

Нелинейный метод наименьших квадратов iconМетод наименьших квадратов. Регрессионный и корреляционный анализ в медицинских исследованиях
Метод наименьших квадратов используется для расчетов параметров функции заданного вида, наилучшим образом отражающую зависимость...
Нелинейный метод наименьших квадратов iconВопрос №4: Предпосылки метода наименьших квадратов, гомоскедастичность, гетероскедастичность, понятие о методе максимального правдоподобия
Метод наименьших квадратов — один из методов регрессионного анализа для оценки неизвестных величин по результатам измерений, содержащих...
Нелинейный метод наименьших квадратов iconТеория математической обработки геодезических измерений. Метод наименьших квадратов

Нелинейный метод наименьших квадратов iconЛабораторная работа по теме: метод наименьших квадратов дисциплина «Численные методы»

Нелинейный метод наименьших квадратов iconВопросы к экзамену по дисциплине «Численные методы»
Метод наименьших квадратов. Нахождение приближающей функции в виде линейной функции
Нелинейный метод наименьших квадратов iconМетод Гаусса в математике Историческая справка
Гаусса), теории электричества и магнетизма, геодезии (разработка метода наименьших квадратов) и многих разделов астрономии
Нелинейный метод наименьших квадратов iconТема фильтры сглаживания. Метод наименьших квадратов
Не перестаю удивляться дерзкой гениальности Стефенсона и братьев Черепановых. Как они отважились построить паровоз, не располагая...
Нелинейный метод наименьших квадратов iconМетод наименьших квадратов
Здесь будет рассмотрена полиномиальная аппроксимация. Это означает, что наша задача состоит в том, что, опираясь на начальные данные...
Нелинейный метод наименьших квадратов iconМатематические методы обработки наблюдений
Данное пособие предназначено для студентов старших курсов астрономического отделения математико-механического факультета спбГУ. В...
Нелинейный метод наименьших квадратов icon3 основы эконометрических методов
Кратко обсудим основные проблемы этой области экономической теории, а затем рассмотрим один из наиболее часто используемых эконометрических...
Разместите кнопку на своём сайте:
ru.convdocs.org


База данных защищена авторским правом ©ru.convdocs.org 2016
обратиться к администрации
ru.convdocs.org