Доклады Российской академии наук. Математика, информатика, процессы управления, 2022, T. 503, № 1, стр. 91-94

1. ВВЕДЕНИЕ

Известно, что многие вычислительные задачи принадлежат к классу некорректных или плохо обусловленных задач. Для решения этих задач широко используются многочисленные методы регуляризации. Одним из наиболее эффективных методов регуляризации продолжает оставаться метод регуляризации А.Н. Тихонова [1].

При решении многих прикладных задач, особенно разреженных большой размерности, например, таких, как обработка сигналов и изображений, анализ больших данных, важную роль играют итерационные методы регуляризации [2, 3]. Достоинством итерационных методов регуляризации является тот факт, что для них роль параметра регуляризации выполняет момент останова итераций, который легко реализуется на основе принципа невязки [2]. Итерационный принцип невязки существенно проще, с вычислительной точки зрения, принципа невязки, используемого для прямых (неитерационных) методов регуляризации.

Одним из важнейших алгоритмов построения итерационных схем регуляризации является неявный метод простых итераций на основе нормальных уравнений [3]. Это связано с тем фактом, что он обладает абсолютной сходимостью, в отличие от явных, и достаточно высокой помехоустойчивостью. В [4] был предложен специальный вариант неявной итерационной схемы на основе сингулярного разложения. Однако высокая вычислительная сложность сингулярного разложения не позволяет достаточно эффективно решать большие, в том числе разреженные, задачи.

В настоящей работе предлагается новый вариант неявной итерационной схемы, основанный на последовательном решении регуляризованных расширенных систем линейных алгебраических уравнений (СЛАУ). Этот подход позволяет существенно повысить скорость сходимости алгоритма без потери вычислительной устойчивости (за счет сохранения неизменным числа обусловленности СЛАУ на каждой итерации) неявной итерационной схемы. Данный результат обеспечивается тем фактом, что число обусловленности расширенной СЛАУ значительно меньше числа обусловленности СЛАУ неявной итерационной схемы на основе нормальных уравнений. Важнейшим достоинством предлагаемого алгоритма также является возможность решать большие разреженные задачи.

2. НОВАЯ НЕЯВНАЯ ИТЕРАЦИОННАЯ СХЕМА

Рассматривается произвольная СЛАУ

где $A \in {{\mathbb{C}}^{{m \times n}}}$, $f \in {{\mathbb{C}}^{m}}$, $u \in {{\mathbb{C}}^{n}}$, $m \geqslant n$ и ${\text{rank}}(A) = n$.

Для построения классического варианта неявной итерационной схемы [2] используется симметризация Гаусса произвольной СЛАУ (1), приводящая к нормальной системе уравнений

где A* – матрица, сопряженная к матрице A.

На основе нормальной системы (2) непосредственно получается классическая форма неявной итерационной схемы (неявный метод простых итераций) для решения произвольных СЛАУ (1):

где $\omega > 0$, ${{E}_{n}}$ – единичная матрица порядка n и $k = 0,\;1,\; \ldots $

Неявная итерационная схема (3) широко используется для получения итерационных алгоритмов регуляризации [2, 3, 5], где номер остановки по k выполняет роль параметра регуляризации.

Ниже приводится новый вариант неявной итерационной схемы.

Преобразуем (3) к виду

или

Пусть

тогда (4) можно записать в виде а (5) в виде

Введем новую переменную ${{y}_{k}} = {{\omega }^{{ - 1}}}{{r}_{k}}$. Объединяя уравнения (6) и (7), получаем расширенную СЛАУ

Из (8) непосредственно следует

Утверждение. Неявный метод простых итераций (3) эквивалентен решению последовательности регуляризованных расширенных систем

Итерационный метод (9) можно рассматривать как специальную вычислительную схему метода (3). Однако при малых α вычисления по итерационной схеме (4) могут, в виду большей численной устойчивости, оказаться намного более целесообразными, чем вычисления по итерационной схеме (3). Кроме того, формирование матрицы A*A может приводить к ее полному заполнению, даже если исходная матрица A является сильно разреженной, что в конечном итоге не позволяет использовать алгоритм (3) для решения больших разреженных задач.

Спектральное число обусловленности СЛАУ на каждой итерации в алгоритме (3) равно

где ${{\sigma }_{1}}$ и ${{\sigma }_{n}}$ – максимальное и минимальное сингулярные числа матрицы $A$ соответственно.

Как показано в [6], собственными значениями матрицы ${{B}_{\omega }}$, ${{\lambda }_{j}}({{B}_{\omega }})$, $j = 1, \ldots ,m + n$, являются числа $ \pm \sqrt {\sigma _{i}^{2} + {{\omega }^{2}}} $, $i = 1,\;2,\; \ldots ,\;\tau $, а также ω и –ω кратностей $m - \tau $ и $n - \tau $ соответственно, где $\tau = {\text{rank}}(A)$ и ${{\sigma }_{1}} \geqslant {{\sigma }_{2}} \geqslant \; \ldots \; \geqslant {{\sigma }_{\tau }} > 0$ – сингулярные числа матрицы A. Для частного случая, когда $A \in {{\mathbb{R}}^{{m \times n}}}$, аналогичный результат приведен (без доказательства) в работе М.А. Саундерса [7].

Так как ${{B}_{\omega }} = B_{\omega }^{*}$ и $\tau = n$, то спектральное число обусловленности матрицы ${{B}_{\omega }}$

где ${{\left\| {\, \cdot \,} \right\|}_{2}}$ – спектральная матричная норма.

Из (10) следует, что если ${{\sigma }_{n}} < \omega < {{\sigma }_{1}}$, то

несмотря на число обусловленности матрицы A.

Из (11) получаем важный практический результат: если ${{\sigma }_{n}} < \omega < {{\sigma }_{1}}$, то число обусловленности матрицы ${{B}_{\omega }}$ не зависит от обусловленности матрицы A и, следовательно, для плохо обусловленных матриц A число обусловленности ${{\mu }_{2}}({{B}_{\omega }}) \ll {{\mu }_{2}}({{G}_{\omega }})$, так как ${{\mu }_{2}}({{B}_{\omega }}) = \sqrt {{{\mu }_{2}}({{G}_{\omega }})} $. Этот факт гарантирует, что неявная итерационная схема (9), на основе расширенных систем, будет обладать более высокой вычислительной устойчивостью по сравнению с итерационной схемой (3), на основе нормальных систем, при одном и том же значении параметра $\omega $.

3. ИТЕРАЦИОННАЯ РЕГУЛЯРИЗИРУЮЩАЯ СХЕМА

Пусть вместо точных исходных данных $(A,f)$ заданы возмущенные исходные данные $(A,\tilde {f})$, где $\left\| {\tilde {f} - f} \right\| \leqslant \delta $. Здесь $\left\| f \right\|$ – евклидова векторная норма.

Для построения итерационной регуляризирующей схемы воспользуемся итерационным алгоритмом (9):

В качестве критерия остановки будет использовано минимальное число k, такое, что

где константа $c > 1$.

Параметр ${{k}_{\delta }}$, момент останова для итерационной схемы (12), играет роль параметра регуляризации на основе принципа невязки [2], а ${{u}_{{{{k}_{\delta }}}}}$ является регуляризованным решением.

Так как $\tilde {f} - A{{u}_{k}} = \omega {{y}_{k}}$, то критерий остановки (13) можно представить в виде

Использование (14) позволяет снизить вычислительные затраты по сравнению с (13).

4. АНАЛИЗ СХОДИМОСТИ И ЧИСЛЕННЫЙ ПРИМЕР

Итерационную схему (3) можно записать в виде

где $\Phi (\omega ) = {{\omega }^{2}}{{({{\omega }^{2}}{{E}_{n}} + A{\text{*}}A)}^{{ - 1}}}$ – итерационная матрица, а ${{g}^{\omega }} = {{({{\omega }^{2}}{{E}_{n}} + A{\text{*}}A)}^{{ - 1}}}b$.

Скорость сходимости итерационного процесса (15) определяется множителем сходимости $\rho (\Phi (\omega ))$, где $\rho ( \cdot )$ – спектральный радиус матрицы. В [4] показано, что $\rho (\Phi (\omega )) = \frac{1}{{1 + {{\lambda }^{2}}}}$, где $\lambda = \frac{{{{\sigma }_{n}}}}{\omega }$.

Очевидно, что ${{\mu }_{2}}(G(\omega )) \to 2$ и ${{\mu }_{2}}(B(\omega )) \to \sqrt 2 $ при $\omega \to \infty $, где $G(\omega ) = {{G}_{\omega }}$ и $B(\omega ) = {{B}_{\omega }}$. Однако при $\omega \to \infty $ спектральный радиус $\rho (\Phi (\omega )) \to 1$, что приводит к сильному замедлению скорости сходимости.

Видно, что параметр ω, с одной стороны, выполняет роль предобусловливателя, но, с другой стороны, требуется компромис при выборе параметра $\omega $ для обеспечения достаточной скорости сходимости алгоритма. Относительно стандартного алгоритма типа (3) замечание о скорости сходимости было сделано в недавней работе [8].

Исследуем сравнение скорости сходимости и чисел обусловленности для алгоритмов (3) и (12).

Пусть параметры ${{\omega }_{1}}$ и ${{\omega }_{2}}$ выбираются из условия ${{\mu }_{2}}(G({{\omega }_{1}})) = {{\mu }_{2}}(B({{\omega }_{2}}))$. Так как μ₂(B_ω) = = $\sqrt {{{\mu }_{2}}({{G}_{\omega }})} $, то очевидно, что ${{\omega }_{2}} < {{\omega }_{1}}$ и, следовательно, $\rho (\Phi ({{\omega }_{2}})) < \rho (\Phi ({{\omega }_{1}}))$. Это означает, что неявная итерационная схема (12), на основе расширенных систем, имеет более высокую скорость сходимости по сравнению с (3), на основе нормальных уравнений, при одном и том же числе обусловленности СЛАУ в обоих итерационных схемах.

Рассмотрим случай, когда $\mu _{2}^{2}(A)\, = \,{{\mu }_{2}}(G(0))\, > \,{{\epsilon }^{{ - 1/2}}}$, где $\epsilon $ – машинный эпсилон. В этом случае для обеспечения гарантированной точности вычислений на каждой итерации неявных итерационных схем (3) и (12) параметры ${{\omega }_{1}}$ и ${{\omega }_{2}}$ должны удовлетворять условиям ${{\mu }_{2}}(G({{\omega }_{1}}))\, \leqslant \,{{\epsilon }^{{ - 1/2}}}$ и ${{\mu }_{2}}(B({{\omega }_{2}}))\, \leqslant \,{{\epsilon }^{{ - 1/2}}}$ соответственно.

Если параметры ${{\omega }_{1}}$ и ${{\omega }_{2}}$ удовлетворяют условию

то в этом случае ${{\omega }_{2}} \ll {{\omega }_{1}}$. Отсюда непосредственно имеем, что $\rho (\Phi ({{\omega }_{2}})) \ll \rho (\Phi ({{\omega }_{1}}))$. Следовательно, при одинаковых числах обусловленности СЛАУ в алгоритмах (3) и (12), алгоритм (12), на основе расширенных линейных систем, будет иметь существенно более высокую скорость сходимости. Этот факт иллюстрируется следующим простым примером.

Численный пример. Рассмотрим СЛАУ Au = f, где

СЛАУ (17) совместна, ее точное решение ${{u}_{*}}\, = \,\left( \begin{gathered} 5 \hfill \\ 2 \hfill \\ \end{gathered} \right)$, а ${{\sigma }_{1}} = 13.19$, ${{\sigma }_{2}} = 3.21 \times {{10}^{{ - 6}}}$ и ${{\mu }_{2}}(A) \approx 4.10 \times {{10}^{6}}$. Для решения СЛАУ (17) использовались алгоритмы (3) и (12). Все вычисления выполнялись в Matlab, где $\epsilon \approx 2.22 \times {{10}^{{ - 16}}}$. В данном примере $\tilde {f} = f$. В алгоритмах (3) и (12) параметры выбирались в соответствии с (16):

Из (18) получаем, что для алгоритма (3) ω₁ = = $0.007745$, а для алгоритма (12) ${{\omega }_{2}} = 3.21 \times {{10}^{{ - 6}}}$. Для решения СЛАУ в алгоритме (3) применялся метод Холецкого, а в (12) – метод LU-разложения. В качестве критерия остановки использовался критерий:

Для алгоритма (12) этот критерий, в соответствии с (14), имеет вид $\omega \left\| {{{y}_{k}}} \right\| \leqslant 1.2 \times \epsilon $.

Результаты вычислений. Алгоритм (3), на основе нормальных уравнений, не позволил вообще выполнить критерий (19). Через $k = 1000$, невязка $\left\| {A{{u}_{{1000}}} - f} \right\| = 1.73$ × 10^–5, а относительная погрешность решения $\frac{{\left\| {{{u}_{{1000}}} - {{u}_{*}}} \right\|}}{{\left\| {{{u}_{*}}} \right\|}} \approx 0.99953$ соответственно.

Для алгоритма (12), на основе расширенных линейных систем, критерий остановки (19) был выполнен при $k = 37$, а относительная погрешность, соответственно, равна

5. ЗАКЛЮЧЕНИЕ

Новый вариант неявной итерационной схемы (9) обладает следующими наиболее важными преимуществами перед классической неявной итерационной схемой (3), основанной на нормальной системе уравнений:

1. Неявная итерационная схема (9), в отличие от неявной схемы (3), позволяет достаточно эффективно решать большие разреженные задачи.

2. При одном и том же числе обусловленности матриц ${{G}_{\omega }}$ и ${{B}_{\omega }}$ неявная итерационная схема (9) позволяет получить существенно более высокую скорость сходимости при решении плохо обусловленных задач. Этот факт особенно важен при решении задач большой размерноcти.