Диссертация (1150844), страница 14
Текст из файла (страница 14)
. . , ⌈/2⌉ do6:Вычислить = ().7:if , − , = 0 then8:Добавить в индекс ограничения вида { − ≥ 0} из (4.16),и перейти к следующей итерации.9:10:11:end ifif , = , thenДобавить в индекс ограничения вида { − ≥ 0} из (4.17),и перейти к следующей итерации.12:end if13:end for14:for = 1, 2, .
. . do15:Решить следующую задачу квадратичного программирования, вычислить значение точки минимума ⋆ :⋆ = arg min (),∈RR = { | T = 0, ∈ },где () определено в (4.11).16:17:18:19:if для точки ⋆ выполняются все условия вида (4.16) и (4.17) thenreturn 0 = ⋆ , 0 = .elseДобавить к рабочему множеству индекс любого ограничения,которое не выполняется.20:21:end ifend for96Заметим, что шаг 5 алгоритма 4.3.3, необходимый для того, чтобы получить корректную начальную точку, совершит не более ⌈/2⌉ итераций вследствие того, что каждая итерация на единицу уменьшает размерность подпространства R , в которой лежит начальная точка ⋆ .Таким образом, при использовании данной эвристики получаем рекурсивный алгоритм, который уменьшает размерность задачи в раз до тех пор, пока не станет достаточно малым для использования первой стратегии.Для решения задачи поиска весов было применено несколько нетривиальных переходов.
Сформулированная в терминах минимизации квадратичногофункционала задача 4.3.1 содержит 2 линейных ограничений, что не позволяет эффективно искать решение на практике. Поэтому была выведена эквивалентная задача 4.3.2, решение которой – минимум из решений /2 задачквадратичного программирования с количеством ограничений порядка . Преимущество формулировки 4.3.2 состоит в том, что перебор задач КП можноостановить заранее, используя разработанный метод остановки, реализованныйв алгоритме 4.3.2. На практике перебор ограничивается уже на первой итерации, вследствие чего решение находится эффективно.4.4.
Поиск весов с помощью минимизации гладкойфункцииВ этом разделе рассмотрим поиск весов при иных предположениях, чемв разделе 4.2. Рассмотрим задачу (4.7), но в случае белого нестационарного2 Tшума, т.е. с W0−1 = Σ = diag(12 , . . . , ) . Как и в разделе 4.2, для упрощения зафиксируем матрицу Σ−1 = I , соответствующую стационарному беломушуму.Для указанного частного случая получим следующую эквивалентную за97дачу:1‖ΣQ, (1 ) − 1 ‖2 .cond R≤1/ 2R⋆nonstat = arg min(4.21)R=diag()Оказывается, данную задачу можно свести к следующей задаче оптимизации в параллелепипеде:R⋆nonstateq⃦⃦2⃦⃦,̃︀1 ⃦ ⟨1 , ΣQ (1 )⟩,̃︀ − 1 ⃦= arg min ⃦ΣQ (1 )⃦ ,̃︀ 2⃦2 ⃦ ‖ΣQ, (1 )‖̃︀max ≤1(4.22)̃︀min ≥где ⟨·, ·⟩ обозначает обычное скалярное произведение.Теорема 4.4.1.
Задачи (4.21) и (4.22) эквивалентны, т.е. R⋆nonstat = R⋆nonstateq .Доказательство. Сделаем замену переменных, введя дополнительный пара̃︀ при этом каждый из элементов ̃︀ ограничен снизу , аметр : пусть = ,сверху 1, таким образом, cond R ≤ 1/. Тогда задача (4.21) принимает вид(︂)︂1,2̃︀ − 1 ‖ .arg min min ‖ΣсQ (1 )(4.23)2̃︀max ≤1̃︀min ≥Применим к полученной задаче принцип Variable projection (см. раздел 1.3.2) поновой переменной , по которой задача минимизации является квадратичной.Мы получим:̃︀⟨1 , ΣQ, (1 )⟩ =,̃︀ 2‖ΣQ, (1 )‖⋆и подстановка ⋆ в (4.23) даёт нам вид (4.22).Замечание 4.4.1.
В формулировке участвует функция (4.22), которая не является квадратичной, но является гладкой в параллелепипеде. Это означает,что для численной минимизации этой функции можно использовать достаточно эффективный метод L-BFGS-B [55].982можно исЗамечание 4.4.2. В случае равных весов 12 = 22 = . . . = пользовать решение задачи (4.21) вместо решения задачи (4.8) с использованием квадратичного программирования. Плюс данного подхода — простота, так как не требуется реализовывать сложные алгоритмы квадратичного программирования, если есть готовая реализация алгоритма L-BFGS-B.Минус подхода — не гарантируется сходимость оптимизационного метода кглобальному минимуму, в отличие от метода квадратичного программирования.4.5.
Быстрая реализация алгоритма Кэдзоу4.5.1. Общая идея быстрой реализацииВ этом разделе мы обсудим быструю реализацию метода итераций Кэдзоу, заданных в (4.1). Оригинальная идея была предложена в работах [56, 5],и представляет из себя следующее: рассмотрим одну итерацию метода Кэдзоу:X+1 = −1 (Πℋ Πℳ (X )), и перепишем результат Πℳ (X ) в виде суммы∑︀ ̂︀̂︀матриц единичного ранга: Πℳ (X ) = =1 X , rank X = 1. Затем мы мо∑︀̂︀ )).
Проекторы Πℳ и Πℋ (вжем вычислить X+1 как X+1 = =1 −1 (Πℋ (Xприменении к матрицам ранга 1) могут быть реализованы эффективно с точкизрения трудоёмкости и потребляемой памяти благодаря использованию быстрого преобразования Фурье. Предлагается расширение оригинальной идеи сослучая фробениусовской нормы на случай матричной нормы, порождённой скалярным произведением (5), в определении которой матрицы L и R являютсяленточными.994.5.2.
Вычисление проектора ΠℳПерейдём к первому шагу итерации метода Кэдзоу — вычислению проекции на множество ℳ . Приведём следующую теорему из [18], которая даётрешение задачи проектирования Πℳ (X ).Теорема 4.5.1. [18, Theorem 1] Рассмотрим следующую задачу: пусть заданаматрица A ∈ R× , требуется вычислить B = Πℳ A по норме, порождённой скалярным произведением (5). Пусть L ∈ R× и R ∈ R× из (5)T×представлены в виде L = OL OT, OR ∈ R× .L , R = OR OR , где OL ∈ R̃︀ = OT AOR , и вычислим её сингулярное разложение: обоВведём матрицу ALзначим первые сингулярных чисел как 1 , . .
. , , сингулярные векторы как∑︀T1′ , . . . , ′ ∈ R , 1′ , . . . , ′ ∈ R . Тогда выполняется B ==1 , гдеT −1 ′−1 ′ = (OTL ) , = (OR ) .Рассмотрим L и R — (21 + 1)- и (22 + 1)-диагональные положительноопределённые матрицы вместе с их разложениями Холецкого OL OTL = L иOR OTR = R, где OL и OR являются (1 + 1)- и (2 + 1)-диагональными нижнетреугольными матрицами соответственно. Далее мы делаем подстановку в̃︀ = Ã︀ = OT A OR при этом мытеорему с A = A = (X ), B = Πℳ (X ), AL∑︀ ̂︀̂︀ = T , т.е.
вполучаем желаемое представление Πℳ (X ) = =1 X с Xвиде суммы матриц ранга 1.Для быстрого вычисления первых сингулярных чисел и векторов матри̃︀ = OT (X )OR можно использовать методы сингулярного разложения сцы ALиспользованием бидиагонализации Ланцоша, подробности см. в [56, 5]. Чтобыиспользовать данный подход, мы должны реализовать быстрое умножение мат̃︀ на произвольный вектор ∈ R справа или ∈ R слева.
Ключеваярицы Aидея предложена в тех же работах [56, 5] и состоит в том, что мы не вычисляем̃︀ явно, а проводим умножение на вектор шаг за шагом, где шаг умноженияA100на (X ) происходит с помощью быстрого преобразования Фурье. Приведёмалгоритм.Алгоритм 4.5.1 (Вычисление OTL (X )OR ). Вход: временной ряд X ,матрицы O , O , вектор .1:̂︀ = OR С.Вычислить вектор 2:̂︀ используя алгоритм [5, Algorithm 1].Вычислить вектор = (X ),3:return Вектор OTL .̃︀ T , так как Ã︀ T =Мы применяем тот же алгоритм и к вычислению AOTR (X )OL .4.5.3. Вычисление проектора Πℋ от суммы матриц ранга 1Теперь перейдём ко второму шагу итерации метода Кэдзоу (4.1) — вычислению Πℋ .Предложение 4.5.1. Для произвольной матрицы X ∈ R× выполнено сле∑︀T−1дующее: Πℋ X = =1 ( ), где = (1 , .
. . , ) , = W , W = L * R, = (⟨X, (1 )⟩L,R , . . . , ⟨X, ( )⟩L,R )T .Доказательство. Рассмотрим постолбчатую векторизацию : R× → R ,переводящую матрицу в вектор. В пространстве R введём взвешенное скалярное произведение с положительно определённой матрицей M ∈ R× такой,что ⟨, ⟩M = ⟨ −1 , −1 ⟩L,R . Таким образом, — изоморфизм.С помощью оператор Πℋ сводится к вычислению взвешенного МНК.Мы можем записать вектор = ([ (1 ) : . . . ( )])†M (X).
Доказательство завершается раскрытием выражения по формуле (1.5), заменой скалярногопроизведения в R на R× и применением теоремы 4.2.1.Применим предложение 4.5.1 для матрицы X =∑︀=1 X̂︀ — результата101̂︀ имеет ранг 1. Получим:применения проектора Πℳ , где матрица X(︃(︃ )︃)︃∑︁∑︁−1−1̂︀X+1 = ΠℋX=W ,=1(4.24)=1̂︀ , (1 )⟩L,R , . . .
, ⟨X̂︀ , ( )⟩L,R )T .где W = L*R, см. следствие 4.2.1, = (⟨X()()Перепишем -й элемент = (1 , . . . , )T в следующем виде:()̂︀ , ( )⟩L,R = tr(LX̂︀ R T ( )) = = ⟨XT′′ T= tr(OL ′ ′T OTR ( )) = ⟨(OL )(OR ) , ( )⟩F .(4.25)Таким образом, мы свели вычисление Πℋ к задаче быстрого вычислениявектора (⟨T , (1 )⟩F ), . . .
, (⟨T , ( )⟩F )T для произвольных ∈ R , ∈ R , решение которой указано в [5, Algorithm 2] (последний шаг 6 намне нужен). Так как W — (21 + 22 + 1)-диагональная, нам надо вычислить(1 + 2 + 1)-диагональное разложение Холецкого OW матрицы W = OW OTW ииспользовать его для быстрого вычисления формулы (4.24).
Заметим, что явно̂︀ = T , вычислять не нужно, потому чтовекторы и , как и матрицу Xони не используются в правой части формулы (4.25).4.5.4. Алгоритм КэдзоуНиже приведем схему быстро реализации алгоритма Кэдзоу.Алгоритм 4.5.2 (Быстрый алгоритм Кэдзоу). Вход: временной ряд X, длинаокна , ленточные матрицы L, R, критерий остановки STOP.1:TВычислить разложения Холецкого матриц L = OL OTL и R = OR OR .2:Вычислить ленточную матрицу W = L * R по формуле (4.5), её разложение Холецкого W = OW OTW.3:Положить = 0.4:Положить X0 = X.1025:6:while не STOP doВычислить , ′ , ′ , = 1, .
. . , — сингулярные числа и векторы матрицы OTL (X )OR , используя быстрый алгоритм сингулярного разложения (например, с использованием бидиагонализации Ланцоша [56, 5])с алгоритмом умножения вектора 4.5.1.7:Вычислить векторы , = 1, .
. . , , из формулы (4.25) с помощью [5,Algorithm 2].8:∑︀−1Вычислить X+1 = (OW )−1 (OT)(W=1 ).9:Положить = + 1.10:end while11:return Y = X .4.5.5. Трудоёмкость быстрой реализацииНиже мы рассмотрим трудоёмкость отдельных шагов алгоритма 4.5.2.Матрицы O и O вычисляются за время (21 ) и (22 ) с использованием (1 ) и (2 ) места в памяти, W вычисляется за время (1 2 log )и с использованием ((1 + 2 ) ) памяти с помощью предложения 4.2.1 и быстрого преобразования Фурье для вычисления свёрток, плюс требуется ((21 +22 ) ) времени и ((1 + 2 ) ) места для вычисления OW .
Одна итерация алго̃︀ требует (( log + 1 + 2 )) времениритма SVD Ланцоша матрицы Aс ( ) дополнительной памяти [56], а вычисление Πℋ требует (( log +1 + 2 ) + (1 + 2 ) ) времени и ( ) памяти.1034.6. Применение метода Кэдзоу к решению задачиаппроксимации временных рядов рядами конечногоранга4.6.1. Выбор метода для поиска весовОбсудим применение алгоритма Кэдзоу (4.1) для решения задачи (3) вместе с методами поиска весов L, R по заданной обратной автоковариационнойматрице Σ−1 = W0 , где W из (3) равна W0 в этом разделе.Мы введём две комбинации методов: назовём Cadzow(quadratic) — методКэдзоу с вычислением весов с помощью алгоритма квадратичного программирования 4.3.1, и Cadzow(box) — метод Кэдзоу с вычислением весов с помощьюминимизации функции (4.22) методом оптимизации в параллелепипеде, например, L-BFGS-B [55].
Проблема в том, что теоретические результаты и алгоритмыполучены только для частных случаев задачи (4.7), поэтому рекомендации ихприменения в более общих случаях требуют объяснения.Рассмотрим следующие сценарии применения:1. W0 = 0 I , 0 > 0 — случай белого шумаВ этом случае решение задачи (4.7) (точнее, её частного случая (4.8)) R*может быть вычислено точно с использованием алгоритма 4.3.1, либо вычислено приближённо с помощью минимизации функции (4.22). Соответственно, стоит использовать метод Cadzow(quadratic), при этом L = I .Метод Cadzow(box) можно рассматривать как более простую в реализации альтернативу при наличии готовой реализации метода оптимизациив параллелепипеде, однако, оптимальность R* не гарантирована.2.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
