Итерационная проекционная схема в реальном времени для решения задачи об общей неподвижной точке и ее приложения
- Авторы: Гибали А1, Теллер Д1
-
Учреждения:
- Ort Braude College
- Выпуск: Том 64, № 4 (2018): Современные проблемы математики и физики
- Страницы: 616-636
- Раздел: Новые результаты
- URL: https://journals.rudn.ru/CMFD/article/view/22278
- DOI: https://doi.org/10.22363/2413-3639-2018-64-4-616-636
Цитировать
Полный текст
Аннотация
В этой работе мы рассматриваем задачу об общей неподвижной точке (CFPP) с демисжимающими операторами и ее частный случай, выпуклую задачу о допустимости (CFP) в вещественных гильбертовых пространствах. Руководствуясь недавними результатами, полученными Ордонесом и др. в работе [35] и в области алгоритмов в реальном времени в общем, например, в [20, 21, 30], где с самого начала нам недоступны целые наборы операторов/множеств, которые затем получаются постепенно, мы предлагаем итерационную схему в реальном времени для решения задач об общей неподвижной точке (CFPP) и выпуклых задач о допустимости (CFP), в которой участвующие операторы/множества появляются со временем. Такая схема способна работать с любыми блоками данных и для любого конечного числа итераций с последовательным переходом к следующему блоку. Схема основана на недавнем результате, описанном в работе Райха и Заласа [37] и известном как процедура модулярного строкового усреднения (MSA). Сходимость схемы следует из [37] и других классических результатов в теории неподвижных точек и области вариационных неравенств, например, [34]. Также в работе представлены вычислительные эксперименты для линейных и нелинейных задач о допустимости в приложении к восстановлению изображений. Они демонстрируют справедливость и потенциальную применимость нашей схемы, например, в условиях реального времени.
Полный текст
ВВЕДЕНИЕ В этой работе мы имеем дело с задачей об общей неподвижной точке (CFPP) и ее частным случаем, выпуклой задачей о допустимости (CFP) в вещественных гильбертовых пространствах H. Даны операторы Ui : H → H при i ∈ I := {1, 2,..., m} с непустыми множествами неподвижных точек. Задача об общей неподвижной точке состоит в нахождении точки x∗ ∈ H, такой, что m x∗ ∈ (] Fix (Ui) . (1.1) i=1 Выпуклая задача о допустимости является частным случаем задачи об общей неподвижной точке. В этом случае имеем m непустых, замкнутых и выпуклых множеств Ci ⊆ H при i ∈ I. Теперь задача заключается в нахождении такой точки x∗ ∈ H, что m x∗ ∈ (] Ci ±= ∅. (1.2) i=1 Qc РОССИЙСКИЙ УНИВЕРСИТЕТ ДРУЖБЫ НАРОДОВ, 2018 616 Очевидно, что если мы выберем Ui = PCi для всех i ∈ I, где PCi обозначает ортогональную проекцию на i-е множество Ci (что будет пояснено далее) в CFPP (1.1), то мы получим CFP (1.2). CFPP и CFP служат основными средствами моделирования при построении многих важных реальных задач, например, при формировании изображений, в сенсорных сетях, при составлении плана лечения лучевой терапией, повышении разрешающей способности и многих других; см., например, [5, 15]. Одна из самых ранних итерационных процедур для решения задач CFPP, см., например, [34], имеет следующий общий вид: выбираем произвольную начальную точку x0 ∈ H; получив текущую итерацию xk, рассчитываем следующую итерацию xk+1 таким образом: xk+1 = T (xk ), (1.3) где оператор T : H→ H фиксирован и зависит от семейства операторов {Ui : H→ H| i ∈ I}. Более общая схема для неподвижной точки позволяет включить семейство операторов k=0 {Tk : H→ H}∞ ; см., например, обобщенный метод Опиала [11, п. 3.6]. Итерационная процедура формулируется следующим образом: выбираем произвольную начальную точку x0 ∈ H; получив текущую итерацию xk, рассчитываем следующую итерацию xk+1 так: xk+1 = Tk (xk ), (1.4) k=0 где семейство операторов {Tk }∞ зависит от {Ui}i∈I и может иметь различные алгоритмические структуры, например: 1. циклическую (с релаксацией): αk ∈ [ε, 2 - ε] при ε > 0: Tk = Ui(k), где i(k) = (k mod m)+ 1; 2. совместную: 1 m 3. композиционную: Tk := m i=1 m Ui; 4. «жадную» (наиболее удаленную): Tk := TT Ui; i=1 Tk := Uik , где ik = argmax dist(·, Fix(Ui)), i∈I здесь dist(·, ·) - это функция расстояния между точкой и множеством. Возвращаясь к выпуклой задаче о допустимости (CFP), хотелось бы обратить внимание на класс проекционных методов. В 1930-х годах Стефан Качмаж [29] и Джанфранко Чиммино [18] представили итерационные проекционные методы для решения систем линейных неравенств Ax � b, где A ∈ Rm×n, x ∈ Rn и b ∈ Rm. Как оказалось, данная задача может быть легко сведена к эквивалентной ей CFP следующим образом: обозначим через Ai и bi i-е строку и элемент A и b, соответственно. Определим понятие полупространства: i := {z ∈ R | A , z � bi} (1.5) и затем получим: H- n i m i Ax � b ⇔ x ∈ (] H-. (1.6) i=1 i Методы Качмажа [29] и Чиммино [18] используют ортогональные проекции (отражения) на полупространства H- последовательно и совместно, соответственно. Этим, собственно, и характеризуется класс проекционных методов, которые, будучи итерационными процедурами, используют проекции различных типов на множества с учетом того факта, что проекция на пересечение множеств являет собой довольно сложную вычислительную задачу, в то время как проекции на отдельные множества намного проще осуществимы. Именно поэтому эти методы успешно применяются во многих реальных приложениях, и даже были названы «Swiss Army knives», см. [6]. В последние десятилетия класс проекционных методов активно развивался, эти методы оказались способными решать общие выпуклые задачи о допустимости (1.2). Также они включают в себя различные алгоритмические структуры, такие как последовательную, совместную, блочноитерационную, усреднение по строкам и т. д., см. [15], а также [10, 11, 17, 22, 23]). Ордонес и др. в работе [35] изучают разреженную и переопределенную систему линейных уравнений Ax = b (A ∈ Rm×n с m ∼ 103 и n ∼ 109) большой размерности, возникающую в области протонной вычислительной томографии (pCT). Исследуя эту задачу, авторы вводят два проекционных метода в режиме реального времени, метод диагонально ослабленных ортогональных проекций (DROP) [2] и метод координатно-усредненных строковых проекций (CARP) [25]. Они экспериментально демонстрируют, что DROP и CARP в режиме реального времени работают гораздо «лучше и быстрее», чем другие проекционные методы, см., например, [36]. Итак, исходя из [35] и результатов, полученных в области алгоритмов в реальном времени, например [20, 21, 30], нашей задачей в этой работе является введение новой итерационной схемы для неподвижной точки типа (1.3) или (1.4) для решения задач об общей неподвижной точке и выпуклых задач о допустимости. Сосредоточим наше внимание на том случае, когда с самого начала нам недоступны целые наборы операторов/множеств, и мы получаем их постепенно. Следовательно, нам необходимо придумать такую итерационную схему в реальном времени, которая будет способна работать с сегментами набора и включать в себя новый набор, когда таковой появится. В CFPP главная идея состоит в том, что операторы Ui при i ∈ I представлены в виде блоков I = I1 ∪ I2 ∪ ... ∪ IM , 1 � M � m, и последовательны по времени. Возьмем процедуру Райха и Заласа [37], модулярное строковое усреднение (MSA), и покажем, как она может быть применена в вышеописанном случае. Некоторые численные эксперименты демонстрируют потенциальную применимость и преимущества предлагаемого метода для линейных и нелинейных задач о допустимости в режиме реального времени. Краткое содержание данной работы таково: в разделе 2 представлены определения и понятия, необходимые в дальнейшем; в разделе 3 вводятся и анализируются новые итерационные схемы в реальном времени для задач о неподвижных точках и задач о допустимости; далее в разделе 4 будет представлено несколько численных примеров для линейных и нелинейных задач о допустимости и будут обоснованы справедливость и потенциальная применимость новой схемы, которая может быть использована, например, для задач в реальном времени; и, наконец, в разделе 5 представлены выводы и направления для дальнейших исследований. 2. ПРЕДВАРИТЕЛЬНЫЕ СВЕДЕНИЯ Всюду далее в этой работе H - это вещественное гильбертово пространство со скалярным произведением ⊗·, ·) и индуцированной нормой 1· 1. Будем использовать записи xk --и xk → x x для обозначения слабой и сильной сходимостей последовательности { k �∞ k=0 к x, соответственно. Вспомним теперь несколько определений и свойств некоторых классов операторов. Их, наряду со многими другими, можно найти, например, в блестящей работе Цегельского [11]. Определение 2.1. Пусть U : H→ H - некоторый оператор. · Множество неподвижных точек оператора U, обозначаемое через Fix(U ), определяется так: Fix(U ) := {x ∈H | U (x) = x}. (2.1) · Оператор U называется срезающим, если для любого x ∈H и любого z ∈ Fix(U ) ⊗z - U (x),x - U (x)) � 0. (2.2) · Оператор U называется нерастягивающим (NE), если для любых x, y ∈H 1U (x) - U (y)1 � 1x - y1. (2.3) · Оператор U называется квази-нерастягивающим (QNE), если для всех (x, q) ∈H × Fix(U ) 1U (x) - q1 � 1x - q1. (2.4) · Оператор U с Fix(U ) ±= ∅ называется ρ-деми-сжимающим (см., например, [32]) при ρ ∈ [-1, 0), если для всех (x, z) ∈H × Fix(U ) 2 2 2 1U (x) - z1 � 1x - z1 - ρ1U (x) - x1 . (2.5) 2 2 · Оператор U называется ρ-сильно квази-нерастягивающим при ρ � 0, если для всех (x, z) ∈ H× Fix(U ) 1U (x) - z1 � 1x - z12 - ρ1U (x) - x1 . (2.6) Если ρ > 0, то оператор U называется сильно квази-нерастягивающим. · При α ∈ [0, ∞] оператор Uα := Id + α(U - Id) называется α-релаксацией оператора U, а α называется параметром релаксации. Легко видеть, что для каждого α ±= 0 Fix(U ) = Fix(Uα). (2.7) · Оператор U называется усредняющим [3] (см. также [9]), если существуют нерастягивающий оператор N : H→ H и число c ∈ (0, 1), такие, что U = (1 - c)Id + cN, (2.8) где Id - единичный оператор в H. · Квази-нерастягивающий оператор U называется деми-замкнутым в точке y ∈ H, если для x любой последовательности { k �∞ k=0 ⊂H имеем xk --x U (xk ) → y =⇒ U (x) = y. (2.9) · Квази-нерастягивающий оператор U называется приближенно стягивающим, если для люk=0 бой ограниченной последовательности {xk }∞ ⊆H имеем lim k→∞ 1U (xk ) - xk 1 = 0 =⇒ lim k→∞ dist(xk, Fix(U )) = 0. (2.10) Подробнее об этом классе операторов можно прочесть, например, в [13]. Следующая теорема является частью нашего анализа и показывает связь между двумя классами операторов, описанных выше. Теорема 2.1. Пусть U : H → H - оператор с неподвижной точкой, и пусть α ∈ (0, 2]. Тогда U - это срезающий оператор, если и только если его Uα является (2 - α)/α-сильно квази-нерастягивающим. Доказательство. См., например, либо [19, Proposition 2.3(ii)], либо [11, Theorem 2.1.39]. Следующий принцип известен как принцип деми-замкнутости [8]. Принцип деми-замкнутости. Пусть H - это вещественное гильбертово пространство, C ⊆H - замкнутое и выпуклое множество, и пусть S : C → H - нерастягивающее отображение. Тогда Id - S деми-замкнут в y ∈ H. Другим важным результатом, необходимым для нашего анализа, является следующее утверждение (см. [14, утверждение 4.1]). Предложение 2.1. Пусть U : H→ H - это квази-нерастягивающий оператор. Тогда верны следующие утверждения: 1. Если U является приближенно стягивающим, то U - Id деми-замкнут в 0; 2. Если dim(H) < ∞ (H конечномерно) и U - Id деми-замкнут в 0, то U является приближенно стягивающим. Определение 2.2. Пусть Ci ⊆H при i ∈ I - это замкнутые и выпуклые множества с непустым пересечением C := (] Ci ±= ∅. i∈I Семейство множеств C := {Ci | i ∈ I} называется ограниченно регулярным, если для любой k=0 ограниченной последовательности {xk }∞ ⊆H имеем lim max dist(xk, Ci) = 0 =⇒ lim dist(xk, C) = 0. (2.11) k→∞ i∈I k→∞ Следующее утверждение (см. [4, утверждение 5.4 (iii), следствие 5.14, следствие 5.22]) дает условия, гарантирующие ограниченную регулярность семейства множеств. Предложение 2.2. Пусть Ci ⊆H при i ∈ I - замкнутые и выпуклые множества с непустым пересечением: m C := (] Ci. i=1 Если выполнено одно из следующих условий: 1. dim(H) < ∞, 2. int (C) ±= ∅, 3. всякое Ci является полупространством, то семейство множеств C := {Ci | i ∈ I} ограниченно регулярно. Теперь вспомним метрические проекции на замкнутые и выпуклые множества. Пусть C ⊂ H. Для каждой точки x ∈ H существует единственная ближайшая точка из C, обозначенная PC (x), и такая, что 1x - PC (x)1 � 1x - y1 для всех y ∈ C. (2.12) Отображение PC : H → C называется метрической проекцией пространства H на C и является нерастягивающим отображением пространства H на C (на самом деле, FNE (следовательно, срезающим), см. [11, теорема 2.2.21]). Отображение PC характеризуется следующими двумя свойствами (см. [24, п. 3]): (следовательно, Fix(PC ) = C) и PC (x) ∈ C (2.13) ⊗x - PC (x) , PC (x) - y) � 0 для всех x ∈ H, y ∈ C, (2.14) и если C - гиперплоскость, то (2.14) превращается в равенство. Другим важным типом проекций является субградиентная проекция (также срезающий оператор, см., например, [11, лемма 4.2.5] и [5]). Такие проекции крайне важны в случаях, когда выпуклое множество C представляется в виде подуровневых множеств выпуклой функции g : H → R, т. е. C := {x ∈H | g(x) � 0}. Теперь, когда мы дали определения различным типам проекций, а также различным классам операторов, необходимо вспомнить два частных типа проекционных методов: последовательный (предложенный Качмажем, также известный как метод последовательных ортогональных проекций (SOP) [1], метод проекции на выпуклые пространства (POCS) [4] и алгебраический способ перестраивания (ART) [26] в линейных случаях) и блочно-итерационные методы (полностью совместные, если есть только один блок, в линейном случае сводятся к методу Чиммино [18]). С этой целью рассмотрим выпуклую задачу о допустимости с непустыми, замкнутыми и выпуклыми множествами Ci ⊆ H при i ∈ I. Следующие определения управляющих последовательностей, ν=0 {i(ν)}∞ , обуславливают порядок, в котором участвуют ортогональные проекции на множества Ci, i ∈ I и, следовательно, указывают структуру алгоритма. Определение 2.3. ν=0 1. Последовательность {i(ν)}∞ называется циклично управляющей, если i(ν) = (ν mod m)+ 1, где m - это количество множеств в (1.2). ν=0 2. Последовательность {i(ν)}∞ называется почти циклично управляющей на множестве I = {1, 2,..., m} если i(ν) ∈ I для всех ν � 0, и существует целое Q � m (называется почти циклической константой), такое, что I ⊆ {i(ν + 1), i(ν + 2),..., i(ν + Q)} для всех ν � 0. (2.15) ν=0 3. Последовательность {i(ν)}∞ называется удаленно управляющей, если она получается указанием i(ν), такого что dist(xν, Ci(ν)) = max{dist(xν, Ci) | i ∈ I}. (2.16) ν=0 4. Последовательность {i(ν)}∞ называется произвольно управляющей, если i(ν) ∈ I выбирается произвольно и определяется независимо в соответствии с заданным распределением вероятности {pi}. Теперь представим последовательный и совместный проекционные методы для решения выпуклых задач о допустимости. Алгоритм 2.1 (Метод SOP). Дано. Пусть x0 ∈H - произвольная начальная точка. Итерационный шаг. Имея текущую итерацию xk, рассчитаем следующую итерацию таким образом: xk+1 = xk + λk (PC i(ν) (xk ) - xk \ , (2.17) где PCi(ν) обозначает ортогональную проекцию на множество Ci(ν), λk ∈ [ε1, 2 - ε2] для всех ν=0 k � 0 и ε1, ε2 > 0. Управляющая последовательность {i(ν)}∞ циклична на I. Для следующего алгоритма необходимо определить некоторые термины. Вектор ω = (ω(i))i∈I называется весовым вектором, если и ω(i) � 0 для всех i ∈ I ω(i) = 1. i∈I Имея весовой вектор ω, можем дать определение выпуклой комбинации: i Pω (x) := ω(i)PC . k=0 Последовательность весовых векторов {ωk }∞ i∈I называется правильной, если для любого i ∈ I существует бесконечно много значений k, для которых ωk (i) > 0. Алгоритм 2.2 (Метод блочного типа). Дано. Пусть x0 ∈H - произвольная начальная точка. Итерационный шаг. Имея текущую итерацию xk, рассчитаем следующую итерацию таким образом: ω xk+1 = xk + λk (P k (xk ) - xk \ , (2.18) k=0 где {ωk }∞ - это правильная последовательность весовых векторов, а λk ∈ [ε1, 2 - ε2] для всех k � 0 и ε1, ε2 > 0. Чтобы проиллюстрировать несколько типов проекционных методов для решения выпуклой задачи о допустимости, ограничимся лишь линейной задачей о допустимости, которая является системой линейных уравнений Ax = b, в которой A ∈ Rm×n, x ∈ Rn и b ∈ Rm. Обозначим через Ai и bi i-е строку и элемент A и b, соответственно, и зададим i-ю гиперплоскость Hi = {z ∈ Rn | Ai, z = bi}. Иллюстрации этих и других методов представлены на рис. 1, взятом из [16]. РИС. 1. Различные проекционные методы для линейного случая (рисунок из [16]). Замечание 2.1. Заметим, что итерации для неподвижной точки (1.3) и (1.4) включают в себя вышеописанные методы. Например, если мы рассмотрим задачу об общей неподвижной точке с Ui = PCi , то мы получим выпуклую задачу о допустимости. Более того, если Tk = Ui(k), где i(k) = (k mod m)+1, то мы получим метод последовательных ортогональных проекций (SOP) (2.1), m а если возьмем только один блок I размера m, то, приняв T = 1 ), Ui, мы получим метод Чиммино в линейном случае и блочный метод в общем случае. m i=1 Далее вспомним две теоремы о неподвижной точке, классическую теорему Опиала [34] и ее обобщение [11, п. 3.6]. Теорема 2.2. Пусть H - вещественное гильбертово пространство, и пусть C ⊂ H - замкнутое и выпуклое множество. Если T : C → C - усредненный оператор с Fix(T ) ±= ∅, то k=0 для любого x0 ∈ C последовательность {xk �∞ , порожденная xk+1 = T (xk ), сходится слабо к точке x∗ ∈ Fix(T ). Приведем обобщенную теорему Опиала, см., например, [11, п. 3.6], позволяющую работать с k=1 семейством операторов {Tk : H→ H}∞ . Теорема 2.3. Пусть C ⊆ H - непустое, замкнутое и выпуклое множество, S : C → H - оператор со множеством неподвижных точек, такой что S - Id - деми-замкнут в 0. Пусть k=1 {Tk : H → H}∞ - асимптотически регулярная последовательность квази-нерастягивающих операторов, таких что Fix(S) ⊆ ∞ (] Fix(Tk ) . k=1 k=0 Пусть {xk �∞ - последовательность, порожденная xk+1 = Tk (xk ), с произвольным x0 ∈ H. k=1 1. Если последовательность операторов {Tk }∞ имеет свойство lim k→∞ 1Tk (xk ) - xk 1 = 0 =⇒ lim k→∞ 1S(xk ) - xk 1 = 0, (2.19) то x { k �∞ k=0 сходится слабо к точке Fix(S). k=1 2. Если H конечномерно и последовательность операторов {Tk }∞ имеет свойство lim k→∞ 1Tk (xk ) - xk 1 = 0 =⇒ lim k→∞ inf 1S(xk ) - xk 1 = 0, (2.20) то x { k �∞ k=0 сходится к точке Fix(S). 3. АЛГОРИТМ В этом разделе мы сосредоточимся на задаче об общей неподвижной точке (CFPP) (1.1) с семейством деми-сжимающих операторов {Ui}i∈I, таких что (] Fix (Ui) ±= ∅. i∈I Ситуация заключается в том, что множество индексов I разбивается на M блоков I = I1 ∪ ... ∪ IM M выбором {mt}t=0 ⊂ Z (где Z - это множество целых чисел), таких что 0 = m0 < m1 < . . . < mM = m, и каждому 1 � t � M соответствует подмножество It := {mt-1 + 1, mt-1 + 2,..., mt}. Это, очевидно, делит семейство операторов {Ui}i∈I на соответствующие группы операторов. Так как наша задача состоит в построении итерационной схемы в реальном времени, сосредоточимся на том случае, когда блоки и соответствующие операторы находятся в нашем распоряжении не с самого начала, а становятся известны постепенно. В недавней работе Ордонеса и др. [35] два метода в реальном времени ((DROP) [2] и (CARP) [25]) для решения систем линейных уравнений Ax = b, где A ∈ Rm×n с m ∼ 103 и n ∼ 109, возникают в области протонной вычислительной томографии (pCT). Недавно Райх и Залас [37] представили процедуру модулярного строкового усреднения (MSA) для решения задачи об общей неподвижной точке в вещественных гильбертовых пространствах. Данная схема является довольно гибкой и позволяет строить вспомогательные операторы Tk, названные модулями, которые могут быть использованы во внутреннем цикле расширенного алгоритма с конечным количеством итераций Nk. Для представления алгоритма введем несколько структур операторов Tk, построенных с помощью семейства операторов {Ui}i∈I с учетом M блоков I = I1 ∪ ... ∪ IM и использующихся во вспомогательном цикле нашего алгоритма. Эти структуры представлены как частные случаи для модулярного строкового усреднения [37]; более подробную информацию (в том числе и исторической справки) см. в работе [37] и имеющейся там библиографии. Определение 3.1. 1. Циклическая (с релаксацией): αk ∈ [ε, 2 - ε] при ε > 0: Tk = Ui(k), где i(k) = (k mod m)+ 1. 2. Выпуклая комбинация: Для весового вектора ωk (i) � 0 для любых i ∈ Ik, таких что ωk (i) = 1, i∈Ik положим 3. Композиционная: 4. Блочная: αk ∈ [ε, 2 - ε] при ε > 0: Tk = ωk (i)Ui. i∈Ik Tk = TT Ui. i∈Ik ⎛ ⎞ Tk = Id + αk ⎝ ωk (i)Ui - Id⎠ . i∈Ik 5. «Жадная» (наиболее удаленная): Tk := Uik , где ik = argmax dist(·, Fix(Ui)). i∈Ik Другими, более частными, структурами, которые могут быть здесь использованы, являются усреднение по строкам и различные типы операторов Дугласа-Рашфорда (см., например, [7]), в случае, если вместо Ui, используется 2Ui - Id. Алгоритм 3.1 (Блочно-итерационная схема в режиме реального времени). Дано. Пусть x0 ∈H - произвольная начальная точка, определено N0 ∈ N (количество итераций), и даны первый блок I1 и соответствующее подмножество операторов {Ui}i∈I1 . Вычислим x1 таким образом: x1 = T0(x0), (3.1) где оператор T0 строится согласно определению 3.1 и может быть циклическим, совместным или композицией {Ui}i∈I1 . Итерационный шаг. Имея текущую итерацию xk и зная Nk ∈ N (количество итераций), вычислим следующую итерацию: xk+1 = Tk (xk ), (3.2) где оператор Tk может быть построен следующим образом. 1. Если k < m: заданы блоки I1, I2,..., Ik и, следовательно, операторы {Ui}i∈I1∪...∪Ik . Тогда Tk могут быть построены с учетом каждого, некоторых или всех операторов {Ui}i∈I1∪...∪Ik в циклической, совместной или композиционной формах, основанных на определении 3.1. 2. Если k � m: заданы сразу все блоки и, следовательно, операторы {Ui}i∈I, и тогда Tk могут быть построены на основании определения 3.1 с учетом всего семейства операторов {Ui}i∈I. 1. Сходимость. Для сходимости нашего алгоритма 3.1 будем считать, что выполнены следующие условия. Условие 3.1. Для всех i ∈ I операторы Ui - деми-сжимающие с Fix(Ui) ±= ∅, и такие, что Ui,α (α-релаксация Ui) являются (2 - α)/α-сильно квази-нерастягивающими. Условие 3.2. I ⊆ Ik ∪ Ik+1 ∪ ... ∪ Ik+s-1 для всех k = 0, 1 ... и некоторого s � m - 1. k=0 Условие 3.3. Последовательность {Nk }∞ , представляющая количество итераций на каждый блок, является ограниченной. В журнале Numerical Algorithms Райх и Залас [37] предложили процедуру модулярного строкового усреднения (MSA) для решения задачи об общей неподвижной точке в вещественных гильбертовых пространствах. Они представили гибкий метод [37, процедура 1.1]) построения вспомогательных операторов Tk, названных модулями, которые могут быть использованы во внутреннем цикле расширенного алгоритма с конечным числом итераций Nk для семейства операторов {Ui}i∈I. Вследствие модульности их схемы, и с учетом условий 3.1-3.3, сходимость нашей итерационной схемы в реальном времени, алгоритма 3.1, вытекает непосредственно из доказательства теоремы 4.1 в работе Райха и Заласа [37], хотя термин «в реальном времени» в этой работе и не упоминается. Следующая теорема является модификацией теоремы [37, теорема 4.1], скорректированной для алгоритма 3.1. Теорема 3.1. Пусть H - вещественное гильбертово пространство, и даны операторы Ui : H → H при i ∈ I, такие что Fix (Ui) ±= ∅. Предположим, что условия 3.1-3.3 выполнены, и k=0 пусть последовательность {xk �∞ сгенерирована алгоритмом 3.1. ∈ 1. Если для каждого i I оператор Ui удовлетворяет принципу деми-замкнутости Опиала, m k=0 то последовательность {xk }∞ сходится слабо к некоторой точке из C = П Fix (Ui) . i=1 2. Если для каждого i ∈ I оператор Ui является приближенно стягивающим, а семейство k=0 C := {Fix(Ui) | i ∈ I} ограниченно регулярно, то последовательность {xk }∞ сходится сильно к некоторой точке из C. Необходимо упомянуть, что в этой работе мы считаем, что для всех i ∈ I операторы Ui являются деми-сжимающим, и, следовательно, по теореме 2.1 для всех i ∈ I и α ∈ (0, 2] мы полагаем Ui,α (α-релаксация Ui) (2 - α)/α-сильно квази-нерастягивающими и, таким образом, срезающими; это, собственно, и используется в алгоритме 3.1. Замечание 3.1. 1. Используя условие 3.1 и теорему 2.1, получаем, что операторы Ui являются срезающими, а структура алгоритма 3.1 может быть произвольного типа согласно определению 3.1. 2. Условие 3.2 означает, что управляющая последовательность почти циклична (определение 2.3 (2)). 3. Условие 3.3 говорит о том, что количество итераций Nk ограничено, что в свою очередь означает, что допустимо любое конечное число промежуточных шагов в каждом блоке. n 1. В случае, если Ui = PCi и H = R , все предположения относительно операторов, касающиеся непрерывности, будь то срезающие, деми-замкнутые или приближенно стягивающие операторы, справедливы. 2. Для задач о допустимости и, в частности, для линейных задач о допустимости произвольная управляющая последовательность оказывается эффективной управляющей последовательностью (это также известно как рандомизированный метод Качмажа, см., например, [31, 33]). Хотя доказательство сходимости не охватывает нашу ситуацию, представляется довольно интересным изучение поведения в теории; это, вероятно, будет исследовано в дальнейших работах. Несмотря на это, мы изучили этот случай в наших численных экспериментах. Другим относящимся к теме результатом, в котором произвольные управляющие последовательности также рассматриваются как класс стохастических алгоритмов (в частности, для задач о допустимости и для вариационных неравенств), является работа Иусема и соавторов [28]. 3. В [12, теорема 4.3 и следствие 4.4] представлено новое обобщение теоремы Опиала с меньшим количеством ограничивающих предположений, чем в теореме 2.3. Поэтому будет интересно наблюдать, как данный результат будет изучен и приложен к результатам этой работы. 4. Ордонес и соавторы в работе [35] представили два проекционных метода в реальном времени, (DROP) [2] и (CARP) [25], для решения разреженных и переопределенных систем линейных уравнений Ax = b большой размерности, в которых A ∈ Rm×n с m ∼ 103 и n ∼ 109, возникающих в области протонной вычислительной томографии (pCT). Эта работа предлагает многообещающие результаты в экспериментальном поведении, но ей, к сожалению, не хватает достаточного математического обоснования. И хотя наш подход не применяется к их схемам, мы способны провести анализ данной гибкой схемы, которая может быть приложена не только к линейным задачам о допустимости. 1. ВЫЧИСЛИТЕЛЬНЫЕ ЭКСПЕРИМЕНТЫ В этом разделе мы будем сравнивать 4 схемы в реальном времени: Чиммино [18] (алгоритм 2.1), Качмажа [29] (алгоритм 2.2), рандомизированный метод Качмажа [31, 33] и «жадный» метод Качмажа для линейных и нелинейных (квадратичных) выпуклых задач о допустимости (CFP) в евклидовых пространствах Rd. Все вычислительные эксперименты были проведены на обычном ноутбуке Lenovo с процессором Intel(R) Core(TM) i5-4200MQ CPU 1,6 ГГц с 8 Гб оперативной памяти. Программа реализована с помощью пакета MATLAB 2017b. Пример 4.1 (Линейные CFP). В этом примере мы рассматриваем разрешение системы линейных уравнений Ax = b, в частности, восстановление тестового изображения x ∈ [0, 1]n «Ленна» (https://en.wikipedia.org/ wiki/Lenna), состоящего из ограниченного количества томографических проекций. Каждый пиксель обозначим через xi ∈ [0, 1], а каждый элемент из b ∈ Rm, названный томографическим измерением, или единичной проекцией, соответствует интегрированному уровню яркости x вдоль единичного луча. Всякий элемент матрицы aij � 0 задает длину пересечения i-го луча с j-м пикселем. Если луч i и пиксель j не пересекаются, то aij = 0, см. рис. 2. Сложение всех уравнений для всех лучей вместе приводит к системе линейных уравнений Ax = b, а размеры таковы, что θ1 A θ2 A = (AT T θnA ... AT )T , и каждая блочная матрица Aθi соответствует своему проекционному углу. РИС. 2. Геометрическая структура параллельного пучка: множество параллельных лучей пропускается через объект по различным направлениям. Это обычно считается одной проекцией. Случай с двумя проекциями проиллюстрирован на рисунке: каждая проекция отвечает измерению вдоль одного луча и соответствует криволинейному интегралу от кусочно-постоянной функции. В наших экспериментах мы используем в MATLAB шаблон paralleltomo.m из пакета AIR Tools [27], который обеспечивает построение томографической матрицы для заданного вектора углов. Размер шкалы яркости изображения Ленны составляет N = 128 (это означает 128 × 128 пикселей), и выберем для каждого угла число параллельных пучков nA = 100. Еще зададим параметр p = round(128√2) = 169. Выбрав параметры таким образом, получаем переопределенную матрицу A размера (nA ∗ p) × (N 2) = 18100 × 16384. В этом случае нашими данными являются строки A и соответствующие элементы из b. Далее разобьем систему Ax = b на 10 подсистем Ajx = bj размера 1810 × 16384 для j = 1,..., 10. Время прибытия для каждого блока фиксировано и должно составлять 1 минуту. Критерием остановки для всех алгоритмов является 1xk+1 - xk 1 � 10-3, а начальная стартовая точка x0 = 0. Во всех алгоритмах, кроме метода Чиммино, мы полагаем параметры релаксации λk ≡ 1, а в методе Чиммино λk ≡ 1,9. Мы тестируем и сравниваем наши схемы на основе изображения «Ленна» (https://en.wikipedia.org/wiki/Lenna). Напомним, что ортогональная проекция точки x ∈ Rd на гиперплоскость H = {z ∈ Rd | ⊗a, z) � β}, где a ∈ Rd (ненулевой) и β ∈ R, рассчитывается следующим образом: ⎧ ⎨ x при ⊗a, x) = β, PH (x) = ⎩ x - ⊗x, a)-β a 1a12 (4.1) при ⊗a, x) ±= β. На рис. 3 представлено сравнение времени работы программы (в секундах) для метода Чиммино в реальном времени и регулярного метода Чиммино по всем данным. Далее на рис. 4-6 изображены графики (слева), позволяющие сравнить время выполнения программы, а также различия восстановленных изображений, полученных с помощью работы алгоритмов в реальном времени и их регулярных аналогов по всем данным, начиная с момента получения этих данных. На графиках ось y, означающая погрешность, определяется как 1Axk - b1. И хотя разница в восстановленных разными методами изображениях едва заметна невооруженному глазу (изображение слева получено с помощью регулярных алгоритмов, в то время как изображение справа - с помощью аналогов в реальном времени), графики демонстрируют, что для получения необходимой аппроксимации лучше применять алгоритмы в реальном времени (порой даже до получения всех данных). РИС. 3. Метод Чиммино с 10-ю свипами и 10-ю блоками. РИС. 4. Сравнение регулярного метода Качмажа и метода Качмажа в реальном времени. Сверху - график сравнения выполнения регулярного метода и метода Качмажа в реальном времени. Изображение слева получено регулярным методом, а изображение справа - его аналогом в режиме реального времени. Пример 4.2 (Квадратичные CFP). Здесь мы генерируем 10 квадратичных задач о допустимости в R1000, т. е. каждое множество представляет собой шар. В каждом из экспериментов мы увеличиваем количество множеств и сравниваем выполнение всех алгоритмов в режиме реального времени и их регулярных аналогов, ожидающих поступления данных. Каждый шар получается выбором центральной точки ci ∈ R1000 с произвольными, равномерно сгенерированными координатами на отрезке [-5, 5]. Выбор радиуса ri := 1ci1 + αi определяется добавлением к расстоянию от 1 1 центра до начала координат 1ci1 произвольного числа, равномерно выбранного из отрезка [0, 0, 1]. 1 1 Что гарантирует, что шар содержит начало координат и, следовательно, приводит к состоятельной CFP. Векторы x0 задаются произвольным выбором координат из отрезка [-10, 10]. Число ограничений (шаров) варьируется от 200 до 20000. Как и в примере 4.1, для каждой CFP мы разбиваем РИС. 5. Сравнение регулярного метода Качмажа и рандомизированного метода Качмажа в реальном времени. Сверху - график сравнения выполнения этих двух методов. Изображение слева получено регулярным методом, а изображение справа - его аналогом в режиме реального времени. число ограничений (шаров) на 10 блоков и устанавливаем правило остановки: 1xk+1 - xk 1 � ε = 10-7. Параметры релаксации λk полагаем равными 1, а в методе Чиммино они были равны 1,9. В этих экспериментах мы можем наблюдать, что с увеличением количества ограничений возрастает и разница в выполнении схем в реальном времени и их регулярных аналогов. Это вновь подчеркивает потенциальную применимость данных методов к задачам в реальном времени. Результаты, представленные на рис. 7, демонстрируют, что с увеличением количества ограничений (шаров) алгоритмы в реальном времени сходятся быстрее, чем их регулярные аналоги. РИС. 6. Сравнение регулярного метода и «жадного» метода в реальном времени. Сверху - график сравнения выполнения регулярного метода и «жадного» рандомизированного метода в реальном времени. Изображение слева получено регулярным методом, а изображение справа - его аналогом в режиме реального времени. Вспомним, как вычисляется ортогональная проекция точки x ∈ Rn на замкнутый шар B(z, r) := {x ∈ Rn | 1x - z1 � r}, где z ∈ Rn и r > 0: ( x при 1x - z1 � r, PB(z,r)(x) = r z + 1x - z1 (x - z) при 1x - z1 > r. (4.2) РИС. 7. Время выполнения программы (в секундах) для 10 свипов и 10 блоков, соответственно, для метода Чиммино, циклического (Качмаж) и произвольного методов для решения квадратичных задач о допустимости при увеличении количества шаров. 2. ЗАКЛЮЧЕНИЕ В данной работе были рассмотрены задача об общей неподвижной точке (CFPP) и выпуклая задача о допустимости (CFP) в вещественных гильбертовых пространствах. Мы изучили ситуацию, в которой нам с самого начала недоступны целые наборы (операторов/множеств), и мы получаем их постепенно. Руководствуясь недавней работой Ордонеса и соавторов [35], мы представили итерационную схему в режиме реального времени, способную работать с любыми блоками данных (операторов/множеств) для любого конечного числа итераций перед переходом к следующему блоку. Доказательство сходимости нашей схемы основано на недавнем результате Райха и Заласа [37], процедуре модулярного строкового усреднения (MSA). Мы также провели вычислительные эксперименты, демонстрирующие, что схема в реальном времени вычисляет решение быстрее в сравнении с тем случаем, когда все данные известны с самого начала. В то время как Ордонес и соавторы в работе [35] сосредоточены на исследовании только линейных систем уравнений без достаточного математического обоснования, мы проводим более общий анализ задач об общей неподвижной точке с полным теоретическим обоснованием. И хотя структуры методов CARP и DROP не включают в себя алгоритмических структур Tk в алгоритме 3.1, мы планируем продолжить изучение в этом направлении и, более того, получить оценки ошибок и скорости сходимости этой новой схемы.×
Об авторах
А Гибали
Ort Braude College
Email: avivg@braude.ac.il
Д Теллер
Ort Braude College
Email: ktui619@gmail.com
Список литературы
- Губин Л. Г., Поляк Б. Т., Райк Е. В. Метод проекции для нахождения общей точки выпуклых множеств// Журн. выч. мат. и мат. физ. - 1967. - 7.- C. 1-24.
- Aharoni R., Censor Y. Block-iterative projection methods for parallel computation of solutions to convex feasibility problems// Linear Algebra Appl. - 1989. - 120. - C. 165-175.
- Baillon J.-B., Bruck R. E., Reich S. On the asymptotic behavior of nonexpansive mappings and semigroups in banach spaces// Houston J. Math. - 1978. - 4. - C. 1-9.
- Bauschke H. H., Borwein J. On projection algorithms for solving convex feasibility problems// SIAM Rev. - 1996. - 38. - C. 367-426.
- Bauschke H. H., Combettes P. L. Convex analysis and monotone operator theory in Hilbert spaces. - Berlin: Springer, 2011.
- Bauschke H. H., Koch V. R. Projection methods: Swiss army knives for solving feasibility and best approximation problems with halfspaces// В сб.: «Infinite Products of Operators and Their Applications. A Research Workshop of the Israel Science Foundation, Haifa, Israel, May 21-24, 2012». - Providence: Am. Math. Soc., 2015. - С. 1-40.
- Borwein J. M., Tam М. K. A cyclic Douglas-Rachford iteration scheme// J. Optim. Theory Appl. - 2014. - 160.- C. 1-29.
- Browder F. E. Fixed point theorems for noncompact mappings in Hilbert space// Proc. Natl. Acad. Sci. USA. - 1965. - 53. - C. 1272-1276.
- Byrne C. L. A unified treatment of some iterative algorithms in signal processing and image reconstruction// Inverse Problems. - 1999. - 20. - C. 1295-1313.
- Byrne C. L. Applied iterative methods. - Wellsely: AK Peters, 2008.
- Cegielski A. Iterative methods for fixed point problems in Hilbert spaces. - Berlin-Heidelberg: Springer, 2012.
- Cegielski A., Reich S., Zalas R. Regular sequences of quasi-nonexpansive operators and their applications// SIAM J. Optim. - 2018. - 28. - C. 1508-1532.
- Cegielski A., Zalas R. Methods for variational inequality problem over the intersection of fixed point sets of quasi-nonexpansive operators// Numer. Funct. Anal. Optim. - 2013. - 34. - C. 255-283.
- Cegielski A., Zalas R. Properties of a class of approximately shrinking operators and their applications// Fixed Point Theory. - 2014. - 15. - C. 399-426.
- Censor Y., Chen W., Combettes P. L., Davidi R., Herman G. T. On the effectiveness of projection methods for convex feasibility problems with linear inequality constraints// Comput. Optim. Appl. - 2012. - 51.- C. 1065-1088.
- Censor Y., Elfving T., Herman G. T. Averaging strings of sequential iterations for convex feasibility problems// В сб.: «Infinite Products of Operators and Their Applications. A Research Workshop of the Israel Science Foundation, Haifa, Israel, March 13-16, 2000». - Amsterdam: North-Holland, 2001. - С. 101-113.
- Censor Y., Zenios S. A. Parallel optimization: theory, algorithms, and applications. - New York: Oxford Univ. Press, 1997.
- Cimmino G. Calcolo approssiomatto per le soluzioni dei sistemi di equazioni lineari// La Ricerca Sci. XVI. Ser. II. - 1938. - 1. - C. 326-333.
- Combettes P. L. Quasi-Feje´rian analysis of some optimization algorithms// В сб.: «Infinite Products of Operators and Their Applications. A Research Workshop of the Israel Science Foundation, Haifa, Israel, March 13-16, 2000». - Amsterdam: North-Holland, 2001. - С. 115-152.
- Das I., Potra F. A. Subsequent convergence of iterative methods with applications to real-time modelpredictive control// J. Optim. Theory Appl. - 2003. - 119. - C. 37-47.
- Diehl M. Real-Time optimization for large scale nonlinear processes. - Heidelberg: Univ. Heidelberg, 2001.
- Escalante R., Raydan M. Alternating projection methods. - Philadelphia: SIAM, 2011.
- Gala´ ntai A. Projectors and projection methods. - Boston-Dordrecht-London: Kluwer Academic Publ., 2004.
- Goebel K., Reich S. Uniform convexity, hyperbolic geometry, and nonexpansive mappings. - New York- Basel: Marcel Dekker, 1984.
- Gordon D., Gordon R. Component-averaged row projections: A robust block-parallel scheme for sparse linear systems// SIAM J. Sci. Comput. - 2005. - 27. - C. 1092-1117.
- Gordon R., Bender R., Herman G. T. Algebraic reconstruction techniques (art) for three-dimensional electron microscopy and X-ray photography// Bull. Am. Math. Soc. - 1970. - 29. - C. 471-481.
- Hansen P. C., Saxild-Hansen M. AIR Tools - a MATLAB package of algebraic iterative reconstruction methods// J. Comput. Appl. Math. - 2012. - 236, № 8. - C. 2167-2178.
- Iusem A., Jofre´ A., Thompson P. Incremental constraint projection methods for monotone stochastic variational inequalities// arXiv:1703.00272v2. - 2017.
- Kaczmarz S. Angena¨herte auflo¨sung von systemen linearer gleichungen// Bull. Int. l’Acad. Polon. Sci. Lett. A. - 1937. - 35. - C. 355-357.
- Karp R. M. On-line algorithms versus off-line algorithms: How much is it worth to know the future?// В сб.: «Proceedings of the IFIP 12th World Computer Congress on Algorithms, Software, Architecture, Information Processing ’92». - 1992. - 1. - С. 416-429.
- Leventhal L., Lewis A. S. Randomized methods for linear constraints: convergence rates and conditioning// Math. Oper. Res. - 2010. - 35. - C. 641-654.
- Ma˘ rus¸ter S¸ t., Popirlan C. On the Mann-type iteration and the convex feasibility problem// J. Comput. Appl. Math. - 2008. - 212. - C. 390-396.
- Needell D. Randomized Kaczmarz solver for noisy linear systems// BIT Numer. Math. - 2010. - 50.- C. 395-403.
- Opial Z. Weak convergence of the sequence of successive approximations for nonexpansive mappings// Bull. Am. Math. Soc. - 1967. - 73. - C. 591-597.
- Ordon˜ ez C. E., Karonis N., Duffin K., Coutrakon G., Schulte R., Johnson R., Pankuch M. A real-time image reconstruction system for particle treatment planning using proton computed tomography (PCT)// Phys. Proc. - 2017. - 90. - C. 193-199.
- Penfold S., Censor Y., Schulte R. W., Bashkirov V., McAllister S., Schubert K. E., Rosenfeld A. B. Blockiterative and string-averaging projection algorithms in proton computed tomography image reconstruction// В сб.: «Biomedical Mathematics: Promising Directions in Imaging, Therapy Planning and Inverse Problems». - Madison: Medical Phys. Publ., 2010. - С. 347-368.
- Reich S., Zalas R. A modular string averaging procedure for solving the common fixed point problem for quasi-nonexpansive mappings in hilbert space// Numer. Algorithms. - 2016. - 72. - C. 297-323.