Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 47
Текст из файла (страница 47)
Начальное присваивание: S̃ = S(t−i ); s̃ = ŝ ti ).Б. m-кратное повторение процедуры «скалярного» обновления:pq := S̃ T h , α := q T q + 1 , γ := 1/(1 + 1/α) ,K := S̃q/α , Ŝ := S̃ − γKq T ,ŝ := s̃ + K(z − hT s̃)с экстраполяцией между повторениями S̃ := Ŝ ,s̃ := ŝ.ŝ(t+В. Завершающее присваивание: S(t+i ) := ŝ,i ) := Ŝ ,Tгде h — j-й столбец матрицы H (ti); z — j-й элемент вектора z(ti )(j = 1, 2, . .
. , m).Замечание 13.17. В рассматриваемой прикладной задаче можетбыть использован также расширенный фильтр второго порядка, как это сделано в работе [95]. Важно отметить, что для придания своему фильтру численной устойчивости авторы [95] также приводят его к квадратно-корневойформе.13.13Пример задачи с мультиколлинеарностьюПусть требуется решить несовместную систему Ax ≈ z при121 121 1A = .. .
. ,. .. .. N 1 11112 2 z = .. . . .. N N1(13.19)по критерию наименьших квадратов: kz−Axk2 = (z−Ax)T (z−Ax) → minx .Этот пример здесь выбран потому, что он иллюстрирует «тяжелый» случай мультиколлинеарности: столбцы матрицы A в (13.19) не просто коллинеарны, они совпадают! Данный пример в его статистической интерпретации эквивалентен решению задачи оптимального оценивания неизвестногопараметра x̂◦N по экспериментальным даннымz = Ax̂◦N + v ,300v ∼ N (0, I) ,(13.20)13.13 Пример задачи с мультиколлинеарностьюгде N (0, I) — нормальное распределение T с нулевым средним значением иединичной матрицей ковариаций: E vv = I.Упражнение 13.5. Докажите, что в данном примере нормальная система AT Ax̂N = AT z относительно МНК-решения x̂N ∈ Rn (здесь n = 2)имеет вид N +1 11 1(13.21)x̂N =11 12с общим решением N +1 11◦◦x̂N = x̂N + γ, x̂N =, ∀γ ∈ R ,(13.22)−114где x̂◦N есть нормальное псевдорешение, — единственное МНК -решение сминимальной нормой среди всех МНК -решений x̂N , γ — любое число.Пригоден ли здесь метод нормальных уравнений? Очевидно — нет, таккак нормальная система (13.21) вырождена.1.
По поводу метода нормальных уравнений. Даже если нормальная система не является вырожденной, она может оказатьсяплохообуслов1 1ленной. В данном примере она вырождена: Λ = N. Кроме1 1того, формирование нормальной системы может потребовать храненияи затем одновременной обработки большого массива данных. Например, при обработке данных в астрономии или спутниковой геодезии N(число экспериментальных данных в векторе z) очень велико — десяткиили сотни тысяч — при том, что число оцениваемых параметров n ввекторе x сравнительно мало — несколько десятков [68]. Единовременная обработка сверхбольших массивов данных A z в таких задачах,очевидно, нецелесообразна. Главное, однако, не в этом, а в особенностяхчисленного решения задач с мультиколлинеарностью.2.
По поводу численного решения задачи. Решение нормальной системыв ситуациях, близких к мультиколлинеарности, любым одновременнымметодом (см. разд. 12, стр. 264) не дает результатов. Так, в данномпримере (13.19) метод Хаусхолдера, примененный к нормальным уравнениям (13.21), приводит их к эквивалентному виду√ √ √N +1 − 2− 2 − 2x̂N =,(13.23)0002с общим решением (13.22), но не позволяет выделить x̂◦N .30113 Устойчивые алгоритмы фильтрацииМожет ли в подобных задачах помочь последовательная обработка, изложенная выше и принципиально возможная как в информационной форме(см. подразд.
11.6, стр. 244), так и в ковариационной форме (см. подразд. 11.7–11.8 и далее подразд. 13.3–13.8)?Ai,Введем обозначения: A1 = [1 1], z1 = [1], aTi = [1 1], Ai+1 =aTizizi+1 =, ζi = i, i = 1, 2, . . . . Тогда AN = A, zN = z. Информационныйζi+1последовательный алгоритм для i = 1, 2, . . .
имеет вид: 000Λi = Λi−1 + ai aTi , Λ0 =,di = di−1 + ai ζi , d0 =.0 00 N (N + 1) 1Тогда ΛN = Λ, dN =, и при i = N получаем нормальную12систему (13.21) в виде ΛN x̂N = dN . Отсюда видно, что при таком Λ0 этотинформационный последовательный алгоритм также не решает проблему.Упражнение13.6.Докажите,что если взять Λ0 = ε2I, то Λi =i + ε2ii(i + 1) 1, и нормальная система имеет вид=2 , di =1ii+ε2 1ε→0i(i+1)i + 1 1 = x̂◦.Λi x̂i = di, следовательно, x̂i = Λ−1−→ii di =142(2i + ε2 ) 1Убедитесь, что прием, примененный в этом упражнении 13.6, фактическисовпадает с ридж-оцениванием [18] (гребневое оценивание), т. е. с методомрегуляризации Тихонова [67, 81].
Проблема, которая здесь остается, — выбормалого ридж-параметра ε2 . Из этого упражнения видно, что хотя рекуррентный информационный алгоритм и прост в реализации, и формально даетрешение, стремящееся к нормальному псевдорешению x̂◦i при ε → 0, для явного его получения не избежать решения нормальной системы. При маломε2 это снова создает проблему. Ридж-оценивание ее не решает. Выход можноискать в рекуррентном ковариационном алгоритме МНК.Упражнение 13.7. Возьмите к качестве первого рекуррентного ковариационного алгоритма МНК скаляризованную форму фильтра Калмана−2(см.
подразд. 13.3) с начальной матрицей ковариации P0 = Λ−10 = ε I прималом ε (в пределе ε → 0). Получите этот алгоритм в виде:30213.13 Пример задачи с мультиколлинеарностьюТаблица 13.1. Скаляризованный фильтр Калмана в задаче регрессионного моделирования с мультиколлинеарностью регрессоровВеличина Зависимость от εαiKiPix̂i2i + ε22(i − 1) + ε2 11212i + εP0i + ε2−i−i i + ε22i + ε2 i(i + 1)122(2i + ε ) 1Предел при ε → 0То же при i = Ni, i≥2i−1 1 12i 1P01 −12 −1 1 i+1 114NN −1 112N 1P01 −12 −1 1 N +1 1= x̂◦N14Для i = 1, 2, .
. . выполнять:αi = aTi Pi−1ai + 1 ,Ki = Pi−1 ai /αi ,(13.24)Pi = Pi−1 − Ki aTi Pi−1 ,x̂i = x̂i−1 + Ki(ζi − aTi x̂i−1) .(13.25)Упражнение 13.8. Докажите результаты работы этого алгоритма, сведенные в табл. 13.1.Таким образом, алгоритм (13.24), (13.25) не требует ни решения систем,ни обращения матриц для задач вида (13.19), (13.20), — независимо от размерностей и конкретных данных, — и при малом ε2 он гарантирует выделение в явном виде решения x̂◦N в составе обшего МНК-решения (13.22).Проблемы выбора малого ε2 не возникает; ε−2 нужно выбирать настолькобольшим, насколько позволяет диапазон вещественных чисел на компьютере, поскольку параметры алгоритма Ki , Pi обычно быстро убывают и стабилизируются [61].Упражнение 13.9.Наилучшие результаты дадут в этой задаче неформулы (13.24), (13.25), а их алгебраический эквивалент, получаемый припереходе от Pi к квадратному корню из Pi .
Аналогично упражнению 13.8,испытайте другие алгоритмы, а именно: стабилизованный фильтр Калмана–Джозефа из подразд. 13.4, фильтр Поттера из подразд. 13.5, LDLT -фильтрБирмана из подразд. 13.7, U DU T -фильтр Бирмана из задания на стр. 311,LLT -фильтр Карлсона из подразд. 13.8 или U U T -фильтр Карлсона из задания на стр. 311.303=λ i→ 2i 11−→∞ + ε 2 ε→−→ 0012i13 Устойчивые алгоритмы фильтрацииξλ1x2}i+14λ20√ i( {z2( i +2i 1 )+ε2)ε→ →1 −2= ε|ξ0¯=∞0ηx1i+14Рис. 13.1. Сечение поверхности критерия качества J(x) = (z − Ax)T (z − Ax) (дляупражнения 13.10) плоскостью уровняИсследование задачи (13.19) будет неполным, если вы не изучите поверхность критерия качества J(x) = (z − Ax)T (z − Ax) и не используете возможность преобразования базиса в пространстве Rn , здесь n = 2 — размерностьоцениваемого вектора x.Упражнение 13.10.
Задайте начальное значение ковариации P0 =−2= Λ−10 = ε I. Возьмите зависимость Pi от ε из табл. 13.1. Найдите собственные значения матрицы Pi . Убедитесь, что они равны:λ1 =1,2i + ε2λ2 =1.ε2Найдите соответствующие собственные векторы. Убедитесь, что они равны: 11v1 =,v2 =.1−1Постройте график на плоскости x = (x1, x2 ) и дайте необходимые обоснования к его построению (рис. 13.1).30413.14 Задание на лабораторный проект № 91ξ1 1Введите преобразование базиса по формуле x? ,=√x.η2 1 −1Убедитесь, что при этом P0? = ε−2I выводится из P0 = ε−2I.
Докажите, чтостандартный ковариационный алгоритм дает на i-м шаге значения√ i(i + 1)21λ1 01????√αi = αi , Ki =,P=,x̄=ii0 λ22i + ε2 02(2i + ε2 ) 0вместо значений во втором столбце табл. 13.1 на стр. 303. Переход к коордиTнатам x? , ξ ηсоответствует методу главных компонент [18]. Видно,что первая компонента выделяется (оценивается)точ√ со все2 возрастающей √¯ностью (по мере i → ∞, ε → 0): ξ = i(i + 1)/( 2(2i + ε )) → (i + 1)/(2 2),а оценка второй компоненты остается равна нулю: η̄ = 0. Проверьте,√ что¯обратное преобразование дает√ правильный результат: x̄1 , (ξ + η̄)/ 2 →→ (i + 1)/4 и x̄2 , (ξ¯ − η̄)/ 2 → (i + 1)/4.13.14Задание на лабораторный проект № 9ВведениеВ данном лабораторном проекте мы предлагаем к изучению современные методы построения численно устойчивых и экономичных по затратамресурсов ЭВМ алгоритмов метода наименьших квадратов (МHК), решающих актуальную задачу последовательного обновления оценок по измерениям.
Как уже отмечалось, эта задача возникает во многих приложениях,например, при оценке состояния (элементов движения) объекта по последовательно поступающим данным наблюдения (см. выше подразд 13.12),при подгонке параметров модели под результаты продолжительных экспериментов, проводимых в физике, астрономии, экономике, бизнесе, социологии,психологии и в других областях, где выявление регрессий, анализ и прогнозирование тенденций опирается на включение в обработку данных наблюдения по мере их поступления, для того, чтобы постепенно идентифицироватьмодель объекта (процесса) или уточнять параметры модели.1. ЗаданиеA.
Спроектировать, отладить и продемонстрировать в действии программу решения несовместной системы уравнений Ax ≈ b, A = A(m, n),30513 Устойчивые алгоритмы фильтрацииm > n, rank A = n, в смысле наименьших квадратов в соответствии с вашимвариантом последовательного алгоритма (см. ниже подразд. 13.15).Б. Оценить результаты решения по трем показателям:(1) погрешность (абсолютная и относительная) решения;(2) затраты основной памяти компьютера на хранение данных, необходимых только для заданного алгоритма;(3) теоретическое и реальное число основных операций компьютера длявыполнения заданного алгоритма.Эти показатели определить в зависимости от следующих параметров задачи:(1) размерность задачи, т.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
