Диссертация (1150844), страница 10
Текст из файла (страница 10)
Обоснованием эвристики послужит следующая лемма.Лемма 3.2.3. Пусть 1 , 2 — две ограниченные вещественные функции, определённые на компакте , и 0 такое, что 1 (0 ) = max∈ 1 () и 1 (0 ) ≤2 (0 ). Определим () = min(1 (), 2 ()). Тогда 0 = arg max∈ ().Доказательство. Предположим противное. Пусть существует такое ^ ∈ , что (^) > (0 ). Тогда 1 (^) ≥ (^) > (0 ) = 1 (0 ) = max∈ 1 (). Получилипротиворечие.(0)Эвристика состоит в следующем: вычислим каждое | ( )|, т.е. | | в(0)точках (0), и выберем 0 = arg min | ( )|.
Будем поддерживать мно()жество : если индекс входит в , то в точке вида будет вычисляться полином . Инициализируем множество = {0 − 1, 0 , 0 + 1}, где сложение и вычитание индексов подразумевается по модулю . Дополнительно()поясним, что () = { | ∈ }. Добавление соседних точек связано с2 -периодичностью()функции min | ( )| по . Дальше найдём решение подзадачиmin | ()|.(3.18)min | ()| = min | ()|,(3.19)^ = arg max−/ ≤</ ∈ ()Если выполнено следующее условие:∈ (^)∈(^)65то при 1 () = min∈ () | ()|, 2 () = () = min∈() | ()| с учётомлеммы 3.2.3, мы получаем, что решение задачи (3.17) найдено, иначе нужно(^)вновь найти индекс 1 , в котором значение | (1 )| достигает минимума, добавить 1 с соседними точками 1 −1 и 1 +1 в множество и повторить процедуруещё раз.
Равенство в условии (3.19) связано с тем, что 2 () ≤ 1 () при любом.Ниже выписан формальный алгоритм решения задачи (3.17) с помощьюэвристики, которая может приводить к значительному уменьшению трудоёмкости.Алгоритм 3.2.5 (Поиск оптимального поворота). Вход: число , вектор коэффициентов ∈ R+1 .(0)1:Вычислить 0 = arg min0≤≤ −1 | ( )|.2:Положить = {0 − 1, 0 , 0 + 1}, = 0.3:while TRUE do4:Численно решить оптимизационную задачу (3.18).5:if выполнено условие (3.19) then6:return ^ — решение оптимизационной задачи (3.17).7:end if8:Положить = + 1.9:Вычислить +1 = arg min | ( )|.10:11:(^)Положить = ∪ {+1 − 1, +1 , +1 + 1}.end whileВ наихудшем случае, алгоритм 3.2.5 добавит во множество все индексы, и гарантированно завершит свою работу, так как условие (3.19) будет всегдаверным. В наилучшем случае, алгоритм потребует всего 2 вычисления функцииmin∈() | ()|, что является большой оптимизацией при больших .
На практике не встречался случай, когда алгоритм 3.2.5 сделал бы больше 3 итераций.66Тем не менее, реальное число итераций может варьироваться в зависимости отколичества корней у полинома ().3.3. Алгоритмы решения задачи3.3.1. Вычисление взвешенных проекторов с заданным базисомКак в модифицированном методе Гаусса-Ньютона (MGN) (3.6), так и вметоде Variable projection Гаусса-Ньютона (VPGN) (3.4) вычисляются проекциивида ΠZ,W , матрица Z лежит в R × или R ×2 . Мы будем считать, что еслиматрица W — (2 + 1)-диагональная, то она представлена своим разложениемХолецкого W = CT C, C — верхнетреугольная с ( + 1) диагональю; если же̂︀ −1 (Ĉ︀ −1 )T , где W−1 =W−1 — (2+1)-диагональная, то мы представляем W = Ĉ︀ T Ĉ︀ представлена аналогично своим разложением Холецкого.CЗамечание 3.3.1. Как показано в разделе 1.3.1, вычисление выражения видаB† сводится к решению задачи обычного МНК, для чего можно использовать QR или SVD разложение матрицы B.Алгоритм 3.3.1 (Вычисление ΠZ,W X).
Вход: матрица Z, временной ряд̂︀ T C.̂︀X, разложение W = CT C или W−1 = C1:2:if W — (2 + 1)-диагональная thenВычислить вектор = CX, матрицу B = CZ.3:end if4:if W−1 — (2 + 1)-диагональная then5:̂︀ −1 )T X, матрицу B = (Ĉ︀ −1 )T Z.Вычислить вектор = (C6:end if7:Вычислить = B† , см.
замечание 3.3.1.8:return Вычислить Z.673.3.2. Вычисление Π(),W XВычисление Π(),W X можно проводить с использованием особенностейподпространства () (см. раздел 3.2.1) и без них, как было предложено в [17].Начнём с алгоритма, предложенного в статье [17]. Вычисление Π(),W Xиспользует формулу (3.7), для которой нужно вычислять Γ() и Γ−1 ().
Мыприводим следующий алгоритм в том виде, как он использован в статье [17], сбыстрым вычислением Γ() и обратной к ней:Алгоритм 3.3.2 (Вычисление решения Π(),W X методом из [17]). Вход:̂︀вектор , временной ряд X, матрица C.1:̃︀ = CQ().̂︀Вычислить ( + + 1)-диагональную матрицу C2:̃︀ T C.̃︀Вычислить (2 + 2 + 1)-диагональную матрицу Γ() = C3:Вычислить ( + + 1)-диагональное разложение Холецкого матрицы Γ()для умножения Γ−1 () справа на вектор.4:return Вычислить Π(),W X по формуле (3.7).Разработанная в разделе 3.2.1 теория позволяет построить другой способвычисления проекции.Алгоритм 3.3.3 (Вычисление решения Π(),W X с помощью свойств ()).Вход: вектор , временной ряд X, разложение W = CT C или W−1 =̂︀ T C.̂︀C1:Вычислить базис Z(), используя алгоритм 3.2.1 или 3.2.3 в зависимостиот требуемой точности.2:return Вычислить ΠZ(),W X, используя алгоритм 3.3.1.Ключевая разница между алгоритмами 3.3.2 и 3.3.3 состоит в замене вычисления образа () оператора Q() на вычисление его ядра ().683.3.3.
Алгоритм VPGNНиже мы приводим алгоритм Гаусса-Ньютона для метода Variable Projection с итерациями (3.4), который был получен в [17] с использованием приципа Variable projection. Заметим, что быстрая реализация предлагаемого в[17] алгоритма вычисления ΠZ(),W X доступна только когда W−1 — (2 +1)-диагональная.Алгоритм 3.3.4 (Алгоритм Гаусса-Ньютона с помощью принципа Variable Projection и вычисления S* из [17] (VPGN)).
Вход: временной ряд X, вектор коэффициентов ОЛРФ 0 , критерий остановки STOP.1:(0)Вычислить = arg max | |.(0)2:Вычислить S0 = S* (( ) ).3:Положить = 0.4:while не STOP do5:()Вычислить JS* (( ) ) согласно формуле (3.8), используя п. 1–3 алгоритма 3.3.2.6:(︁)︁†()(X − S ), используя алгоритм 3.3.1.Вычислить Δ = JS* (( ) )W7:()Найти = arg min0≤≤1 ‖X − S* (( ) + Δ )‖W , используя численныйметод минимизации на отрезке.(+1)()= ( ) + Δ , +1= .8:Положить ( )9:Вычислить S+1 = S* (( ) ).10:(+1)Положить = + 1.11:end while12:return Временной ряд S .Заметим, что вычисление S* (( ) ) = Π(( ) ),W можно производить однимиз алгоритмов 3.3.2 или 3.3.3. В [17] был использован алгоритм 3.3.2.693.3.4. Модифицированный алгоритм Гаусса-Ньютона MGNНиже мы приводим схему предлагаемого модифицированного метода ГауссаНьютона с итерациями в форме (3.6).Алгоритм 3.3.5 (Модифицированный алгоритм Гаусса-Ньютона (MGN)).Вход: временной ряд X, вектор коэффициентов ОЛРФ 0 , критерий остановки STOP.1:(0)Вычислить = arg max | |.(0)2:Вычислить S0 = S* (( ) ).3:Положить = 0.4:while не STOP do5:Вычислить Z((() )2 ), используя алгоритм 3.2.2 или 3.2.4 в зависимости от требуемой точности для вычисления проекции.6:Вычислить Δ = M† QT (() )Π((() )2 ),W (X − S ), используя алгоритм3.3.1.7:()Найти = arg min0≤≤1 ‖X − S* (( ) + Δ )‖W , используя численныйметод минимизации на отрезке.(+1)()= .= ( ) + Δ , +18:Положить ( )9:Вычислить S+1 = S* (( ) ).10:(+1)Положить = + 1.11:end while12:return Временной ряд S .Подробно возможный критерий остановки STOP рассматривается в разделе 5.1.2.3.3.5.
Сравнение вычислительной трудности алгоритмовСравним варианты алгоритмов 3.3.4 и 3.3.5.70Сравниваемые алгоритмыМы рассматриваем три фиксированных алгоритма:1. Метод Variable Projection Гаусса-Ньютона (VPGN) с вычислением проектора через алгоритм 3.3.2 и шага через алгоритм 3.3.4.2. Полуустойчивый метод Variable Projection Гаусса-Ньютона (S-VPGN) свычислением проектора через алгоритм 3.3.3 и шага через алгоритм 3.3.4.3. Устойчивый модифицированный метод Гаусса-Ньютона (MGN) с использованием алгоритма 3.3.3 и вычислением шага через алгоритм 3.3.5.Замечание 3.3.2.
Алгоритм MGN допускает более устойчивую реализациюв случае сигналов, управляемых ОЛРФ с кратными корнями у характеристического многочлена (в частном случае, для полиномиальных сигналов). Этообъясняется тем, что предложенный метод MGN работает с матрицами счислом обусловленности ( ) в методе (см.
теорему 3.2.1), где — наибольшая степень корня характеристического полинома, в то время, как методVPGN оперирует с матрицами с числом обусловленности ( 2 ) [17, Section6.2].Численное сравнение упомянутых алгоритмов приведено в разделе (5.1.2).Сравнение трудоёмкости алгоритмов 3.3.2 и 3.3.3Найдём вычислительную сложность алгоритмов при → ∞. Заметим,что трудоёмкость поиска оптимального поворота (оптимальное 0 ) в шаге 1алгоритма 3.2.1 в асимптотику алгоритма 3.3.3 не включена.∙ W−1 — (2 + 1)-диагональная.71Асимптотическая сложность алгоритма 3.3.3 — ( 2 + 2 + log ):алгоритмы 3.2.1, 3.2.3 — ( log + 2 ), вычисление проекции через алгоритм 3.3.1 — ( 2 + 2 ), МНК через QR-разложение требует( 2 ) времени.Асимптотическая сложность алгоритма 3.3.2 — ( 2 + 2 ).∙ W — (2 + 1)-диагональная, > 0.Асимптотическая сложность алгоритма 3.3.3 — ( 2 + 2 + log ),оценка времени работы такая же, как и в предыдущем случае.Нет реализации алгоритма 3.3.2 метода с линейной по асимптотикойиз-за того, что Γ(A) — не ленточная.Сравнение трудоемкости алгоритмов 3.3.5 и 3.3.4Ниже мы приведём асимптотики трудоёмкости одной итерации алгоритма3.3.5 и метода 3.3.4.∙ W−1 — (2 + 1)-диагональнаяАсимптотическая сложность итерации алгоритма 3.3.5 — ( 2 + 2 + log ): алгоритмы 3.2.2, 3.2.4 — ( log + 2 ), вычисление проекции через Алгоритм 3.3.1 — ( 2 + 2 ); МНК через QR-разложениетребует ( 2 ) времени.Асимптотическая сложность итерации алгоритма 3.3.4 — ( 2 + 2 ):учитывая, что разложение Холецкого матрицы Γ( ) требуется считатьодин раз за итерацию за ( 2 + 2 ), умножение Γ−1 ( ) на векторстоит ( ( + )).∙ W — (2 + 1)-диагональная, > 0.72Асимптотическая сложность итерации алгоритма 3.3.5 — ( 2 + 2 + log ), оценка времени работы такая же, как и в предыдущем случае.Нет реализации алгоритма 3.3.4 с линейной по асимптотикой из-за того,что Γ(A) — не ленточная.Таким образом, если W−1 — (2 + 1)-диагональная, то трудоёмкость предложенного метода MGN выше на ( log ), чем у метода VPGN.
Однако, еслиW — (2 + 1)-диагональная (а это случай шума-авторегрессии, что являетсявполне естественным предположением), то трудоёмкость метода MGN на порядок ниже.73Глава 4Матричный подход к задаче аппроксимациивременного ряда рядами конечного ранга (методКэдзоу)В этой главе мы рассмотрим алгоритм численного решения задачи HankelSLRA (4) — метод итераций Кэдзоу [1]. Определим следующую последовательность матриц, начинающуюся с заданной матрицы X:X0 = X,X+1 = Πℋ Πℳ X , ≥ 0,(4.1)где ℋ = ℋ, — пространство ганкелевых матриц размера × , ℳ ⊂ R×— множество матриц ранга, не превосходящего , L ∈ R× , R ∈ R× —положительно определённые матрицы весов Πℋ и Πℳ — проекторы на соответствующие множества по матричной норме ‖ · ‖L,R , порождённой скалярнымпроизведением (5).В этой главе мы исследуем несколько вопросов, связанных с методом Кэдзоу.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
