Неявная итерационная схема на основе алгоритма псевдообращения и её применения
- Авторы: Жданов А.И.1,2, Сидоров Ю.В.1
-
Учреждения:
- Самарский государственный технический университет
- Филиал ФГБОУ ВО «СамГТУ» в г. Новокуйбышевске
- Выпуск: Том 28, № 1 (2024)
- Страницы: 117-129
- Раздел: Математическое моделирование, численные методы и комплексы программ
- URL: https://journals.eco-vector.com/1991-8615/article/view/462798
- DOI: https://doi.org/10.14498/vsgtu2026
- EDN: https://elibrary.ru/HIGWRZ
- ID: 462798
Цитировать
Полный текст
Аннотация
Предложена новая версия неявной итерационной схемы, для реализации которой требуются лишь матрично-векторные вычислительные процедуры. Это делает предлагаемую вычислительную схему потенциально высокоэффективной для решения широкого класса задач большой размерности на современных высокопроизводительных вычислительных платформах, например Nvidia Cuda. Показано, что предлагаемые алгоритмы могут быть использованы для решения плохо обусловленных линейных систем и задач наименьших квадратов, а также для построения итерационных алгоритмов регуляризации. Приводятся результаты вычислительных экспериментов, подтверждающие эффективность предлагаемых вычислительных алгоритмов.
Полный текст
Введение
Неявные итерационные схемы широко применяются для решения различных задач вычислительной математики, таких как плохо обусловленные и сингулярные системы линейных алгебраических уравнений (СЛАУ) большой размерности [1], задач наименьших квадратов, а также в итерационных методах регуляризации решений приближенных СЛАУ [2]. Известно, что одним из существенных преимуществ неявных итерационных схем перед явными является их абсолютная сходимость и достаточно высокая вычислительная устойчивость, но эти схемы не лишены недостатков.
В [2] предложена вычислительная схема неявного итерационного метода на базе сингулярного разложения (SVD-разложения). Однако применение этой вычислительной схемы приводит к существенному увеличению вычислительной работы и усложнению алгоритмической реализации итерационного алгоритма. Кроме того, использование в алгоритме SVD-разложения приводит к тому, что этот алгоритм не может быть адаптирован для эффективной обработки разреженных матриц. В конечном итоге все это серьезно осложняет практическое применение данного алгоритма для решения задач большой размерности на современных высокопроизводительных вычислительных платформах.
В [3] предложена еще одна неявная итерационная схема на основе расширенных линейных систем. Эта схема представляет большой интерес для решения плохо обусловленных больших разреженных линейных систем и задач наименьших квадратов, а также итерационных алгоритмов регуляризации. Однако для ее применения требуется использовать известные прямые методы решения («решатели») СЛАУ.
В настоящей работе предлагается новый вариант неявного метода простых итераций с использованием итерационного алгоритма псевдообращения Бен– Израэля [4]. Этот вариант алгоритма имеет полностью итерационный характер и обладает алгоритмическая простотой. Для реализации этого метода требуются лишь матрично-векторные вычислительные процедуры, при его применении нет необходимости использования дополнительных стандартных «решателей» СЛАУ. Этот факт делает его потенциально более эффективным, чем указанные методы [2, 3], для решения задач большой размерности на современных высокопроизводительных вычислительных платформах, например Nvidia Cuda.
С вычислительной точки зрения наиболее затратной составляющей предлагаемого метода является итерационный метод Бен– Израэля. Однако учитывая специальную структуру матрицы с использованием «регуляризирующего» параметра, который одновременно выполняет роль предобуславливателя в задаче, данный метод, как показывают многочисленные вычислительные эксперименты, обладает высокой скоростью сходимости и вычислительной устойчивостью.
1. Постановка задачи
Рассмотрим задачу решения СЛАУ вида
\[ \begin{equation}
Au=f,
\end{equation} \tag{1} \]
где $A\in\mathbb{R}^{m \times n}$, $m \geq n$, $\operatorname{rank} A=n$.
Псевдорешение исходной СЛАУ (1) определяется как $u_*=A^{+}f$, где $A^+$ — псевдообратная матрица Мура– Пенроуза [4].
Отметим, что $ A^{+}=(A^\top A)^{-1}A^\top$, поскольку предполагается, что матрица $A$ имеет полный столбцовый ранг, т.е. $\operatorname{rank} A=n$, и, следовательно, является левой обратной матрицей.
В [2] показано, что в общем виде неявный метод простых итераций можно представить выражением
\[ \begin{equation}
u_{k+1}=\mathop{\operatorname{argmin}}_{u\in\mathbb{R}^{n}}
\left\|\left(
\begin{array}{c}
A \\
\sqrt{\alpha}E_n \\
\end{array}
\right)u-
\left(
\begin{array}{c}
f \\
\sqrt{\alpha}u_k \\
\end{array}
\right)
\right\|_2^2,
\end{equation} \tag{2} \]
где
\[ \begin{equation*}
\left(\begin{array}{c}
A \\
\sqrt{\alpha}E_n \\
\end{array}\right)\in\mathbb{R}^{(m+n)\times n} , \qquad
\left(\begin{array}{c}
f \\
\sqrt{\alpha}u_k \\
\end{array}\right)\in\mathbb{R}^{m+n},
\end{equation*} \]
$E_n$ — единичная матрица порядка $n$, $\|\cdot\|_2$ — евклидова векторная норма, ${\alpha>0}$, $k=0, 1, \dots$.
Из (2) видно, что различные варианты реализации неявного метода простых итераций отличаются лишь методом решения на каждой итерации линейной задачи наименьших квадратов:
\[ \begin{equation}
\min_{u\in\mathbb{R}^{n}}\left\|\left(
\begin{array}{c}
A \\
\sqrt{\alpha}E_n \\
\end{array}
\right)u-
\left(
\begin{array}{c}
f \\
\sqrt{\alpha}u_k \\
\end{array}
\right)
\right\|_2^2.
\end{equation} \tag{3} \]
В зависимости от методов решения задачи (3) получаются различные варианты неявных итерационных схем.
Первый подход основан на применении метода нормальных уравнений. В этом случае получаем классическую форму неявного метода простых итераций, которая была предложена J.D. Riley в [5]:
\[ \begin{equation}
(A^\top A+\alpha E_n)u_{k+1}=A^\top f+\alpha u_k.
\end{equation} \tag{4} \]
Второй подход для решения задачи наименьших квадратов (3) предполагает использование метода, основанного на QR-разложении.
Третий подход основан на SVD-разложении. Этот подход рассмотрен в работе [2].
В [3] показано, что неявная итерационная схема (4) решения задачи (3) эквивалентна решению последовательности регуляризованных расширенных систем
\[ \begin{equation}
\left(
\begin{array}{cc}
\omega E_m & A \\
A^{*} & -\omega E_n \\
\end{array}
\right)
\left(
\begin{array}{c}
y_{k+1} \\
u_{k+1} \\
\end{array}
\right)=
\left(
\begin{array}{c}
f \\
-\omega u_k \\
\end{array}
\right),
\end{equation} \tag{5} \]
где $y_k=\omega^{-1}(f-Au_k)$, $\omega=\sqrt{\alpha}$.
Итерационная схема (5) имеет существенно более высокую скорость сходимости по сравнению с неявной итерационной схемой (4) и значительно меньшее число обусловленности вычислительной задачи. Важнейшим достоинством итерационной схемы (5) является возможность решать большие разреженные СЛАУ. Это обусловлено тем, что для решения расширенных систем (5) можно использовать хорошо известные «решатели», адаптированные для решения разреженных СЛАУ.
2. Новый вариант неявного итерационного метода
В настоящей работе предлагается новый вариант неявного метода простых итераций с использованием итерационного алгоритма псевдообращения Бен–Израэля.
Пусть
\[ \begin{equation*}
\left\|\left(\begin{array}{c}
A \\
\sqrt{\alpha}E_n \\
\end{array}
\right)u-
\left(
\begin{array}{c}
f \\
\sqrt{\alpha}u_k \\
\end{array}
\right)
\right\|_2^2=
\|A_{\alpha}u-f^{\alpha}_{k}\|_2^2,
\end{equation*} \]
где
\[ \begin{equation*}
A_\alpha=\left(\begin{array}{c}
A \\
\sqrt{\alpha}E_n \\
\end{array}\right)\in\mathbb{R}^{(m+n)\times n} , \qquad
f^{\alpha}_{k}=\left(\begin{array}{c}
f \\
\sqrt{\alpha}u_k \\
\end{array}\right)
\in\mathbb{R}^{m+n}.
\end{equation*} \]
Как показано в [2], итерационную схему (2) можно представить в виде
\[ \begin{equation}
u_{k+1}=A_{\alpha}^{+}f^{\alpha}_{k}.
\end{equation} \tag{6} \]
Пусть $A^+_\alpha=(U_\alpha\vdots V_\alpha)\in\mathbb{R}^{n \times (m+n)}$, где $U_\alpha\in\mathbb{R}^{n \times m}$, $V_\alpha\in\mathbb{R}^{n \times n}$. Тогда неявный метод простых итераций (6) можно записать в виде
\[ \begin{equation}
u_{k+1}=U_\alpha f+\sqrt{\alpha}V_\alpha u_k.
\end{equation} \tag{7} \]
Для последовательности $(u_k)_{k=0}^\infty$, определяемой итерационной формулой (7) при произвольном начальном $u_0$, справедливо $\lim_{k \to \infty}\|u_k-u_*\|_2=0$, где $u_*$ — псевдорешение СЛАУ (1).
Так как матрица $A_\alpha$ всегда имеет полный столбцовый ранг, т.е. ${\operatorname{rank} A_\alpha=n}$, для вычисления $A_\alpha^+$ можно воспользоваться итерационным методом Бен–Израэля [4]:
\[ \begin{equation}
X_{i+1}=(2E_n-X_iA_\alpha)X_i,
\end{equation} \tag{8} \]
где $X_i\in\mathbb{R}^{n \times (m+n)}$, $i=0, 1, 2, \dots$.
Итерационный метод Бен–Израэля (8) имеет квадратичную скорость сходимости [4] и, как отмечено в [6, 7], высокую численную устойчивость.
Пусть $X_i=(Y_i \ \vdots \ Z_i )$, где $Y_i\in\mathbb{R}^{n \times m}$, $Z_i\in\mathbb{R}^{n \times n}$. Тогда (8) можно записать в виде
\[ \begin{equation}
X_{i+1}=\bigl[2E_n-(Y_iA+\sqrt{\alpha}Z_i)\bigr]X_i
\end{equation} \tag{9} \]
или
\[ \begin{equation*}
(Y_{i+1} \ \vdots \ Z_{i+1})=(2E_n-Y_iA-\sqrt{\alpha}Z_i)(Y_i \ \vdots \ Z_i).
\end{equation*} \]
Если $X_0=\beta A_\alpha^\top=\beta(A^\top \ \vdots \ \sqrt{\alpha}E_n)$ или $Y_0=\beta A^\top$, $Z_0=\beta\sqrt{\alpha}E_n$ и
\[ \begin{equation*}
0<\beta<\frac{2}{\sigma_1^2(A_\alpha)}=\frac{2}{\sigma_1^2(A)+\alpha},
\end{equation*} \]
где $\sigma_1(A_\alpha)$ и $\sigma_1(A)$ — максимальные сингулярные числа матриц $A_\alpha$ и $A$ соответственно, то последовательность (8) при $i\to\infty$ сходится к $A_\alpha^{+}$ [4].
Так как $\sigma_1^2=\|A\|_2^2\leqslant\|A\|_F^2$, можно принять, что
\[ \begin{equation*}
\beta=\frac{1.8}{\|A\|_F^2+\alpha},
\end{equation*} \]
где $\|A\|_F$ — сферическая матричная норма (норма Фробениуса), $\|A\|_2$ — спектральная матричная норма [8, § 2.3].
Критерий остановки итерационного процесса (9): $i=0, 1, \dots,i(\delta)$, где
\[ \begin{equation*}
i(\delta)=\min \Bigl\{i\in\mathbb{N}: \frac{\|X_{k+1}-X_k\|_\infty}{1+\|X_k\|_\infty}\leqslant\delta \Bigr\},
\end{equation*} \]
$\delta$ — априори заданное число для остановки итерационного процесса (8) по относительной скорости сходимости, $\|X\|_\infty$ — матричная норма, подчиненная кубической векторной норме ($\infty$-норме) [8, § 2.3].
Тогда (7) принимает вид
\[ \begin{equation}
\tilde{u}_{k+1}=\sqrt{\alpha}V_\delta\tilde{u}_{k}+g_\delta,
\end{equation} \tag{10} \]
где $X_\delta=X_{i_{(\delta)}}$, а $U_\delta=Y_{i_{(\delta)}}$, $V_\delta=Z_{i_{(\delta)}}$ и $g_\delta=U_\delta f$.
3. Вычислительные эксперименты
Проведем сравнение предлагаемого в статье итерационного алгоритма (10) для решения задач вида (1) с методами, основанными на вычислении псевдообратной матрицы $A^+$, т.е. когда решения находятся по формуле $u_*=A^{+}f$.
В первом случае псевдообратная матрица $A^+$ вычисляется на основе хорошо известного алгоритма сингулярного разложения (SVD-разложения).
Во втором случае псевдообратная матрица $A^+$ вычисляется известным итерационным методом Бэн– Израэля. В случае, когда матрица $A$ задачи (1) — невырожденная квадратная матрица, для вычисления обратной матрицы применяется итерационный метод Шульца [9], который является частным случаем итерационного алгоритма псевдообращения Бен– Израэля.
Все дальнейшие вычисления выполняются на Python 3.7.4.
3.1. Решение плохо обусловленных совместных СЛАУ
Рассмотрим решение плохо обусловленных совместных СЛАУ вида (1) с квадратными матрицами $A\in\mathbb{R}^{n \times n}$.
Для первого вычислительного эксперимента матрица $A$ системы (1) генерировалась с помощью функции deriv2 из пакета RTools [10], вектор точного решения брался в виде $u_{\mathrm{ext}}=\left(1,\; 2\;, \ldots,\; n \right)^\top$, а вектор правой части определялся по формуле $f=Au_{\mathrm{ext}}$. Спектральное число обусловленности матрицы $A$ определяется по формуле $\kappa_2 (A ) = {\sigma_1}/{\sigma_n}$, где $\sigma_1$ и $\sigma_n$ — максимальное и минимальное сингулярные числа матрицы $A$ соответственно.
Решение задачи (1) находилось с помощью следующих подходов:
- решение вычислялось как $u_{\mathrm{svd}}=A^{+}f$, где матрица $A^{+}$ вычислялась на основе SVD-алгорима;
- нахождение обратной матрицы $A^{-1}$ с помощью итерационного алгоритма Шульца и вычисление решения задачи (1);
- вычисление решения задачи (1) с помощью предлагаемого в статье алгоритма (10).
Критерий остановки итерационных алгоритмов Шульца и Бен– Израэля:
\[ \begin{equation*}
\frac{ \Vert X_{k+1}-X_k \Vert_{\infty}}{1+ \Vert X_k \Vert_{\infty}}\leqslant 10^{-7}.
\end{equation*} \]
Критерий остановки итерационного алгоритма (10):
\[ \begin{equation*}
\mathrm{Err} =
\frac{\Vert u_{k+1}-u_k \Vert_{\infty}}{1+ \Vert u_k \Vert_{\infty}}\leqslant 10^{-16}.
\end{equation*} \]
При использовании итерационного алгоритма (10) для рассматриваемой задачи параметр $\omega=\sqrt{\alpha}$ выбирался равным $\sigma_n/2$, $\sigma_n$, $2 \sigma_n$, $3\sigma_n$.
В табл. 1 и на рис. 1 представлены результаты вычислений для задачи (1), когда $n=512$, $\kappa_2 (A ) = 3.19\cdot 10^5$, $\sigma_n= 3.17\cdot 10^{-7}$, полученные при использовании итерационного алгоритма (10). Здесь и далее $\kappa_2(A_{\alpha} )$ — спектральное число обусловленности матрицы $A_{\alpha}$, $\mathrm{it}_1$ — число итераций алгоритма (10), $\mathrm{it}_2$ — число итераций алгоритма Бен– Израэля, $\mathrm{RErr} = \frac{\Vert u_{\mathrm{ext}}-u_k\Vert_{2}}{ \Vert u_{\mathrm{ext}} \Vert_{2}}$ — относительная ошибка вычислений.
$\omega$ | $\kappa_2 (A_{\alpha} )$ | $\mathrm{it}_1$ | $\mathrm{it}_2$ | $\mathrm{RErr}$ |
$\sigma_n/2$ | $2.85\cdot 10^5$ | 23 | 41 | $1.90 \cdot 10^{-11}$ |
$\sigma_n$ | $2.25\cdot 10^5$ | 53 | 40 | $1.88 \cdot 10^{-11}$ |
$2\sigma_n$ | $1.43\cdot 10^5$ | 151 | 39 | $1.52 \cdot 10^{-11}$ |
$3\sigma_n$ | $1.01\cdot 10^5$ | 309 | 38 | $2.16\cdot 10^{-11}$ |
Рис. 1. Скорость сходимости итерационного алгоритма (10) при решении плохо обусловленных совместных СЛАУ (1)
[Figure 1. Convergence rate of the iterative algorithm (10) when solving ill-conditioned systems of linear equations (1)]
Отметим, что для этой же задачи для SVD-алгоритма получена относительная ошибка вычислений, равная $1.62 \cdot 10^{-10}$, а для алгоритма Шульца при числе итераций, равном 41, получена относительная ошибка вычислений, равная $5.75 \cdot 10^{-8}$.
Данный вычислительный эксперимент демонстрирует, что предлагаемый в работе итерационный алгоритм (10) по точности превосходит метод на основе итерационного алгоритма Шульца и не уступает методам, основанным на сингулярном разложении, применяемым для решения подобных плохо обусловленных совместных СЛАУ. Итерационный метод Бен–Израэля вычисления псевдообратной матрицы для предложенного способа выбора параметра $\omega$, имеет достаточно высокую скорость сходимости и не сильно зависит от параметра $\omega$.
Рассматриваемый в статье итерационный алгоритм (10) при выборе параметра $\omega$, близким к минимальному сингулярному числу $\sigma_n$ матрицы $A$, позволяет вычислять решения плохо обусловленных задач вида (1) с достаточно высокой скоростью и точностью. Увеличение параметра $\omega$ для задач данного класса приводит к значительному замедлению итерационного процесса, что, в свою очередь, может вызвать потерю точности получаемого решения.
Как отмечалось ранее, одними из преимуществ предлагаемого алгоритма (10) являются алгоритмическая простота и возможность реализации итерационного процесса только с помощью матрично-векторных операций.
3.2. Решение плохо обусловленных линейных задач наименьших квадратов
Рассмотрим плохо обусловленную линейную задачу наименьших квадратов
\[ \begin{equation}
A u =f,
\end{equation} \tag{11} \]
где
\[ \begin{equation*}
A = \left(\begin{array}{ccccc}
1 & 1 & 1 & 1 & 1 \\
10^{-8} & 0 & 0 & 0 & 0 \\
0 & 10^{-8} & 0 & 0 & 0 \\
0 & 0 & 10^{-8} & 0 & 0 \\
0 & 0 & 0 & 10^{-8} & 0 \\
0 & 0 & 0 & 0 & 10^{-8} \\
\end{array}\right),
\end{equation*} \]
$u = (1,\; 1,\; 1,\; 1,\; 1 )^\top$, $f = A u + r_0$, $r_0\in\operatorname{null} (A^\top )$.
Учитывая, что вектор невязки $r_0\in\operatorname{null} (A^\top )$, псевдорешение $u_*=A^{+}f$ системы (11) равно вектору $u$. Спектральное число обусловленности матрицы $A$ задачи (11) равно $\kappa_2(A)=2.24\cdot 10^8$.
Решение задачи (11) с помощью псевдообратной матрицы, вычисленной с применением SVD-алгоритма, получено с относительной ошибкой вычислений, равной $9.42 \cdot 10^{-10}$. Решение этой задачи итерационным алгоритмом Бен–Израэля при числе итераций, равном 60, получено с относительной ошибкой вычислений, равной $1.46 \cdot 10^{-9}$. При этом в качестве критерия остановки итерационного алгоритма Бен–Израэля использовалось условие
\[ \begin{equation*}
\frac{ \Vert X_{k+1}-X_k \Vert_{\infty}}{1+ \Vert X_k \Vert_{\infty}}\leqslant 10^{-7}.
\end{equation*} \]
В табл. 2 представлены результаты вычислений для задачи (11), полученные с помощью алгоритма (10). Здесь $\sigma_5$ — минимальное сингулярное значение матрицы $A$ для задачи (11). При этом в качестве критерия остановки итерационного процесса использовалось условие
\[ \begin{equation*}
\frac{\Vert u_{k+1}-u_k\Vert_{\infty}}{1+\Vert u_k \Vert_{\infty}}\leqslant 10^{-16} .
\end{equation*} \]
$\omega$ | $\kappa_2 (A_{\omega} )$ | $\mathrm{it}_1$ | $\mathrm{it}_2$ | $\mathrm{RErr}$ |
$\sigma_1$ | $\approx 1$ | 64 | 7 | $5.98\cdot 10^{-15}$ |
$\sigma_1 /100$ | $\approx 1\cdot 10^2$ | 7 | 18 | $ 2.67 \cdot 10^{ -16}$ |
$\sigma_5$ | $\approx 1\cdot 10^8$ | 30 | 59 | $3.67 \cdot 10^{-8}$ |
Анализ результатов вычислительного эксперимента демонстрирует, что при решении плохо обусловленных линейных задач наименьших квадратов предлагаемый алгоритм (10) по точности не уступает известным методам решения этой задачи, а соответствующий выбор параметра $\omega$ позволяет получать более точные решения. Алгоритм Бен–Израэля для данного вычислительного эксперимента, как и для предыдущего, демонстрирует достаточно высокую скорость сходимости и близкую к линейной зависимость от параметра $\omega$.
3.3. Решение задач c возмущенными исходными данными
В рассмотренных выше вычислительных экспериментах предложенный итерационный алгоритм (10) применялся для нахождения решения задач вида (1) с точными исходными данными $(A,\; f)$. Очень часто при решении различных прикладных задач вместо точных данных $(A,\; f)$ имеются возмущенные исходные данные $(A,\; \tilde{f})$. Наиболее широко применяемыми на практике подходами для решения задач вида (1) c возмущенными исходными данными $(A,\; \tilde{f})$ являются методы регуляризации.
Рассматриваемый итерационный метод (10) при выборе критерия остановки, в котором номер итерации играет роль параметра регуляризации (правило остановки по невязке), можно отнести к итерационным алгоритмам регуляризации.
Проиллюстрируем регуляризирующие свойства итерационного алгоритма (10) на следующей задаче решения СЛАУ:
\[ \begin{equation}
A u = f,
\end{equation} \tag{12} \]
где
\[ \begin{equation*}
A = \frac{1}{2}\left(
\begin{array}{cc}
1 & 1\\
1+10^{-8} & 1-10^{-8}
\\
\end{array}\right), \quad
f = (1, \; 1)^\top.
\end{equation*} \]
Максимальное и минимальное сингулярные числа матрицы $A$ для задачи (12) соответственно равны $\sigma_{1} = 1$ и $\sigma_{2} = 5\cdot 10^{-9}$, а спектральное число обусловленности задачи $\kappa_2(A) = 2\cdot 10^8$. Рассматриваемая задача относится к плохо обусловленным вычислительным задачам регрессионного анализа.
Вектор $u = (1,\; 1 )^\top$ является точным решением задачи (12).
Рассмотрим задачу (12) с небольшим возмущением в векторе правой части:
\[ \begin{equation}
A u = \tilde{f},
\end{equation} \tag{13} \]
где $\tilde{f} = f + (0.01 , \; 0)^\top$.
Решая систему (13) известными методами, например с помощью процедуры solve библиотеки NumPy, получим $u_* = (-10^6, \; 10^6)^\top$.
В работе [11] для решения возмущенной задачи (13) предлагалось использовать усеченное сингулярное разложение (TSVD). Для возмущенной задачи (13) в [11] параметр усечения был принят равным $\mathrm{tol}=10^{-7}$ и получено решение $u_{\mathrm{tsvd}} = (1.0050,\; 1.0049 )^\top$.
Найдем решение возмущенной задачи (13) с помощью итерационного алгоритма (10) со следующими параметрами: $\delta = 0.01$, а $\omega = \sqrt{\alpha}$ со значениями $\sigma_1$, $\sigma_1/2$, $\sigma_1/5$.
Критерий остановки итерационного процесса примем в виде
\[ \begin{equation}
k_{\delta} = \min \bigl\lbrace k\in \mathbb{N} \; : \; \Vert Au_k - \tilde{f}\Vert\leqslant \tau\delta \bigr\rbrace,
\end{equation} \tag{14} \]
где $\tau \geqslant 1$ — параметр демпфирования (в вычислительном эксперименте использовалось значение $\tau=1.01$). Выражение (14) представляет собой итерационный принцип невязки Тихонова.
Результаты вычислений представлены в табл. 3 и на рис. 2.
$\omega$ | Number of iterations, $k$ | ${\Vert u - u_k \Vert}/{ \Vert u \Vert}$ | $u_k$ |
$\sigma_1$ | 8 | $1.07\cdot 10^{-3}$ | $(1.0011,\;1.0011 )^\top$ |
$\sigma_1/2$ | 4 | $3.39\cdot 10^{-3}$ | $(1.0033,\;1.0033 )^\top$ |
$\sigma_1/5$ | 2 | $3.51\cdot 10^{-3}$ | $(1.0035,\;1.0035 )^\top$ |
Рис. 2. Скорость сходимости итерационного процесса при решении возмущенной задачи (13); горизонтальная линия соответствует значению $\tau \delta$
[Figure 2. Convergence rate of the iterative process when solving the perturbed problem (13); the horizontal line corresponds to the value $\tau \delta$]
Полученные результаты позволяют сделать вывод, что предложенный алгоритм (10) с критерием остановки по итерационному принципу невязки (14) по точности не уступают методу регуляризации, основанному на сингулярном разложении. При этом важно отметить, что для решения задачи итерационной регуляризации можно использовать оценку максимального сингулярного числа $\sigma_1$, например $\omega_1=\|A\|_F$, которая может быть легко вычислена в отличие от оценки минимального сингулярного числа $\sigma_n$.
Заключение
Одним из преимуществ предлагаемого в статье итерационного алгоритма является его простая алгоритмическая структура. Основными операциями алгоритма являются матрично-векторные произведения, для которых существуют эффективные программные реализации на современных высокопроизводительных платформах. Например, данный алгоритм может быть реализован с использованием высокопроизводительных библиотек Intel® Math Kernel Library [12] для многоядерных/многопроцессорных систем, CUDA Toolkit [13] для систем с графическими ускорителями общего назначения или адаптирован к конкретной вычислительной платформе.
Предлагаемый в статье алгоритм позволяет получать высокоэффективные алгоритмы итеративной регуляризации, основанные на итерационном принципе невязки (discripancy principle [14]).
Конкурирующие интересы. Конкурирующих интересов не имеем.
Авторская ответственность. Мы несем полную ответственность за предоставление окончательной версии рукописи в печать. Окончательная версия рукописи нами
одобрена.
Финансирование. Исследование выполнялось без финансирования.
Об авторах
Александр Иванович Жданов
Самарский государственный технический университет; Филиал ФГБОУ ВО «СамГТУ» в г. Новокуйбышевске
Автор, ответственный за переписку.
Email: zhdanovaleksan@yandex.ru
ORCID iD: 0000-0001-6082-9097
https://www.mathnet.ru/person41724
доктор физико-математических наук, профессор; профессор; каф. прикладной математики и информатики; профессор; каф. электроэнергетики, электротехники и автоматизации технологии процессов
Россия, 443100, Самара, ул. Молодогвардейская, 244; 446200, Новокуйбышевск, ул. Миронова, 5Юрий Вячеславович Сидоров
Самарский государственный технический университет
Email: linuxboy2007@gmail.com
ORCID iD: 0000-0002-8138-9200
https://www.mathnet.ru/person114787
кандидат физико-математических наук; доцент; каф. прикладной математики и информатики
Россия, 443100, Самара, ул. Молодогвардейская, 244Список литературы
- Sun L.,Wei Y., Zhou J. On an iterative method for solving the least squares problem of rankdeficient systems // Int. J. Comp. Math., 2015. vol. 92, no. 3. pp. 532–541. DOI: https://doi.org/10.1080/00207160.2014.900173.
- Zhdanov A. I. Implicit iterative schemes based on singular decomposition and regularizing algorithms // Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2018. vol. 22, no. 3. pp. 549–556. EDN: PJITAX. DOI: https://doi.org/10.14498/vsgtu1592.
- Жданов А. И. Неявная итерационная схема на основе расширенных линейных систем // Докл. РАН. Матем., информ., проц. упр., 2022. Т. 503. С. 91–94. EDN: COMJBF. DOI: https://doi.org/10.31857/S2686954322020205.
- Ben–Israel A., Greville T. N. E. Generalized Inverses: Theory and Applications. New York: Springer, 2003. xv+420 pp. DOI: https://doi.org/10.1007/b97366.
- Riley J. D. Solving systems of linear equations with a positive definite, symmetric, but possibly ill-conditioned matrix // Math. Comp., 1955. vol. 9. pp. 96–101. DOI: https://doi.org/10.1090/S0025-5718-1955-0074915-1.
- Esmaeili H., Erfanifar R., Rashidi M. An efficient Schulz-type method to compute the Moore–Penrose inverse // Int. J. Industr. Math., 2018. vol. 10, no. 2. pp. 221–228.
- Toutounian F., Soleymani F. An iterative method for computing the approximate inverse of a square matrix and the Moore–Penrose inverse of a non-square matrix // Appl. Math. Comp., 2013. vol. 224. pp. 671–680. DOI: https://doi.org/10.1016/j.amc.2013.08.086.
- Golub G. H. van Loan C. F. Matrix Computations. Baltimore: Johns Hopkins Univ., 2013. xxiv+756 pp.
- Schulz G. Iterative Berechung der reziproken Matrix // ZAMM, 1933. vol. 13, no. 1. pp. 57–59. DOI: https://doi.org/10.1002/zamm.19330130111.
- Hansen P. C. REGULARIZATION TOOLS: A Matlab package for analysis and solution of discrete ill-posed problems // Numer. Algor., 1994. vol. 6. pp. 1–53. DOI: https://doi.org/10.1007/BF02149761.
- Pillonetto G., Chen T., Chiuso A., et al. Regularized System Identification: Learning Dynamic Models from Data. Cham: Springer, 2022. xxiv+377 pp. DOI: https://doi.org/10.1007/978-3-030-95860-2.
- Wang E., Zhang Q., Shen B. et al. Intel math kernel library / High-Performance Computing on the Intel® Xeon Phi™. Cham: Springer, 2014. pp. 167–188. DOI: https://doi.org/10.1007/978-3-319-06486-4_7.
- Fatica M. CUDA toolkit and libraries / 2008 IEEE Hot Chips 20 Symposium (24–26 August 2008). Stanford, CA, USA, 2008. pp. 1–22. DOI: https://doi.org/10.1109/HOTCHIPS.2008.7476520.
- Zare H., Hajarian M. Determination of regularization parameter via solving a multiobjective optimization problem // Appl. Num. Math., 2020. vol. 156. pp. 542–554. DOI: https://doi.org/10.1016/j.apnum.2020.05.021.
Дополнительные файлы
