Журнал вычислительной математики и математической физики, 2021, T. 61, № 5, стр. 723-743

2.1. Преобразование Шэнкса

Запись уравнения (2) для индексов $n,\; \ldots ,\;n + k$ позволяет определить коэффициенты ${{a}_{i}}$ как решение системы линейных уравнений, а затем вычислить $t$ при помощи (2). Конечно, если последовательность $({{S}_{n}})$ не принадлежит ядру, все еще возможно записать и решить систему уравнений, определяющую коэффициенты ${{a}_{i}}$ и далее из уравнений (2) определить так называемое преобразование Шэнкса, обычно обозначаемое ${{e}_{k}}({{S}_{n}})$, так как результат будет зависеть и от $n$, и от $k$. В результате получим

где $\Delta $ – разностный оператор дифференцирования “справа”, определяемый $\Delta {{S}_{i}} = {{S}_{{i + 1}}} - {{S}_{i}}$ для всех $i$.

Решая разностные уравнения (1) через запись характеристического многочлена и поиск всех его корней, мы получаем следующий результат [11] (см. также [12] с более детальными комментариями).

Теорема 1. Необходимым и достаточным условием для ${{e}_{k}}({{S}_{n}}) = S$ при всех $n$ является условие, что последовательность $({{S}_{n}})$ удовлетворяет

где ${{r}_{i}} \ne 0$ для $i = 1,\; \ldots ,\;p$, ${{A}_{i}}$, ${{B}_{i}}$ и ${{C}_{i}}$ – многочлены порядка $n$ такие, что если ${{d}_{i}}$ равняется степени ${{A}_{i}}$ плюс $1$ для всех $i = 1,\; \ldots ,\;p$ и максимуму из степеней ${{B}_{i}}$ и ${{C}_{i}}$ плюс $1$ для $i = p = 1, \ldots ,q$, тос обязательным условием, что $m = - 1$, если нет члена ${{\delta }_{{in}}}$ (символ Кронекера).

Коэффициенты полиномов ${{A}_{i}}$, ${{B}_{i}}$ и ${{C}_{i}}$ определяются начальными условиями в уравнении (1). Из этого результата видно, что многие итерационные методы, используемые в численном анализе и прикладной математике, порождают последовательности такого вида или близкие к нему, что объясняет важность метода экстраполяции Шэнкса.

2.2. $\varepsilon $-Алгоритм

Несколько преобразований последовательности могут быть выражены как отношение двух определителей. Таким образом, для нахождения преобразованных членов исходной последовательности должны быть получены алгоритмы, избегающие численного вычисления этих определителей. Рекурсивная процедура, $\varepsilon $-алгоритм, для реализации преобразования Шэнкса была предложена в 1956 г. Питером Винном [13]. Организация этой процедуры очень проста:

при $\varepsilon _{{ - 1}}^{{(n)}} = 0$ и $\varepsilon _{0}^{{(n)}} = {{S}_{n}}$ для $n = 0,\;1,\; \ldots $, и это выполнено для всех $k$ и $n$, Таким образом, значения $\varepsilon _{{2k + 1}}^{{(n)}}$ – это промежуточные результаты и не представляют интереса для наших целей.

Значения $\varepsilon _{k}^{{(n)}}$ обычно изображаются в виде двумерного $\varepsilon $-массива, где нижний индекс $k$ сохраняется вдоль столбцов, а верхний индекс $n$ сохраняется вдоль нисходящей диагонали:

В процедуре (4) $\varepsilon $-алгоритма обрабатываются значения, расположенные в четырех вершинах следующего ромба:

В случае, когда в знаменателе (4) два значения оказались равными, вычислительная процедура должна быть остановлена (данное явление принято называть обвалом). Но если значения в знаменателе почти равны, возникает неустойчивость при вычислении очередного члена (называемая почти–обвалом). Для разрешения этих недостатков алгоритма Винн, см. [14], также предложил конкретные правила для скалярного $\varepsilon $-алгоритма, которые позволяют повысить устойчивость вычислений.

2.3. Аппроксимации Паде

$\varepsilon $-Алгоритм, описанный в п. 2.2, связан с известными аппроксимациями Паде, используемыми, в частности, для суммирования сходящихся и расходящихся рядов [6], [15].

Пусть $f(z) = \sum\nolimits_{i = 0}^\infty {{{c}_{i}}} {{z}^{i}}$ – формальный степенной ряд. Аппроксимация Паде для $f$ – это рациональная функция, разложение в степенной ряд которой по возрастающим степеням $z$ является близким к ряду $f$ настолько, насколько это возможно, т.е. с точностью до суммы $p + q$ включительно, где $p$ – степень его числителя, а $q$ – степень $q$ его знаменателя.

Обозначим такую аппроксимацию через

где коэффициенты можно получить из следующих выражений: при условии, что ${{c}_{i}} = 0$ для $i < 0$.

Линейная система, порожденная последними $q$ уравнениями, содержит $q + 1$ неизвестную ${{b}_{0}},\; \ldots ,\;{{b}_{q}}$ и, следовательно, существует ее нетривиальное решение. Положив ${{b}_{0}} = 1$, мы можем получить значения ${{b}_{1}},\; \ldots ,\;{{b}_{q}}$, решая линейную систему, предполагая ее невырожденность. Зная значения ${{b}_{i}}$, из первых $p + 1$ уравнений мы можем вычислить ${{a}_{i}}$.

Выполнено

при ${{S}_{n}}(z) = \sum\nolimits_{i = 0}^n {{{c}_{i}}} {{z}^{i}}$.

Применяя $\varepsilon $-алгоритм к частичным суммам $f$ при $\varepsilon _{0}^{{(n)}} = {{S}_{n}}(z)$, можно заметить, сопоставив предыдущую формулу с (3), что

Аппроксиманты со степенью знаменателя большей, чем степень числителя, могут быть получены через применение $\varepsilon $-алгоритма к частичным суммам обратного ряда $g$ функции $f$, формально определяемого через $f(z)g(z) = 1$ и через свойство ${{[p{\text{/}}q]}_{f}}(z){{[q{\text{/}}p]}_{g}}(z) = 1$. Ряд $g$ существует только тогда, когда ${{c}_{0}} \ne 0$.

3.1. Случай векторных и матричных последовательностей

В работе [16] Винн предложил расширение своего $\varepsilon $-алгоритма для случая векторных и матричных последовательностей. В данном случае правила (4) остаются прежними, а степень $ - 1$ соответствует получению обратной матрицы. Ядро преобразования снова можно записать в форме (1). Позже А. Салам (см. [17]) предложил доказательство, что данный результат все еще верен, если коэффициенты ${{a}_{i}}$ являются квадратными матрицами, и произведение возникает либо слева, либо справа. Алгоритм может быть обобщен и для случая прямоугольных матриц через замену операции обращения матриц вычислением псевдообратных, но результат, касающийся ядра, не остается верным. Подробнее случай $\varepsilon $-алгоритма Винна для матричных последовательностей рассмотрен в [18].

В случае векторных последовательностей Винн определил обращение вектора $u$ через ${{u}^{{ - 1}}} = u{\text{/}}(u,u)$. Векторы $\varepsilon _{{2k}}^{{(n)}}$ тоже могут быть выражены, как и в скалярном случае, через отношение определителей, но размерности $2k + 1$ вместо исходной $k + 1$, как было доказано П.Р. Грейвс-Моррисом и К. Дженкинсом [9], [19]. Применив алгебру Клиффорда, Дж.Б. МакЛеод доказал, что ядро векторного $\varepsilon $-алгоритма также содержит векторы формы (1) [20], где коэффициенты ${{a}_{i}}$ являются вещественными числами. Позже А. Салам доказал, что ядром векторного $\varepsilon $-алгоритма является непосредственно (1), но с коэффициентами ${{a}_{i}}$, принадлежащими алгебре Клиффорда, ассоциированной с ${{\mathbb{R}}^{p}}$ [21], [22], где $p$ – размерность векторов последовательности ${{S}_{n}}$. Аналогичная (3) формула по-прежнему имеет место, но определители необходимо заменить на конструкторы (designants), которые обобщают их в некоммутативной алгебре [17], [23].

3.2. Топологическое преобразование Шэнкса

Столкнувшись с трудностями, возникающими в алгебраических теориях для векторных и матричных $\varepsilon $-алгоритмов, К. Брежинский решил заново заняться этой задачей с начальных этапов [24]. Пусть $({{S}_{n}})$ – последовательность элементов топологического векторного пространства $E$ над $\mathbb{C}$, или просто над $\mathbb{R}$ (топология необходима, чтобы уметь рассматривать вопросы сходимости). Задачей было построить преобразование последовательности $({{S}_{n}}) \mapsto \left\{ {({{e}_{k}}({{S}_{n}})} \right\}$ такое, что для всех последовательностей $({{S}_{n}})$, удовлетворяющих (1) с $t = S$, при фиксированном $k$ и с коэффициентами ${{a}_{i}} \in \mathbb{C}$, ${{e}_{k}}({{S}_{n}}) = {{a}_{0}}{{S}_{n}} + {{a}_{1}}{{S}_{{n + 1}}} + \; \ldots \; + {{a}_{k}}{{S}_{{n + k}}} = S \in E$ для всех $n$. Записав (1) для индексов $n$ и $n + 1$ и вычислив разность, получим

Следующей задачей стало вычислить коэффициенты ${{a}_{i}}$. Для этого полученное выражение было преобразовано в выражение над $\mathbb{C}$. Пусть $y$ является элементом алгебраического сопряженного пространства $E{\text{*}}$ для $E$, являющегося пространством линейных функционалов над $E$, и обозначим через $\left\langle { \cdot , \cdot } \right\rangle $ соотношение двойственности между $E{\text{*}}$ и $E$. Применение этого соотношения к выражениям, описанным выше, позволяет получить Записав это выражение для индексов от $n$ до $n + k - 1$, добавив условие нормировки и решив соответствующую линейную систему для ${{a}_{i}}$, мы получим (первое) топологическое преобразование Шэнкса, определяемое через:

Второе топологическое преобразование Шэнкса, обозначаемое через ${{\tilde {e}}_{k}}({{S}_{n}})$, можно определить, сделав замену первой строки числителя на ${{S}_{{n + k}}},\; \ldots ,\;{{S}_{{n + 2k}}}$. Если $E = \mathbb{C}$, два данных преобразования редуцируются до скалярного преобразования Шэнкса (3).

3.3. Топологические $\varepsilon $-алгоритмы

Следующей задачей стало получение рекурсивного алгоритма для реализации полученного преобразования. Этот алгоритм также должен был быть сокращен до (4) в скалярном случае. Ключевым моментом стало определение обращения, которое ограничивается обычным преобразованием для скаляров и векторным в форме, использованной Винном в его векторном $\varepsilon $-алгоритме. Решением оказалось определение обращения пары $(y,u) \in E{\text{*}} \times E$ через ${{y}^{{ - 1}}} = {u \mathord{\left/ {\vphantom {u {\left\langle {y,u} \right\rangle }}} \right. \kern-0em} {\left\langle {y,u} \right\rangle }} \in E$ и ${{u}^{{ - 1}}} = {y \mathord{\left/ {\vphantom {y {\left\langle {y,u} \right\rangle }}} \right. \kern-0em} {\left\langle {y,u} \right\rangle }} \in E{\text{*}}$. Такое обращение удовлетворяет всем свойствам, которым должна удовлетворять данная операция: обратный элемент может быть вычислен, а отношение двойственности элемента со своим обратным равно единице. Это определение обращения привело к тому, что правило этого, так называемого, топологического $\varepsilon $-алгоритма пришлось разделить на два правила: одно для членов с еще четным нижним индексом, принадлежащих $E$, и второе для элементов с нечетным нижним индексом, которые находятся в $E{\text{*}}$ и являются промежуточными результатами, как в скалярном случае. Эти наблюдения приводят к (первому) топологическому $\varepsilon $-алгоритму (TEA1). Второй алгоритм (TEA2) позволяет реализовать рекурсивно второе преобразование. Доказано, что $\hat {\varepsilon }_{{2k}}^{{(n)}} = {{\hat {\varepsilon }}_{k}}({{S}_{n}})$, как определено уравнением (6), а для второго преобразования и алгоритма, что $\tilde {\varepsilon }_{{2k}}^{{(n)}} = {{\tilde {e}}_{k}}({{S}_{n}})$ [24].

Правила алгоритмов TEA были довольно сложными (два разных правила, они требуют хранения элементов $E$ и $E{\text{*}}$, отношение двойственности было рекурсивно использовано в правилах). Таким образом, было трудно реализовать их для линейных функционалов общего вида, и в прошлом эти алгоритмы использовались только в простом случае, когда $E$ совпадает с ${{\mathbb{R}}^{p}}$ (или ${{\mathbb{C}}^{p}}$), $y$ был вектором, и отношение двойственности сводится к внутреннему произведению между действительными или комплексными векторами. В 2014 г. авторы этой работы осознали, спустя 40 лет после публикации [24], что соотношение, записанное ранее, позволяет существенно упростить правила двух топологических $\varepsilon $-алгоритмов [25], что привело к формулировке упрощенных топологических $\varepsilon $-алгоритмов (STEA). Каждый из двух этих упрощенных алгоритмов может быть записан в 4 эквивалентных формах. Приведем только правило второго (STEA2), которое еще более эффективно и включает работу только с величинами с четными нижними индексами:

при $\tilde {\varepsilon }_{0}^{{(n)}} = {{S}_{n}}$, и где значения $\varepsilon _{k}^{{(n)}}$ вычисляются с использованием скалярного $\varepsilon $-алгоритма Винна, примененного к последовательности $(\varepsilon _{0}^{{(n)}} = \left\langle {y,{{S}_{n}}} \right\rangle )$.

По сравнению с оригинальными топологическими $\varepsilon $-алгоритмами упрощенные используют только одно треугольное правило вместо двух более сложных, также они позволяют использовать только элементы $E$ с нижними четными индексами, отсутствует необходимость вычислять или хранить элементы $E{\text{*}}$, соотношение двойственности используется только на этапе инициализации скалярного $\varepsilon $-алгоритма, требования по количеству используемых ячеек памяти были существенно снижены (что важно в случае работы с векторными и матричными последовательностями), а устойчивость алгоритмов была повышена с использованием конкретных правил Винна для скалярного $\varepsilon $-алгоритма [14]. Больше деталей и информацию о программной реализации можно найти в работе [26].

3.4. Другие преобразования Шэнкса

Как мы видим теперь, существуют различные способы использования уравнения (5) для вычисления коэффициентов ${{a}_{i}}$ в преобразовании. Всегда предполагается, что сумма коэффициентов ${{a}_{i}}$ равна 1.

Вместо использования одного элемента $y \in E{\text{*}}$ при меняющемся $n$ мы можем зафиксировать $n$ и выбрать $k$ линейно-независимых функционалов ${{y}_{i}} \in E{\text{*}}$. Эта процедура приводит к методу минимальной полиномиальной экстраполяции (MMPE) [24], который записывается через:

Для этого преобразования мы используем те же обозначения ${{e}_{k}}({{S}_{n}})$, что и для топологического преобразования Шэнкса, даже если оно отличается от него. Это преобразование может быть реализовано рекурсивно через $S\beta $-алгоритм согласно Л. Джейбилоу [27] (см. также [28]).

Метод полиномиальной экстраполяции (MPE) [29] и метод экстраполяции с редуцированным рангом (RRE) [30], [31] можно напрямую получить из выражения MMPE, подставив ${{y}_{i}} = \Delta {{S}_{{n + i - 1}}}$ для первого и ${{y}_{i}} = {{\Delta }^{2}}{{S}_{{n + i - 1}}}$ для второго. Однако оба этих метода применимы лишь в случае векторных последовательностей.

Существует и несколько других обобщений, упомянем, например, [18], [32] и книги [1], [4], [5].

4.1. Метод ускорения (AM)

Пусть $({{S}_{n}})$ – последовательность скаляров, векторов, матриц или тензоров, сходимость которой необходимо ускорить с помощью одного из $\varepsilon $-алгоритмов. Предположим, что для фиксированного значения $k$ мы хотим вычислить последовательность $(\varepsilon _{{2k}}^{{(n)}})$. Члены последовательности $({{S}_{n}})$ определяются один за другим в соответствующей программной реализации, и четные столбцы $\varepsilon $-массива вычислены настолько далеко, насколько это возможно. Начиная с ${{S}_{0}}$ и ${{S}_{1}}$, ни один член не может быть получен. Тогда, вводя ${{S}_{2}}$, можно вычислить $\varepsilon _{2}^{{(0)}}$. С использованием дополнительно еще ${{S}_{3}}$ мы можем получить $\varepsilon _{2}^{{(1)}}$. Далее с ${{S}_{4}}$, можно получить $\varepsilon _{4}^{{(0)}}$. Таким образом, нисходящая “лестница” вдоль главной диагонали $\varepsilon $-массива получается до тех пор, пока, введя член ${{S}_{{2k}}}$, мы не сможем вычислить $\varepsilon _{{2k}}^{{(0)}}$. Далее каждый новый член ${{S}_{{2k + 1}}},{{S}_{{2k + 2}}},\; \ldots $ позволяет вычислить следующий член в столбце $2k$, т.е. $\varepsilon _{{2k}}^{{(1)}},\varepsilon _{{2k}}^{{(2)}},\; \ldots $. Если вместо вычисления при фиксированном значении $k$ последовательности $(\varepsilon _{{2k}}^{{(n)}})$ требуется вычислить элементы главной нисходящей диагонали или “лестницу” вдоль нее, достаточно дать большое значение $k$ и прекратить вводить новые члены последовательности для преобразования. Данная стратегия называется методом ускорения (AM). Фигура 1 демонстрирует сопутствующий алгоритм, а фиг. 2 приводит пример данной стратегии.

Поясним далее, как построить таблицу, см. фиг. 2, по последовательности $({{S}_{n}})$. Простейший способ построения $\varepsilon $-массива – это вычислить столбцы один за другим, начиная с первого столбца, используя далее для всех столбцов (кроме первого) значения из двух предыдущих. С точки зрения требований по количеству используемых ячеек памяти и количества действий данная процедура может быть слишком затратной в случае, если членами последовательности являются векторы или матрицы. Поэтому обычно при реализации алгоритмов используется прием бегущего окна, который известен благодаря Винну [14], использовавшему данную технику в своем скалярном $\varepsilon $-алгоритме.

Конечно, эта техника может использоваться и в других алгоритмах, в частности, в топологическом $\varepsilon $-алгоритме или в его упрощенных версиях. Идея состоит в том, чтобы вести работу от возрастающих диагоналей, т.е. мы вычисляем каждую новую возрастающую диагональ массива, используя предыдущую диагональ (или две предшествующие в случае первого топологического $\varepsilon $-алгоритма). Большее число деталей сообщено в работе [26] и в пользовательской инструкции соответствующего программного обеспечения [33].

4.2. Метод рестарта (RM)

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

и затем ускорить его сходимость с помощью AM. Однако можно использовать и другую стратегию стеффенсенского типа, известную под названием метода рестарта (RM). Начиная с заданного ${{x}_{0}}$, мы вводим ${{u}_{0}} = {{x}_{0}}$ и проводим некоторое количество итераций Пикара, базовых итераций, ${{u}_{{i + 1}}} = F({{u}_{i}})$ для $i = 0,\; \ldots ,\;2k - 1$ для фиксированного значения $k$. Далее применение $\varepsilon $-алгоритма к этим $2k + 1$ элементам генерирует ${{x}_{1}} = \varepsilon _{{2k}}^{{(0)}}$. Базовые итерации после этого можно начать уже с ${{u}_{0}} = {{x}_{1}}$ и продолжать далее, что порождает последовательность $({{x}_{n}})$, которая при некоторых предположениях сходится к неподвижной точке $x$. Фигура 3 демонстрирует соответствующий алгоритм. В векторном случае $F:{{\mathbb{R}}^{p}} \mapsto {{\mathbb{R}}^{p}}$, Х. Ле Ферран [34] доказал, что при выборе $k = p$ использование любого из топологических $\varepsilon $-алгоритмов в RM генерирует последовательность, которая при некоторых предположениях сходится квадратично к неподвижной точке $x \in {{\mathbb{R}}^{p}}$ для отображения $F$. Данная процедура называется обобщенным методом Стеффенсена (GSM), так как при $p = 1$ процедура оказывается непосредственно методом Стеффенсена [35]. Для метода RRE подробный анализ представлен в работе Сиди [36].

4.3. Программные реализации

Несмотря на то что во многих публикациях в литературе утверждается, что преобразования последовательностей и соответствующие алгоритмы могут быть полезными методами для получения или улучшения численных решений широкого круга задач численного анализа и прикладной математики, потенциальных пользователей часто отговаривают использовать их из-за сложности их реализации простым и оптимизированным способом. По этой причине в 1991 г. в дополнение к [1] на языке $FORTRAN$ 77 была опубликована библиотека (процедуры и демонстрационные программы). Недавно библиотека $EPS{\text{fun}}$ [33] на языке ${\text{M}}ATLAB$ была внедрена в библиотеку с открытым программным кодом Netlib. Она содержит необходимые функции по реализации $\varepsilon $-алгоритмов Винна, оригинальных топологических $\varepsilon $-алгоритмов и их упрощенных версий. Напомним, что оригинальные и упрощенные топологические $\varepsilon $-алгоритмы требуют определения элемента $y \in E{\text{*}}$. В оригинальных алгоритмах $y$ непосредственно используется в одном из правил, в то время как в упрощенных данный элемент играет роль только на этапе инициализации скалярного $\varepsilon $-алгоритма, который применяется к последовательности $\left( {\left\langle {y,{{S}_{n}}} \right\rangle } \right)$. Демонстрационные программы, шаблоны проектов (как для методов AM, так и для RM), а также пользовательские инструкции являются частью данного программного обеспечения.

5.1. Скалярный $\varepsilon $-алгоритм

Начнем с демонстрации двух случаев применения скалярного $\varepsilon $-алгоритма. В первом случае рассматривается задача об ускорении сходимости конкретной последовательности. Во втором случае алгоритм применяется к рядам Фурье.

5.1.1. Скалярные последовательности. Применим скалярный $\varepsilon $-алгоритм к частичным суммам ${{f}_{{2k}}}$ ряда

который сходится для $\left| z \right| \leqslant 1$, кроме $z = - 1$.

Для $z = 1$ имеем $ln2 = 0.6931471805599453\; \ldots $. Для $z = 2$ ряд расходится, и мы имеем $ln3 = 1.098612288668110\; \ldots $. Аппроксимации Паде приводят к результатам, содержащимся в табл. 1 (напомним, что $\varepsilon _{{2k}}^{{(0)}} = {{[k{\text{/}}k]}_{f}}(z)$).

Мы можем видеть, что $\varepsilon $-алгоритм позволяет получить суммы для расходящихся последовательностей (рядов, в данном конкретном случае). Это следует из факта, что значения $\varepsilon _{{2k}}^{{(n)}}$ связаны с аппроксимациями Паде и что эти аппроксимации позволяют получить аналитическое продолжение для некоторых рядов вне области сходимости, что и демонстрирует приведенный пример.

5.1.2. Эффект Гиббса. Далее мы приведем численные примеры, показывающие, что $\varepsilon $-алгоритм позволяет локализовать разрывы в рядах Фурье, ускорить их сходимость и уменьшить выбросы, возникающие из-за эффекта Гиббса.

Рассмотрим ряд Фурье общего вида (без ограничений общности, предположим, что ряд является периодическим на отрезке $[0,2\pi ]$ и обладает нулевым постоянным членом):

Складывая этот ряд с сопряженным:

как с мнимой частью, мы получаем комплексный ряд Фурье:

Применим теперь $\varepsilon $-алгоритм к частичным суммам $F(t)$ и далее возьмем вещественную часть значений $\varepsilon _{{2k}}^{{(n)}}$. Данная процедура известна под названием комплексного $\varepsilon $-алгоритма, а получаемые значения обозначаются $_{C}\varepsilon _{{2k}}^{{(n)}}$. В случае, если скалярный $\varepsilon $-алгоритм применяется к частичным суммам $S(t)$ (как в п. 2.3), мы считаем, что речь идет о вещественном $\varepsilon $-алгоритме и обозначаем получаемые значения $_{R}\varepsilon _{{2k}}^{{(n)}}$.

В качестве примера рассмотрим ряд Фурье на отрезке $[0,\pi ]$:

На фиг. 4а–в сплошная линия соответствует точному результату, пунктирная представляет результаты использования $\varepsilon $-алгоритма. На фиг. 4г мы демонстрируем уровень ошибки ${{S}_{9}}(t)$ (что соответствует частичной сумме из первых 9 членов $S(t)$) относительно точных значений при значении $a = 0.98$. График $_{R}\varepsilon _{4}^{{(0)}}$ демонстрирует, что алгоритм позволяет выявить почти все разрывы.

График значений ошибки демонстрирует улучшения, полученные в результате использования комплексного $\varepsilon $-алгоритма. Эта идея была предложена Винном [37], затем приведена в работе [38], где $\varepsilon $-алгоритм использовался для локализации разрывов, и окончательно оформлена в [39].

5.2. Работа с системами уравнений

Применим теперь ранее описанные алгоритмы к задачам о решении систем линейных и нелинейных уравнений.

5.2.1. Метод Качмарца. Метод Качмарца – итерационный метод решения систем линейных уравнений [40]. Хотя история этого алгоритма началась в 1937 г., в последнее время к нему возрождается интерес из-за его хорошей структуры для параллельной реализации. Другим достоинством этого подхода является его сходимость для всех случаев, как было доказано [41]. Тем не менее сходимость обычно оказывается медленной, и несколько процедур для ее ускорения можно найти в литературе.

Одна итерация циклического метода Качмарца для решения линейной системы $Ax = b$ размерности $N$ состоит в следующих действиях:

Обозначим через ${{a}_{i}} = {{A}^{{\text{т}}}}{{e}_{i}}$ вектор-столбец, сформированный из $i$-го столбца $A$, а через ${{b}_{i}}$ – компоненту с номером $i$ вектора правой части $b$. Тогда вычисление каждого из векторов ${{p}_{i}}$ в шаге итерации не требует вычисления умножения матрицы на вектор, что приводит к:

Данное замечание позволяет легко получить параллельную реализацию алгоритма.

Обозначим далее ${{P}_{i}} = I - A{{\alpha }_{i}}e_{i}^{{\text{т}}}$ и ${{Q}_{i}} = {{A}^{{ - 1}}}{{P}_{i}}A$ при ${{\alpha }_{i}} = {{{{A}^{{\text{т}}}}{{e}_{i}}} \mathord{\left/ {\vphantom {{{{A}^{{\text{т}}}}{{e}_{i}}} {{{{\left\| {{{A}^{{\text{т}}}}{{e}_{i}}} \right\|}}^{2}}}}} \right. \kern-0em} {{{{\left\| {{{A}^{{\text{т}}}}{{e}_{i}}} \right\|}}^{2}}}}$. Выполнено ${{x}_{{n + 1}}} - x = Q({{x}_{n}} - x)$ при $Q = {{Q}_{N}}\; \cdots \;{{Q}_{1}}$, что есть ${{x}_{n}} - x = {{Q}^{n}}({{x}_{0}} - x)$. Данное выражение демонстрирует, что итерационная последовательность $({{x}_{n}})$ принадлежит к ядру топологического преобразования Шэнкса. Так, в теории его применение может привести к получению точного решения системы. Однако, так как $N$ обычно (очень) велико, преобразование нельзя использовать на практике, хотя его можно использовать в целях ускорения сходимости, как это делается в работе [42].

Мы рассматриваем $parter$ матрицу $A$, $N = 1000$, $\varkappa (A) \simeq 4.2306$, – тёплицеву матрицу с сингулярными числами, близкими к $\pi $. Фигура 5 демонстрирует сравнение метода Качмарца с MMPE, MPE, RRE, топологическим $\varepsilon $-алгоритмом и с векторным $\varepsilon $-алгоритмом при $k = 5$.

Мы видим, что все методы достигают хорошей точности с преимуществом векторного $\varepsilon $-алгоритма. Более того, его сходимость оказывается более гладкой. В приведенном примере старшее собственное значение матрицы $A$ равно $0.8732178$, а следующее за ним – $0.3170877$. Таким образом, согласно теоретическим результатам, общим для всех этих методов, хорошее ускорение может наблюдаться даже при $k = 1$.

5.2.2. Интегральные уравнения. Рассмотрим следующее нелинейное интегральное уравнение Фредгольма II рода с ядром $K$:

Будем считать, что выполнены все обычные условия, гарантирующие существование единственного решения.

Стандартным методом решения этого уравнения является приближение интегрального оператора с помощью квадратурных формул:

где $x_{0}^{{(p)}},\; \ldots ,\;x_{p}^{{(p)}}$ – это $p + 1$ точка на отрезке $[a,b]$, а верхний индекс $p$ соответствует зависимости от номера выбранной точки. Напомним, что веса $w_{j}^{{(p)}}$ строго положительны, а их сумма равна $b - a$. Таким образом, уравнение аппроксимируется через

Далее мы приближаем решение ${{u}_{p}}$ методом коллокаций для точек $t_{i}^{{(p)}} = x_{i}^{{(p)}}$ при $i = 0,\; \ldots ,\;p$. При фиксированном значении $p$ обозначим для удобства ${{t}_{i}} = {{x}_{i}} = t_{i}^{{(p)}} = x_{i}^{{(p)}}$, ${{f}_{i}} = f({{t}_{i}})$, $w_{i}^{{(p)}} = {{w}_{i}}$ и определим далее аппроксимации ${{u}_{i}}$ для ${{u}_{p}}({{t}_{i}})$, $i = 0,\; \ldots ,\;p$, как решение системы из $p + 1$ нелинейных уравнений:

Для решения этой системы для начала мы используем итерационный метод Пикара:

или, в более общем случае схему релаксации: где $\alpha $ – параметр, который необходимо найти, а $u_{i}^{{(0)}}$, $i = 0,\; \ldots ,\;p$ – начальное приближение решения в точках ${{t}_{i}}$.

Затем эти итерации вводятся в один из методов, полученных из преобразования Шэнкса, таких как MPE, MMPE, RRE, векторный или упрощенный топологический $\varepsilon $-алгоритмы. Полный математический анализ этой методологии и алгоритм описаны в работе [43]. Программы на языке ${\text{M}}ATLAB$ с несколькими примерами находятся в свободном доступе на сайте Matlab File Exchange.

В качестве примера мы рассмотрим интегральное уравнение

чьим решением является $u(x) = \sqrt x $. Из преобразований Шэнкса здесь используются первый и второй упрощенные топологические $\varepsilon $-алгоритмы, каждый с использованием AM и RM стратегий.

Результаты использования стратегии AM при параметрах $\alpha = 0.2$, $p = 31$, $y = {{(1,\; \ldots ,\;1)}^{{\text{т}}}}$ и $2k = 10$ представлены на фиг. 6а, результаты при $2k = 4$ для стратегии RM представлены на фиг. 6б.

Нашей целью было не конкурировать со сложными методами, которые можно найти в литературе, а показать, что простая процедура, сопровождаемая методом ускорения, может быть весьма эффективной.

5.2.3. Нелинейные алгебраические уравнения. Для ускорения сходимости итерационного метода ${{x}_{{n + 1}}} = {{x}_{n}} + \alpha f({{x}_{n}})$ при решении системы линейных или нелинейных уравнений $f(x) = 0$ мы можем воспользоваться методом ускорения (AM) или методом рестарта (RM).

Рассмотрим далее нелинейную систему:

с решением ${{( - 1,\;1,\;0)}^{{\text{т}}}}$.

Начиная с точки ${{x}_{0}} = 0 \in {{\mathbb{R}}^{3}}$, мы наблюдаем результаты, представленные на фиг. 7а для AM при $\alpha = 0.2$ и $\varepsilon _{4}^{{(n)}}$. При $\alpha = 0.1$ наблюдаются аналогичные результаты. Для GSM ($2k = 6$) при $\alpha = 0.1$ использование STEA2 позволяет получить результаты лучше, чем STEA1, что показано на фиг. 7б. Ступенчатая структура графика на фиг. 7б наглядно демонстрирует квадратичную сходимость.

Квадратичная сходимость GSM также наблюдается для MMPE, MPE и RRE, что доказано в работе [44].

5.2.4. Матричные уравнения. Рассмотрим симметричное матричное уравнение Cтейна, известное также как уравнение Ляпунова с дискретным временем:

$F \in {{\mathbb{R}}^{{m \times s}}}$, $s \ll m$, при собственных значениях $A$, находящихся внутри единичного круга.

Известен итерационный метод, позволяющий получить численное решение (см. [25]):

$n = 0,\;1,\; \ldots $, при ${{S}_{0}} = 0$.

Результаты на фиг. 8 соответствуют матрице $moler$ размерности 500 для матрицы $A$, нормированной на скалярное значение такое, что ее спектральный радиус стал равен 0.9. Матрица $moler$ – симметричная положительно-определенная матрица с одним малым собственным значением. Ее элементы определяются как ${{a}_{{ij}}} = min(i,j) - 2$ и ${{a}_{{ii}}} = i$. Матрица $F$ – $parter$ матрица размерности $500 \times 30$. Данные две матрицы принадлежат набору примеров $gallery$. При $2k = 6$ для упрощенных топологических $\varepsilon $-алгоритмов, где дуальное произведение с функционалом $y$ соответствует следу матрицы, мы применяем стратегию AM.

Сплошная линия на фиг. 8 соответствует итерационному методу, штриховая – алгоритму STEA1, а штрихпунктирная – методу STEA2.

5.3. Матричные функции

Рассмотрим задачу вычисления:

где $A$ – вещественная матрица.

Данный ряд сходится, если спектральный радиус матрицы $A$ строго меньше единицы.

В данном примере мы покажем, что упрощенные топологические $\varepsilon $-алгоритмы позволяют ускорить сходимость частичных сумм данного ряда, а также получить сходимость к приближенному значению в случае расходимости исходного ряда.

Выберем квадратную случайную матрицу $B$ размерности $50$, значение $r$ и определим

Таким образом, спектральный радиус $A$ равен $r$.

Используя упрощенные топологические $\varepsilon $-алгоритмы при $\left\langle {y, \cdot } \right\rangle = \operatorname{tr} ({\kern 1pt} \cdot {\kern 1pt} )$, при $\varepsilon _{8}^{{(n)}}$ мы получаем результаты для частичных сумм ряда, представленные на фиг. 9, соответствующие $r = 0.9$ (случай сходимости, на фиг. 9а) и $r = 1.2$ (случай расходимости, на фиг. 9б). На графиках изображены значения евклидовой нормы погрешности с использованием функции $logm$ для языка программирования ${\text{M}}ATLAB$ (которая позволяет вычислить основной матричный логарифм).

6.1. Тензорные ${{\ell }^{p}}$-собственные пары

Так как $T$ переводит ${{\mathbb{R}}^{n}}$ в себя, мы можем ввести концепцию собственных пар для ${\mathbf{T}}$ через отображение $T$. Вещественное число $\lambda \in \mathbb{R}$ будем считать ${{\ell }^{p}}$-собственным значением ${\mathbf{T}}$ [45] (при $p > 1$), соответствующим ${{\ell }^{p}}$-собственному вектору ${\mathbf{x}} \in {{\mathbb{R}}^{n}}$, если

где ${{\left\| {\mathbf{x}} \right\|}_{p}}$ соответствует обычно ${{\ell }^{p}}$-норме: и отображение ${{\Phi }_{p}}\,:{{\mathbb{R}}^{n}} \to {{\mathbb{R}}^{n}}$ определяется поэлементно ${{\Phi }_{p}}{{({\mathbf{x}})}_{i}} = {{\left| {{{x}_{i}}} \right|}^{{p - 2}}}{{x}_{i}} = {\text{sign}}({{x}_{i}}){{\left| {{{x}_{i}}} \right|}^{{p - 1}}}$ при $i = 1,\; \ldots ,\;n$. Если $p \ne d$, то ${{\ell }^{p}}$-собственные векторы нельзя определить с точностью до скалярного множителя. По этой причине мы добавляем условие нормировки ${{\left\| {\mathbf{x}} \right\|}_{p}} = 1$.

Тензор ${\mathbf{T}} = ({{t}_{{{{i}_{1}}, \ldots ,{{i}_{d}}}}})$ будем называть неотрицательным и писать ${\mathbf{T}} \geqslant 0$, если ${{t}_{{{{i}_{1}}, \ldots ,{{i}_{d}}}}}0$ для всех ${{i}_{j}} = 1,\; \ldots ,\;n$ и $j = 1,\; \ldots ,\;d$. Далее можно ввести несколько обобщений понятия матричного спектрального радиуса для тензоров (см., например, [46]). В этой работе мы используем следующее определение:

Доказано, что для кубических неотрицательных тензоров ${\mathbf{T}}$ таких, что $T(1)$ поэлементно положителен, где $1$ – вектор из единиц, при $p > d$ выполнено ${{r}_{p}}({\mathbf{T}}) > 0$, и существует единственный поэлементно положительный ${\mathbf{u}} \in {{\mathbb{R}}^{n}}$ такой, что ${{\left\| {\mathbf{u}} \right\|}_{p}} = 1$ и $T({\mathbf{u}}) = {{r}_{p}}({\mathbf{T}}){{\Phi }_{p}}({\mathbf{u}})$.

Определим дополнительно:

Для вычисления максимальной ${{\ell }^{p}}$-собственной пары тензора ${\mathbf{T}}$ степенной метод, возможно, лучший из известных подходов. Начнем с ${{{\mathbf{x}}}_{0}} \in {{C}_{ + }}({\mathbf{T}})$, где ${{C}_{ + }}({\mathbf{T}})$ – конус неотрицательных векторов, имеющих такой же шаблон нулей, как $T(1)$. Выберем $p > d$, точность $\varepsilon > 0$ и $q = p{\text{/}}(p - 1)$ (сопряженный показатель). Степенной метод, записанный ниже, позволяет получить последовательность $({{\lambda }_{k}})$, сходящуюся к требуемому собственному значению, а последовательность $({{{\mathbf{x}}}_{k}})$ сходится к соответствующему собственному вектору при любом выборе начального условия ${{{\mathbf{x}}}_{0}} \in {{C}_{ + }}({\mathbf{T}})$:

Данный алгоритм может работать достаточно медленно, что делает его неприменимым для решения реальных прикладных задач, например, для работы с тензорами, соответствующими сетевым данным. Таким образом, для ускорения сходимости последовательности $({{{\mathbf{x}}}_{k}})$ мы используем алгоритм STEA2 с техникой рестарта.

На графиках ниже (фиг. 10) мы отображаем результаты вплоть до $30$-й итерации для поэлементной невязки собственного вектора

как для исходной последовательности, так и для экстраполированной, и мы “подчеркиваем” крупными точками каждый рестарт внешнего цикла, например, невязку, сгенерированную итерациями экстраполяции. Мы фиксируем $p = d + {{10}^{{ - 5}}}$ и устанавливаем ${{{\mathbf{x}}}_{0}}$ случайным вектором. С выбранным $p$ мы можем вычислить приближение H-собственных пар (соответствующих $p = d$), так как в данном случае сходимость не гарантируется.

Линейный функционал ${\mathbf{y}}$ обновляется в конце каждого внешнего цикла как ${\mathbf{y}} = {{{\mathbf{\tilde {e}}}}_{h}}({{{\mathbf{x}}}_{0}})$ (для первого шага экстраполяции мы выбираем ${\mathbf{y}} = {{{\mathbf{x}}}_{0}}$).

В табл. 2 показаны четыре примера реальных сложных сетей, описываемых с помощью неотрицательных тензоров ${\mathbf{T}}$ третьего порядка, порожденных данными из работы [47].

Больше подробностей, информацию о других алгоритмах и примерах для данного типа задач можно найти в [48].

6.2. Вычисление полилинейного PageRank

Недавно идея PageRank была обобщена для случая многомерных марковских цепей [49]. Исходное многомерное распределение PageRank аппроксимируется так называемым полилинейным PageRank, а задача поиска истинного многомерного стационарного распределения заменяется задачей, связанной с вычислением симметричного тензора ранга один.

В рамках предложенной модели требуется решить следующую задачу о неподвижной точке:

где ${\mathbf{T}}$ – стохастический тензор и ${\mathbf{v}}$ – поэлементно положительный стохастический вектор, и где по определению из предыдущего раздела ${\mathbf{T}}{{{\mathbf{x}}}^{{d - 1}}} = T({\mathbf{x}})$.

Для решения данной задачи в работе [50] рассмотрены два метода решения задачи о неподвижной точке: метод простой итерации со сдвигами (кратко SFPM) [49] в формулировке, адаптированной для конкретного случая вычисления полилинейного PageRank, а также многомерный степенной метод со сдвигами [51] и дополнительно метод внутренних и внешних итераций из работы [49]. Скорость сходимости может быть очень низкой. Таким образом, мы применили топологический $\varepsilon $-алгоритм с техникой рестарта для последовательностей ${{{\mathbf{x}}}_{\ell }}$, порожденных этими методами. Многочисленные численные эксперименты на синтетических и реальных данных демонстрируют улучшение работы перечисленных методов из-за применения методов экстраполяции. Мы демонстрируем результаты для нескольких задач из работы [49] при $d = 3$, $n = 4$ или $n = 6$, и где $y \in {{\mathbb{R}}^{n}}$ в STEA2 – последний экстраполированный член с предыдущего цикла. На фиг. 11 демонстрируются графики изменения относительной первой нормы невязки.