Алгоритм Качмажа: различия между версиями

[непроверенная версия]

[досмотренная версия]

Текущая версия на 15:43, 12 июля 2022

Алгоритм Качмажа (S. Kaczmarz)
Последовательный алгоритм
Последовательная сложность	[math]k \cdot O(n)[/math], где [math]k[/math] - количество итераций, [math]n[/math] - количество неизвестных.
Объём входных данных	[math]m \cdot n + m + n[/math], где [math]m[/math] - количество уравнений в системе.
Объём выходных данных	[math]n[/math].

Основные авторы описания: Е.С. Козлова, А.А. Иванов.

Синонимы названия алгоритма: метод алгебраической реконструкции, Algebraic Reconstruction Technique, Projection onto Convex Sets.

1 Свойства и структура алгоритмов

1.1 Общее описание алгоритма

Метод алгебраической реконструкции - полностью последовательный строковый (row-action method) метод с длинной историей, богато описанный в литературе. В классическом виде метод был впервые сформулирован и предложен польским математиком С. Качмажем (Stefan Kaczmarz) для решения систем линейных алгебраических уравнений с квадратными и невырожденными матрицами, а в последствии независимо был использован Р. Гордоном и др. (R. Gordon; R. Bender; Herman) для реконструкции изображений по их проекционным данным в задаче рентгеновской компьютерной томографии.

В более общем случае последовательность, полученная по методу итераций С. Качмажа, также сходится к решению и для переопределенных, но обязательно совместных систем линейных алгебраических уравнений. На каждой итерации алгоритма используется только одно уравнение системы линейных алгебраических уравнений, но в результате пересчитывается каждая компонента приближенного решения на данной итерации.

Далее представим наивный вывод итерационного метода следуя общим представлениям. Запишем систему линейных алгебраических уравнений с невырожденной матрицей коэффициентов [math]A \cdot u=f,\, A\in \mathbb{R}^{n\times n},\, u\in \mathbb{R}^{n},\, f\in \mathbb{R}^{n}[/math]. Проследим процедуру минимизации с использованием метода оптимального координатного спуска, описанную, например, в ^[1]. Рассмотрим функцию ошибки: [math]\Phi\left(u\right)=\Vert u_{*}-u\Vert^{2}[/math], где [math]u_{*}[/math] - единственное решение системы в смысле метода наименьших квадратов.

Будем искать решение системы на некотором линейном подпространстве [math]\Lambda\left(\varphi_{1},\ldots,\varphi_{p}\right)[/math] в виде итерационной последовательности: [math]u^{0}=\Theta,\, u^{k+1}=u^{k}\,+\, C_{k}\varphi_{j\left(k\right)}[/math], где [math]j\left(k\right)\in\left\{ 1,2,\ldots,p\right\}[/math], а [math]C_{k}[/math] - некоторый неизвестный коэффициент. Легко заметить, что [math]\Phi\left(u^{k+1}\right)=\Vert u_{*}-u^{k+1}\Vert^{2}=\Vert u_{*}-u^{k}-C_{k}\varphi_{j\left(k\right)}\Vert^{2}[/math]

Обозначим вектор [math]z^{k}=u_{*}-u^{k}[/math] - вектор ошибок приближения решения системы уравнений на шаге с номером [math]k[/math]. Предположим также, что [math]\Vert\varphi_{i}\Vert=1,\,\forall i=\overline{1,p}[/math] и [math]p=m[/math]. С учетом этого: [math]\Phi\left(u^{k+1}\right)=\Vert z^{k}-C_{k}\varphi_{j\left(k\right)}\Vert^{2}=\Vert z^{k}\Vert^{2}-2C_{k}\left(z^{k},\varphi_{j\left(k\right)}\right)+C_{k}^{2}[/math].

Продифференцируем последнее уравнение по [math]C_{k}[/math] и приравняем к нулю, в результате получим, что [math]C_{k}=\left(z^{k},\varphi_{j\left(k\right)}\right)[/math].

Возвращаясь к системе уравнений, пусть вектора в линейной оболочке [math]\Lambda[/math] выбраны как [math] \varphi_{i}=\frac{a_{i}}{\Vert a_{i}\Vert}[/math], тогда [math]\left(z^{k},a_{j\left(k\right)}\right)=f_{j\left(k\right)}-a_{j\left(k\right)}^{T}u^{k}[/math] и последовательность приобретает законченный вид. Основная формула алгоритма:

[math]u^{k+1}=u^{k}+\frac{f_{j\left(k\right)}-\left(a_{j\left(k\right)},u^{k}\right)}{\Vert a_{j\left(k\right)}\Vert^{2}}a_{j\left(k\right)}[/math].

Мы показали достаточно наивное построение алгоритма, которое не доказывает ни факт сходимости, ни тем более к решению задачи - это предмет доказательства теоремы о сходимости при определенном задании отображения [math]j\left(k\right)[/math] (при построении мы не накладывали на него никаких ограничений). Однако данное построение позволяет проследить идею о минимизации ошибки решения (а не невязки системы) и связь алгоритма Качмажа с методом координатного спуска по векторам [math]\varphi_{i}[/math].

Все дальнейшие модификации данной формулы, а, соответственно, и самого алгоритма, основаны на:

введения параметра релаксации,
экзотических способах задания отображения [math]j\left(k\right)[/math],
выборе специальных критериев останова итераций,
переходе к блочным модификациям.

1.1.1 Библиография

Стефан Качмаж - польский математик, родился 1895 г. в городе Самбор (Sambor) Галиция (Galicia), Австро-Венгрия (Austria-Hungary), ныне Украина. Стефан Качмаж был профессором математики на факультете машиностроения университета Яна Казимира (Jan Kazimierz University) в городе Львов (Lwów) с 1919 по 1939 год, там же сотрудничал с Стефаном Банахом. Предложенный им итерационный метод решения системы линейных алгебраических уравнений послужил основой для многих современных технологий в обработке изображений, в том числе в задачах компьютерной томографии ^[2].

Применение метода С. Качмажа или так называемого метода алгебраической реконструкции наиболее известно в задачах реконструктивной компьютерной томографии высокого разрешения. Для интерпретации данных, полученных в результате неразрушающего сканирования некоторого объекта, с 60-70-х годов применяются различные аппаратные и программные реализации итерационного метода, разработанного польским математиком С. Качмажем (S. Kaczmarz) в 1937 году^[3]; этот метод был использован Г. Хаусфилдом (G.N. Hounsfield) - конструктором первого томографа «ЭМИ-сканнер»^[4] (нобелевская премия по физиологии и медицине в 1979 году).

Так же данный метод можно часто встретить под следующими названиями: Algebraic Reconstruction Technique (ART)^[5]^[6]^[7], Projection onto Convex Sets (PCS)^[8], Kaczmarz Algorithm ^[9] или метод ортогональных проекций.

Конструированию различных алгоритмов на основе итераций С. Качмажа посвящено существенное количество работ: согласно исследованию польского математика и библиографа А. Сигелски (A. Cegielski), на начало 2010 года статья С. Качмажа цитируется более, чем в четырехстах реферируемых и значимых публикациях ученых из различных областей знаний. В Советской и Российской научных школах исследования итераций С. Качмажа представлены в литературе такими известными учеными как: В.Н. Ильин, Г.П. Васильченко, А.А. Светлаков и другие, например, в работах^[10]^[11]. Общие идеи о проекционных методах также представлены в классической работе Л. Г. Гурина, Б. Т. Поляка, Э. В. Райка^[12], где впервые предлагается и обосновывается идентичный алгоритм для решения систем линейных алгебраических неравенств и уравнений.

В 2008 году американскими математиками Т. Штромером (T. Strohmer) и Р. Вершининым (R. Vershynin) в ^[13] была предложена рандомизированная модификация итерационного метода, в которой последовательность приближений к решению для совместных и переопределенных систем линейных алгебраических уравнений сходится с экспоненциальной скоростью.

Множество модификаций идеи метода ортогональных проекций нашли свое практическое применение в различных областях знаний, связанных с обработкой и интерпретацией больших объемов данных: при реконструкции синограмм в медицине, биологии, геологии; при решении задач диагностики плазмы, кристаллографии^[14], дефектоскопии и микроскопии, параллельных вычислений и распределенного анализа^[15] и др.

1.2 Математическое описание алгоритма

1.2.1 Математическое описание

Запишем совместную (в общем случае - переопределенную) систему линейных алгебраических уравнений:

[math]Au=f,\, A\in \mathbb{R}^{m\times n},\, u\in \mathbb{R}^{n},\, f\in \mathbb{R}^{m},[/math]

тогда итерации согласно алгоритму определяются последовательностью

[math]u^{k+1}=u^{k}+\lambda_{k}\frac{f_{j\left(k\right)}-\left(a_{j\left(k\right)},u^{k}\right)}{\Vert a_{j\left(k\right)}\Vert^{2}}a_{j\left(k\right)},[/math]

где

[math]A=\left(a_{1},\, a_{2},\,\ldots\,,\, a_{m}\right)^{T}\in \mathbb{R}^{m\times n},\, u\in \mathbb{R}^{n},\, f=\left(f_{1},\, f_{2},\,\ldots\,,\, f_{m}\right)\in \mathbb{R}^{m},[/math]

а последовательность проектирования в классическом варианте алгоритма определяется как

[math]j\left(k\right)=\left(k\ \bmod\ m\right)+1,\, k=\overline{0,\infty},[/math]

где [math]k[/math] - номер итерации, [math]u^{0}[/math] - начальное приближение, а [math]\lambda_{k} \in \left(0,2\right)[/math] - параметр релаксации (Рисунки 1 и 2 построены для случая, когда [math]\lambda_{k}=1,\,\forall k[/math]).

1.2.2 Геометрическая интерпретация

Метод алгебраической реконструкции достаточно просто формулируется с использованием его геометрической интерпретации. Рассмотрим совместную систему уравнений с двумя неизвестными:

[math]A=\begin{pmatrix}3 & 2\\
2 & 3
\end{pmatrix},\text{ }f=\begin{pmatrix}1\\
2
\end{pmatrix}[/math]

и несовместную переопределенную, полученную добавлением к этой системе дополнительного уравнения [math]-u_{1}+u_{2}=1.5[/math]. Каждое уравнение из системы уравнений можно интерпретировать как гиперплоскость в пространстве координат, а решение совместной системы - как точку пересечения всех этих гиперплоскостей. При использовании метода алгебраической реконструкции для решения системы линейных алгебраических уравнений поиск приближенного решения осуществляется в направлениях, перпендикулярных к гиперплоскостям, каждая из которых описывается отдельным уравнением решаемой системы (см. Рис.1).

Рисунок 1. К геометрической интерпретации метода алгебраической реконструкции: слева - случай совместной системы, справа - несовместной системы

Интерактивная диаграмма для совместного случая подготовлена авторами статьи здесь (Notes about S. Kaczmarz algorithm - consistent case), а для несовместного случая здесь (Notes about S. Kaczmarz algorithm - inconsistent case). Данные диаграммы позволяют регулировать параметры системы уравнений и начальное приближение, что позволяет наглядно демонстрировать сходимость итерационной последовательности.

Описанные результаты наводят на очевидную идею - проекционная последовательность влияет на скорость сходимости метода ортогональных проекций. Обратите внимание на Рис.2.

Рисунок 2. Выбор последовательности проектирования в алгоритме С. Качмажа: К выбору последовательности ортогональных проекций: слева - последовательность проектирования начинается от первого к последнему уравнению и в обратном порядке - от последнего к предпоследнему (symART^[16]^[17]^[18]), а справа - последовательность строится циклическим способом (от последнего к первому уравнению и т.д.)

Как показано на Рис. 2 - итерационная последовательность, в случае Рис. 2 (справа) сходится вообще за конечное число итераций.

1.3 Вычислительное ядро алгоритма

Вычислительное ядро алгоритма Качмажа состоит из двух частей:

предварительный расчет норм всех строк матрицы и задание начального приближения к решению,
серия итераций, на каждой из которых приближенное решение пересчитывается с использованием одной операций SAXPY (scalar multiplication and vector addition) ^[19] и одного скалярного произведения векторов.

Часть предварительного расчета необходима для оптимизации вычислений внутри итерационного процесса, для этого требуется выделение дополнительной памяти ([math]m[/math] числовых ячеек памяти) для хранения норм строк матрицы.

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

достижение необходимой точности;
достижение максимального числа итераций.

Оценку сходимости последовательности можно проводить, например, по невязке (на что будут потрачены дополнительные вычислительные ресурсы), или контролируя норму разности приближенных решений между итерациями. Количество сложений, умножений и сравнений прямо пропорционально количеству итераций.

1.4 Макроструктура алгоритма

Опуская предварительный расчет, базовой макрооперацией алгоритма является операция SAXPY. В подготовительной части алгоритма рассчитываются нормы строк матрицы (заранее могут быть рассчитаны и номера проекций). Далее последовательно в цикле производится запись (перезапись) предыдущего приближения решения, а далее путем обращения к базовым макрооперациям алгоритма: SAXPY и векторного произведения, вычисляется следующее приближение решения. Выход из цикла осуществляется путем проверки достижения заданной точности (сравнение отклонения двух полученных приближений). Исходя из вышеописанного, можно сказать, что параллелизация алгоритма ограничена его итерационной природой (каждое следующее приближение зависит от предыдущего), а эффективность параллельных реализаций операции SAXPY и скалярного произведения увеличивается с увеличением размерности задачи (фактически, с увеличением количества неизвестных в системе уравнений [math]n[/math]).

1.5 Схема реализации последовательного алгоритма

Инициализируем вектор решения первым приближением: [math]u = u^{0}[/math]. В качестве первого приближения зачастую используется нулевой вектор.
Рассчитываем нормы строк матрицы системы уравнений: [math]{\Vert A \Vert}[/math].
Вычисляем очередной номер проекции: [math]j\left(k\right)=\left(k\ \bmod\ m\right)+1,\, k=\overline{0,\infty}[/math].
Для целей обеспечения устойчивости, проверяем, что норма соответствующей проекции не близка к нулю и продолжаем вычисления, иначе переходим к новой итерации (пункт 3).
Запоминаем предыдущее приближение решения и вычисляем новое: [math]u^{k+1}=u^{k}+\lambda_{k}\frac{f_{j\left(k\right)}-\left(a_{j\left(k\right)},u^{k}\right)}{\Vert a_{j\left(k\right)}\Vert^{2}}a_{j\left(k\right)}[/math].
Если решение вычислено с заданной точностью, то выходим из цикла; иначе переходим к новой итерации (пункт 3).

1.6 Последовательная сложность алгоритма

Вычислительная сложность алгоритма Качмажа определяется сложностью реализации операции SAXPY - [math]O(n)[/math]. Затраты на предварительный расчет обычно не включают в сложность самого алгоритма. В итоге, общая сложность алгоритма Качмажа зависит от количества итераций и равна [math]k \cdot O(n)[/math], где [math]n[/math] - количество неизвестных в системы линейных алгебраических уравнений. При этом, в определенных случаях, можно сделать верхние оценки необходимого количества итераций для нахождения приближенного решения с заданной точностью.

К важным особенностям необходимо отнести то, что алгоритмическая сложность каждой итерации не зависит от ее номера, а общий объем памяти, необходимый для реализации алгоритма, линейно зависит от [math]n[/math] (количества неизвестных), но не от [math]m[/math] (количества измерений).

1.7 Информационный граф

Исходя из описания алгоритма, можно построить информационный граф алгоритма (см. Рис. 3). В качестве вершин графа выступают следующие операции:

Задание начального приближения: [math]U = Zeros()[/math] (см. пункт 1 раздела 1.5).
Вычисление норм строк матрицы: [math]nrma2 = Norms(A)[/math] (см. пункт 2 раздела 1.5).
Вычисление номера проекции (вычисляется в цикле): [math]j = Index()[/math] (см. пункт 3 раздела 1.5).
Сохранение предыдущего приближения решения (вычисляется в цикле): [math]U_{prev} = U[/math].
Расчет нового приближения решения (вычисляется в цикле): [math]U = Calc(U,j,nrma2)[/math] (см. пункт 5 раздела 1.5).

Граф алгоритма является параметризованным (зависит от входных данных) и не детерминированным (в смысле, что количество вычислений зависит от начальных данных).

Рисунок 3. Информационный граф алгоритма Качмажа.

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

1.8 Ресурс параллелизма алгоритма

Алгоритм Качмажа, как любой итерационный алгоритм, является последовательным и в строковой версии алгоритма ресурс для параллелизма оказывается небольшим - только за счет параллельных реализаций операции SAXPY и скалярного произведения.

Однако, существуют различные блочные версии алгоритма Качмажа, где проектирование осуществляется не на конкретную гиперплоскость, а на их пересечение в рамках выбранного блока. Такие модификации имеют существенно больший ресурс для распараллеливания и связаны с трудоемкими задачами псевдообращения и задачей выбора оптимального размера блока (в общем случае не предполагается, что система разбивается на блоки равных размеров, что еще больше усложняет задачу).

1.9 Входные и выходные данные алгоритма

Входными данными алгоритма являются:

матрица системы уравнений [math]A\in \mathbb{R}^{m\times n}[/math],
вектор правых частей [math]f\in \mathbb{R}^{m}[/math],
требуемая точность [math]\epsilon\gt 0[/math],
а так же максимальное количество итераций [math]K_{max}\in \mathbb{N}[/math].

Объем входных данных составляет [math]m \cdot n + m + 2[/math] ячеек памяти.

Выходными данными алгоритма является вектор решения системы уравнений [math]u\in \mathbb{R}^{n}[/math]. Объем выходных данных составляет [math]n[/math] ячеек памяти.

Важно отметить, что для начала итераций (в случае, когда этап предварительных расчетов опускается) вся система не нужна, достаточно всего одного уравнения, остальные можно подавать в "решатель" последовательно. Таким образом, в задачах обработки результатов измерений нет необходимости запоминать в оперативной памяти все результаты измерений, достаточно результаты каждого измерения подавать в "решатель" и по завершению итерации - забывать их (особенно интересно, в случае, когда провести измерение повторно оказывается дешевле, чем сохранить его результаты). Как считают авторы, данный факт играл решающую роль в выборе алгоритма для решения задач компьютерной томографии, в которых измерения (срезы синограммы) появляются последовательно, а повторение измерения лишь увеличивает лучевую нагрузку на сканируемый объект и само время сканирования (но, объем лучевой нагрузки, естественно, является критичным для биологических объектов).

1.10 Свойства алгоритма

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

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

Алгоритм недетерминирован по структуре вычислительного процесса: число итераций, а значит и число операций, зависят от входных данных. С другой стороны, алгоритм детерминирован в том смысле, что выдаёт всегда одинаковый результат при одинаковых входных данных (естественно, исключая рандомизированные модификации алгоритма Качмажа).

2 Программная реализация алгоритма

Ниже представлен прототип программной реализации алгоритма на языке MatLab. Большое количество вариантов реализации и различных модификаций представлено в этих пакетах ^[20]^[21].

function [ u ] = ART(A, f, maxit, tol)
%ART1 - алгоритм алгебраической реконструкции (С. Качмажа) с циклическим обходом
%   A, f - исходные данные задачи
%   tol - чувствительность правила останова
	[m,n] = size(A); u = zeros(n,1);
	% Нормы строк матрицы СЛАУ рассчитываются заранее
	nrma2 = full(abs(sum(A.*A',1)));
	for k = 1:1:maxit     
		j = mod(k-1,m) + 1;
		if (nrma2(j) <= 0), continue, end;
		u_prev = u;
		u = u + (f(j) - (A(j,:)*u))/nrma2(j) * A(j,:)';
		if (norm(u_prev-u) <= tol), break, end;
	end;
end;

2.1 Особенности реализации последовательного алгоритма

В качестве критерия завершения итераций используют наиболее популярные правила останова (но не всегда эффективные): критерий останова по невязке [math]\left\Vert A\cdot u_{k}-f\right\Vert \lt \epsilon[/math] или [math]\left\Vert u_{k-1}-u_{k}\right\Vert \lt \epsilon[/math] (как в представленной реализации). Последний критерий может слишком рано завершать итерации, в виду того, что часто алгоритм имеет медленную скорость сходимости.

2.2 Возможные способы и особенности параллельной реализации алгоритма

Векторность операции SAXPY может обеспечить хорошую пространственную локальность алгоритму Качмажа.

Важно, при реализации учитывать структуру хранения двумерных массивов в оперативной памяти для выбранной архитектуры. К примеру, в языках С/С++ двумерный массив (матрица) из трех строк и четырех колонок будет сохранен в памяти в следующей последовательности: (1,1),(1,2),(1,3),(1,4),(2,1),(2,2),(2,3),(2,4),(3,1),(3,2),(3,3),(3,4) - элементы матрицы выписываются последовательно, по строкам.

Однако, при реализации на языке FORTRAN эта же матрица будет сохранена в следующей последовательности: (1,1),(2,1),(3,1),(1,2),(2,2),(3,2),(1,3),(2,3),(3,3),(1,4),(2,4),(3,4) - элементы матрицы выписываются последовательно, по столбцам ^[22]. При реализации на FORTRAN-подобных языках, матрицу [math]A[/math] из входных данных следует хранить в оперативной памяти в транспонированном виде, таким образом существенно возрастает пространственная локальность векторных операций, в которых участвуют строки матрицы системы линейных алгебраических уравнений.

Как уже упоминалось выше, параллелизация алгоритма Качмажа возможна за счет распараллеливания реализаций макроопераций алгоритма. Основными операциями, требующими распараллеливания, являются:

операция скалярного произведения векторов;
операция SAXPY;
вычисление евклидовой нормы вектора.

При реализации параллельной версии алгоритма на системах с общей памятью параллельную программу можно разделить на следующие этапы.

Распределение памяти, загрузка исходных данных, уточнение параметров параллелизации.
Предрасчет проекционной последовательности (для рандомизированной реализации, оказывается эффективнее сгенерировать случайную последовательность из [math]j(k)[/math] заранее).
Вычисление приближения решения.
Сохранение результатов.

Для параллельных реализаций алгоритма достаточно иметь соответсвующие реализации для BLAS (Basic Linear Algebra Subprograms) операций. К примеру, одна из модификаций алгоритма Качмажа для решения задачи регуляризации А.Н. Тихонова представлена в ^[23] (описание данной модификации выходит за рамки содержания данной публикации). В этой работе использовались функции

cblas_ddoti - скалярное произведение векторов,
cblas_daxpyi - реализация операции SAXPY,

из библиотеки базовых математических операций Intel MKL. В статье отмечается, что за счет распараллеливания макроопераций представленная реализация позволила уменьшить время решения задачи от 2 до 4 раз по сравнению с последовательной версией.

Реализации базовых операций линейной алгебры существуют и для графических процессоров, например, библиотека cuBLAS для архитектуры параллельных вычислений от NVIDIA. Как отмечено в работе ^[24], параллельная реализация алгоритма Качмажа на CUDA архитектуре оказывается неэффективной за счет операций первичной инициализации и копирования между различными ОЗУ.

Наиболее перспективными архитектурами для реализации алгоритма Качмажа, по мнению авторов, являются ARM архитектуры (Advanced RISC Machine, Acorn RISC Machine, усовершенствованная RISC-машина). Например, существенное ускорение можно получить в реализациях для процессоров, построенных по SIMD (single instruction, multiple data — одиночный поток команд, множественный поток данных) принципу. Стоит отметить технологию ARM NEON, которая может быть применена при разработке приложений для iPhone устройств, в которых, по причине ограниченного заряда батареи и в целом вычислительных ресурсов мобильного устройства, остро стоят задачи оптимизации вычислений при обработке снимков с камеры телефона или видео-ряда в режиме реального времени.

2.3 Результаты прогонов

2.4 Выводы для классов архитектур

Как было сказано выше, ресурсом параллелизма в алгоритме являются векторные операции: SAXPY и скалярное произведение векторов. Данные макрооперации имеют высокую степень параллелизма и могут быть реализованы для различных классов архитектур. Однако факт того, что вычисления выполнятюся многократно в цикле, требующем полного вектора решения на этапе проверки достижения необходимой точности, говорит о целесообразности применения архитектур с общей памятью.

Так же стоит отметить, что графические процессоры хорошо зарекомендовали себя для реализации матричных и векторных операций, в частности той же операции SAXPY, например, в ^[25]. Но, затраты, связанные с операциями инициализации, обращения и копирования между основным ОЗУ и ОЗУ графического устройства, нивелируют какое-либо преимущество, получаемое за счет ускорения времени выполнения одной итерации алгоритма.

Авторы считают, что векторные процессоры, построенные по SIMD принципам, находятся вне конкуренции для реализации алгоритма Качмажа, например, архитектура ARM NEON.

3 ЧАСТЬ. Дополнительная информация

3.1 Классификация проекционных последовательностей

Под проекционной последовательностью понимается порядок обхода уравнений системы в процессе итераций алгоритма Качмажа. А рассматриваемая последовательность задается с использованием правила [math]j\left(k\right)[/math], где [math]k\in \mathbb{N}_{0}[/math], которое позволяет вычислить номер используемого уравнения на итерации с номером [math]k[/math]. Правило [math]j\left(k\right)[/math] следует рассматривать как сюръективное отображение множества [math]\mathbb{N}_{0}[/math] на множество [math]\left\{ 1,2,\ldots,m\right\}[/math], где [math]m[/math] - количество уравнений в системе линейных алгебраических уравнений [math]Au=f,\ A\in \mathbb{R}^{m\times n},\ u\in \mathbb{R}^{n},\ f\in \mathbb{R}^{m}[/math].

Отображения такого вида могут быть как детерминированными, так и случайными, рандомизированными. В последнем случае [math]j\left(k\right)[/math] является дискретной случайной величиной с заданным законом распределения. Например, к детерминированным отображениям необходимо отнести:

правило циклического обхода [math]j\left(k\right)=\left(k\ \bmod\ m\right)+1,\, k=\overline{0,\infty}[/math];
метод четной перестановки - эвристический метод для задач компьютерной томограции из ^[26];
симметричное правило (symmetric ART), в котором последовательность проектирования на гиперплоскости идет в начале в прямом порядке от первого уравнения к последнему, а затем в обратном порядке - от последнего - к первому.

Существенно важной является классификация рассматриваемых правил по их зависимости от данных. Например, результат из ^[13], где [math]j\left(k\right)[/math] является дискретной случайной величиной с заданным законом распределения [math]P\left(j\left(k\right)=i\right)=\frac{\left\Vert a_{i}\right\Vert _{2}^{2}}{\left\Vert A\right\Vert _{F}^{2}}[/math], где [math]\left\Vert \cdot\right\Vert _{F}[/math] - матричная норма Фробениуса и [math]\left\Vert a_{i}\right\Vert[/math] - евклидова норма строки матрицы системы уравнений с номером [math]i=\overline{1,m}[/math]. Этот результат качественно отличается от всех известных рандомизированных способов наличием зависимости от матрицы системы - от исходных данных задачи.

Рисунок 3. Классификация правил для задания проекционной последовательности в алгоритме C. Качмажа

3.2 Рандомизированная версия алгоритма и скорость сходимости

При решении системы линейных алгебраических уравнений сверхбольшой размерности использование всех уравнений системы может оказаться избыточным. Более того, для учета всех уравнений системы, полученных в результате некоторых измерений, может не быть технической возможности, что обусловлено ограниченными ресурсами ЭВМ. Существуют подходы, связанные с обработкой данных решаемой задачи, цель которых состоит в отборе только информативных результатов измерений.

Например, при решении задачи идентификации оптических искажающих систем ^[27] следует использовать не весь объем измерений, а только тот, что выполнен для, в определенном смысле, информативных фрагментов изображений. Формально, это соответствует удалению части уравнений из системы, которые оказываются линейно зависимыми и отрицательным образом влияют на обусловленность задачи ^[28], что, в свою очередь, приводит к крайне медленной скорости сходимости известных итерационных алгоритмов.

Таким образом, основная цель рандомизации итерационных проекционных алгоритмов решения системы линейных алгебраических уравнений состоит в том, чтобы обеспечить скорость сходимости (а значит и вычислительную сложность), независимую от количества уравнений в системе (количества измерений) ^[29], но при решении не использовать систему целиком, а только отдельные уравнения, выбранные определенным случайным образом.

Наиболее простым и очевидным способом случайную последовательность [math]j\left(k\right)[/math] определяют с использованием дискретного равномерного закона распределения [math]P\left(j\left(k\right)=i\right)=\frac{1}{m}[/math].

Таким образом, в процессе итераций все уравнения системы будут участвовать в решении задачи примерно одинаковое число раз. Данный способ задания [math]j\left(k\right)[/math] обширно используют в качестве базового для демонстрации преимущества новых подходов в работах ^[13] ^[30], посвященных рандомизации итерационных алгоритмов. В указанных работах авторы произвели настоящую революцию в рандомизации алгоритма Качмажа, указав закон распределения, определяющий последовательность [math]j\left(k\right)[/math] таким образом, что итерационная последовательность имеет экспоненциальную скорость сходимости. Их результат формулируется в виде следующей важной теоремы:

Пусть [math]Au=f\ A\in \mathbb{R}^{m\times n},\ u\in \mathbb{R}^{n},\ f\in \mathbb{R}^{m}[/math] - совместная переопределенная [math]\left(m\geq n\right)[/math] система линейных алгебраических уравнений, а [math]u_{0}[/math] - некоторое начальное приближение к ее решению, и для рекуррентной формулы алгоритма Качмажа значения последовательности [math]j\left(k\right)[/math] определяются с использованием закона распределения: [math] P\left(j\left(k\right)=i\right)=\frac{\left\Vert a_{i}\right\Vert _{2}^{2}}{\left\Vert A\right\Vert _{F}^{2}}, [/math] где [math]\left\Vert \cdot\right\Vert _{F}[/math] - матричная норма Фробениуса.

Тогда последовательность приближенных решений в алгоритме Качмажа сходится к решению системы уравнений со средней ошибкой, которая на каждой итерации оценивается с использованием неравенств

[math] M\left[\left\Vert u^{k}-u^{*}\right\Vert _{2}^{2}\right]\leq\left(1-\kappa\left(A\right)^{-2}\right)^{k}\left\Vert x^{0}-x^{*}\right\Vert _{2}^{2}, [/math]

[math] \left(1-2k\cdot\kappa\left(A\right)^{-2}\right)\left\Vert x^{0}-x^{*}\right\Vert _{2}\leq M\left[\left\Vert u^{k}-u^{*}\right\Vert _{2}\right]. [/math]

В приведенных неравенствах величина [math]\kappa\left(A\right)=\left\Vert A\right\Vert _{F}\left\Vert A^{+}\right\Vert _{2}[/math] является одной из оценок [math]1\leq\frac{\kappa\left(A\right)}{\sqrt{n}}\leq k\left(A\right)[/math] спектрального числа обусловленности матрицы системы [math]k\left(A\right)=\left\Vert A\right\Vert _{2}\left\Vert A^{+}\right\Vert _{2}[/math], предложенной в ^[31] (здесь [math]A^{+}[/math] - псевдообратная матрица, в то же время левая обратная матрица для [math]A[/math]).

До публикации работ Штромера (Strohmer, T.) не существовало оценок скорости сходимости алгоритма Качмажа в терминах числа обусловленности матрицы системы уравнений. Данная теорема также наиболее точно демонстрирует цель рандомизации - показать, что скорость сходимости не зависит от [math]m[/math] - количества уравнений в системе, но зависит от порядка системы [math]n[/math].

3.3 Алгоритм Качмажа с четной перестановкой

Данная эвристика оказывается, пожалуй, одной из самых изящных и определяет отображение [math]j\left(k\right)[/math] следующим образом. В работе ^[32] предлагается алгоритм MLS (multilevel scheme) для переупорядочивания строк в матрице таким образом, чтобы каждое следующее проецирование в итерационной последовательности осуществлялось на гиперплоскость, наиболее ортогональную к текущей (только для матриц в задачах томографии). Допустим, что имеется 8 проекций (в матрице системы - всего восемь строк), снятых под углами в 0, 22.5, 45 ,67.5, 90, 112.5, 135 и 157.5 градусов, тогда период последовательности [math]j\left(k\right)[/math] выглядит как: [math]{0,4,2,6,1,5,3,7}[/math].

Когда количество проекций (строк в системе) является степенью двойки, то [math]j\left(k\right)[/math] совпадет с последовательностью, которая используется для реализации алгоритма «на месте» одномерного быстрого преобразования Фурье, например в ^[33].

Формально, данная модификация эквивалентна исходному алгоритму с циклическим правилом обхода, который применяется к системе уравнений. [math]D \cdot A \cdot u=D \cdot f[/math] , где [math]D\in R^{m\times m}[/math] суть - перестановочная матрица.

В чем же состоит изящность данного способа задания отображения [math]j\left(k\right)[/math]? Правило четной перестановки, было специально разработано для решения задачи восстановления данных по измерениям, выполненным с использованием компьютерного томографа с параллельной схемой сканирования. Правило построено его авторами из интуитивных соображений о структуре самих данных, а не их значений.

4 Литература

↑ Бахвалов Н. С., Жидков Н. П., Кобельков Г. Г. Численные методы. — 8-е изд. — М.: Лаборатория Базовых Знаний, 2000.
↑ Natterer, Frank (2001), "V.3 Kaczmarz's method", The Mathematics of Computerized Tomography, Classics in Applied Mathematics 32, SIAM, p. 128, ISBN 9780898714937.
↑ Kaczmarz, S. Angenäherte Auflösung von Systemen linearer Gleichungen [Текст] / S. Kaczmarz // Bulletin International de l'Académie Polonaise des Sciences et des Lettres / Classe des Sciences Mathématiques et Naturelles. Série A, Sciences Mathématiques. — 1937. — Vol. 35. — С. 355-357.
↑ Hounsfield, G. N. A Discussion on Recent Developments in Medical Endoscopy and Related Fields [Текст] / G. N. Hounsfield // Proceedings of the Royal Society of London. Series B, Biological Sciences. — 1977. — Vol. 195, N. 1119. — С. 281-289.
↑ Gordon, R. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography [Текст] / R. Gordon, R. Bender, G. T. Herman // Journal of theoretical biology. — 1970. — Vol. 29, Issue 3. — С. 471–481.
↑ Micke, A. The Mathematics of Computerized Tomography [Текст] / A. Micke, F. Natterer // ZAMM - Zeitschrift für Angewandte Mathematik und Mechanik (Journal of Applied Mathematics and Mechanics). — 1987. — Vol. 67, Issue 11. — С. 580. —ISBN 3-519-02103-X and 0-471-90959-9.
↑ Herman, G. T. Fundamentals of computerized tomography: Image reconstruction from projection [Текст] : 2nd ed. / G. T. Herman. — Springer, 2009. — 300 с.
↑ Sezan, K. M. Applications of convex projection theory to image recovery in tomography and related areas [Текст] / K. M. Sezan, H. Stark // Image Recovery: Theory and Application / H. Stark, editor. — Academic Press, 1987. — С. 415-462.
↑ Cegielski, A. Bibliography on the Kaczmarz method [Электронный ресурс] : Faculty of Mathematics, Computer Science and Econometrics University of Zielona Gora ; Poland, March, 28th, 2010 / Режим доступа: http://www.uz.zgora.pl/~acegiels/Publikacje-Kaczmarz.pdf
↑ Ильин, В. П. Об итерационном методе Качмажа и его обобщениях [Текст] / В. П. Ильин // Сибирский журнал индустриальной математики. — 2006. — Т. 9, вып. 3. — С. 39–49.
↑ Г. П. Васильченко, А. А. Светлаков, “Проекционный алгоритм решения систем линейных алгебраических уравнений большой размерности”, Ж. вычисл. матем. и матем. физ., 20:1 (1980), 3–10.
↑ Л. Г. Гурин, Б. Т. Поляк, Э. В. Райк, “Методы проекций для отыскания общей точки выпуклых множеств”, Ж. вычисл. матем. и матем. физ., 7:6 (1967), 1211–1228
↑ ^13,0 ^13,1 ^13,2 Strohmer, T. A Randomized Kaczmarz Algorithm with Exponential Convergence [Текст] / T. Strohmer, R. Vershynin // Journal of Fourier Analysis and Applications. — 2009. — Vol. 15, Issue 2. — С. 262-278.
↑ Marks, L. D. A feasible set approach to the crystallographic phase problem [Текст] / L. D. Marks, W. Sinkler, E. Landree // Acta Crystallographica Section A. — 1987. — Vol. 55, Issue 4. — С. 601-612.
↑ JElble, J. M. GPU computing with Kaczmarz’s and other iterative algorithms for linear systems [Текст] / J. M. Elble, N. V. Sahinidis, P. Vouzis // Parallel Computing. — 2010. — Vol. 36, Issues 5-6. — С. 215–231.
↑ Bjorck, A. Accelerated projection methods for computing pseudoinverse solutions of systems of linear equations [Текст] / A. Bjorck, T. Elfving // BIT Numerical Mathematics. — 1979. — Vol. 19, Issue 2. — С.145-163.
↑ Nikazad, T. Algebraic Reconstruction Methods [Текст] : PhD. Thesis, Linköping University / T. Nikazad. — Linköping, 2008.
↑ Kamath, G. Component-Average Based Distributed Seismic Tomography in Sensor Networks [Текст] / G. Kamath, L. Shi, W.-Z. Song // IEEE International Conference on Distributed Computing in Sensor Systems. — 2013. — С. 88-95.
↑ Mark Harris. Six Ways to SAXPY https://devblogs.nvidia.com/parallelforall/six-ways-saxpy/
↑ Hansen, P. C. AIR Tools – A MATLAB package of algebraic iterative reconstruction methods [Электронный ресурс] / Режим доступа: http://www2.compute.dtu.dk/~pcha/AIRtools/
↑ Ivanov, A. A. Regularization Kaczmarz Tools Version 1.4 for Matlab, Matlabcentral Fileexchange [Электронный ресурс] / Режим доступа: http://www.mathworks.com/matlabcentral/fileexchange/43791
↑ How FORTRAN Stores Two-Dimensional Arrays in Memory, https://support.microsoft.com/en-us/kb/27780
↑ А.И. Жданов, Ю.В. Сидоров. Параллельная реализация рандомизированного регуляризованного алгоритма Качмажа / Компьютерная оптика. - 2015. - Т.39(4). - С. 536-541
↑ Elble J. M., Sahinidis N. V., Vouzis P. GPU computing with Kaczmarz’s and other iterative algorithms for linear systems //Parallel computing. – 2010. – Т. 36. – №. 5. – С. 215-231.
↑ А.Ю. Суравикин, В.В. Коробицын Оценка производительности связывания CUDA с графическим API на примере задачи SAXPY / Математические структуры и моделирование. - 2010ю - Т.21. - С. 96-103.
↑ Guan H. A projection access order for speedy convergence of ART (algebraic reconstruction technique): a multilevel scheme for computed tomography [Текст] / Guan H., R. Gordon // Physics in Medicine and Biology. — 1994. — Vol. 39, N. 11. — С. 2005-2022.
↑ Фурсов, В.А. Идентификация оптических искажающих систем с отбором информативных фрагментов изображений [Текст] / В.А. Фурсов // Компьютерная оптика. — 1995. —Вып. 15, ч.1. — С. 79–89.
↑ Фурсов, В.А. Введение в идентификацию по малому числу наблюдений. [Текст] / В.А. Фурсов. — М. : Изд-во МАИ, 1991. —36 с.
↑ Sabelfeld, K Stochastic iterative projection methods for large linear systems [Текст] / K. Sabelfeld, N. Loshchina // Monte Carlo Methods and Applications. — 2010. — Vol. 16, Issue 3-4. — С. 343-359.
↑ Strohmer, T. A randomized solver for linear systems with exponential convergence [Текст] / T. Strohmer, R. Vershynin // Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques. — Springer; 2006. С. 499-507.
↑ Demmel, J. The probability that a numerical analysis problem is difficult [Текст] / J. Demmel // Mathematics of Computation. . — 1988. — Vol. 50, Issue 182. — С. 449–480.
↑ Guan H. A projection access order for speedy convergence of ART (algebraic reconstruction technique): a multilevel scheme for computed tomography [Текст] / Guan H., R. Gordon // Physics in Medicine and Biology. — 1994. — Vol. 39, N. 11. — С. 2005-2022.
↑ Cooley, J.W. An algorithm for the machine calculation of complex Fourier series [Текст] / J.W. Cooley, J. W. Tukey // Mathematics of Computation. — 1965. — Vol. 19, N. 90. — С. 297–301.

[1] Бахвалов Н. С., Жидков Н. П., Кобельков Г. Г. Численные методы. — 8-е изд. — М.: Лаборатория Базовых Знаний, 2000.

[2] Natterer, Frank (2001), "V.3 Kaczmarz's method", The Mathematics of Computerized Tomography, Classics in Applied Mathematics 32, SIAM, p. 128, ISBN 9780898714937.

[3] Kaczmarz, S. Angenäherte Auflösung von Systemen linearer Gleichungen [Текст] / S. Kaczmarz // Bulletin International de l'Académie Polonaise des Sciences et des Lettres / Classe des Sciences Mathématiques et Naturelles. Série A, Sciences Mathématiques. — 1937. — Vol. 35. — С. 355-357.

[4] Hounsfield, G. N. A Discussion on Recent Developments in Medical Endoscopy and Related Fields [Текст] / G. N. Hounsfield // Proceedings of the Royal Society of London. Series B, Biological Sciences. — 1977. — Vol. 195, N. 1119. — С. 281-289.

[5] Gordon, R. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography [Текст] / R. Gordon, R. Bender, G. T. Herman // Journal of theoretical biology. — 1970. — Vol. 29, Issue 3. — С. 471–481.

[6] Micke, A. The Mathematics of Computerized Tomography [Текст] / A. Micke, F. Natterer // ZAMM - Zeitschrift für Angewandte Mathematik und Mechanik (Journal of Applied Mathematics and Mechanics). — 1987. — Vol. 67, Issue 11. — С. 580. —ISBN 3-519-02103-X and 0-471-90959-9.

[7] Herman, G. T. Fundamentals of computerized tomography: Image reconstruction from projection [Текст] : 2nd ed. / G. T. Herman. — Springer, 2009. — 300 с.

[8] Sezan, K. M. Applications of convex projection theory to image recovery in tomography and related areas [Текст] / K. M. Sezan, H. Stark // Image Recovery: Theory and Application / H. Stark, editor. — Academic Press, 1987. — С. 415-462.

[9] Cegielski, A. Bibliography on the Kaczmarz method [Электронный ресурс] : Faculty of Mathematics, Computer Science and Econometrics University of Zielona Gora ; Poland, March, 28th, 2010 / Режим доступа: http://www.uz.zgora.pl/~acegiels/Publikacje-Kaczmarz.pdf

[10] Ильин, В. П. Об итерационном методе Качмажа и его обобщениях [Текст] / В. П. Ильин // Сибирский журнал индустриальной математики. — 2006. — Т. 9, вып. 3. — С. 39–49.

[11] Г. П. Васильченко, А. А. Светлаков, “Проекционный алгоритм решения систем линейных алгебраических уравнений большой размерности”, Ж. вычисл. матем. и матем. физ., 20:1 (1980), 3–10.

[12] Л. Г. Гурин, Б. Т. Поляк, Э. В. Райк, “Методы проекций для отыскания общей точки выпуклых множеств”, Ж. вычисл. матем. и матем. физ., 7:6 (1967), 1211–1228

[strohmer-13] 13,0 ^13,1 ^13,2 Strohmer, T. A Randomized Kaczmarz Algorithm with Exponential Convergence [Текст] / T. Strohmer, R. Vershynin // Journal of Fourier Analysis and Applications. — 2009. — Vol. 15, Issue 2. — С. 262-278.

[14] Marks, L. D. A feasible set approach to the crystallographic phase problem [Текст] / L. D. Marks, W. Sinkler, E. Landree // Acta Crystallographica Section A. — 1987. — Vol. 55, Issue 4. — С. 601-612.

[15] JElble, J. M. GPU computing with Kaczmarz’s and other iterative algorithms for linear systems [Текст] / J. M. Elble, N. V. Sahinidis, P. Vouzis // Parallel Computing. — 2010. — Vol. 36, Issues 5-6. — С. 215–231.

[16] Bjorck, A. Accelerated projection methods for computing pseudoinverse solutions of systems of linear equations [Текст] / A. Bjorck, T. Elfving // BIT Numerical Mathematics. — 1979. — Vol. 19, Issue 2. — С.145-163.

[17] Nikazad, T. Algebraic Reconstruction Methods [Текст] : PhD. Thesis, Linköping University / T. Nikazad. — Linköping, 2008.

[18] Kamath, G. Component-Average Based Distributed Seismic Tomography in Sensor Networks [Текст] / G. Kamath, L. Shi, W.-Z. Song // IEEE International Conference on Distributed Computing in Sensor Systems. — 2013. — С. 88-95.

[19] Mark Harris. Six Ways to SAXPY https://devblogs.nvidia.com/parallelforall/six-ways-saxpy/

[airtools-20] Hansen, P. C. AIR Tools – A MATLAB package of algebraic iterative reconstruction methods [Электронный ресурс] / Режим доступа: http://www2.compute.dtu.dk/~pcha/AIRtools/

[regkaz-21] Ivanov, A. A. Regularization Kaczmarz Tools Version 1.4 for Matlab, Matlabcentral Fileexchange [Электронный ресурс] / Режим доступа: http://www.mathworks.com/matlabcentral/fileexchange/43791

[22] How FORTRAN Stores Two-Dimensional Arrays in Memory, https://support.microsoft.com/en-us/kb/27780

[Sidorov_Jdanov-23] А.И. Жданов, Ю.В. Сидоров. Параллельная реализация рандомизированного регуляризованного алгоритма Качмажа / Компьютерная оптика. - 2015. - Т.39(4). - С. 536-541

[24] Elble J. M., Sahinidis N. V., Vouzis P. GPU computing with Kaczmarz’s and other iterative algorithms for linear systems //Parallel computing. – 2010. – Т. 36. – №. 5. – С. 215-231.

[saxpy1-25] А.Ю. Суравикин, В.В. Коробицын Оценка производительности связывания CUDA с графическим API на примере задачи SAXPY / Математические структуры и моделирование. - 2010ю - Т.21. - С. 96-103.

[26] Guan H. A projection access order for speedy convergence of ART (algebraic reconstruction technique): a multilevel scheme for computed tomography [Текст] / Guan H., R. Gordon // Physics in Medicine and Biology. — 1994. — Vol. 39, N. 11. — С. 2005-2022.

[27] Фурсов, В.А. Идентификация оптических искажающих систем с отбором информативных фрагментов изображений [Текст] / В.А. Фурсов // Компьютерная оптика. — 1995. —Вып. 15, ч.1. — С. 79–89.

[28] Фурсов, В.А. Введение в идентификацию по малому числу наблюдений. [Текст] / В.А. Фурсов. — М. : Изд-во МАИ, 1991. —36 с.

[29] Sabelfeld, K Stochastic iterative projection methods for large linear systems [Текст] / K. Sabelfeld, N. Loshchina // Monte Carlo Methods and Applications. — 2010. — Vol. 16, Issue 3-4. — С. 343-359.

[30] Strohmer, T. A randomized solver for linear systems with exponential convergence [Текст] / T. Strohmer, R. Vershynin // Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques. — Springer; 2006. С. 499-507.

[31] Demmel, J. The probability that a numerical analysis problem is difficult [Текст] / J. Demmel // Mathematics of Computation. . — 1988. — Vol. 50, Issue 182. — С. 449–480.

[32] Guan H. A projection access order for speedy convergence of ART (algebraic reconstruction technique): a multilevel scheme for computed tomography [Текст] / Guan H., R. Gordon // Physics in Medicine and Biology. — 1994. — Vol. 39, N. 11. — С. 2005-2022.

[33] Cooley, J.W. An algorithm for the machine calculation of complex Fourier series [Текст] / J.W. Cooley, J. W. Tukey // Mathematics of Computation. — 1965. — Vol. 19, N. 90. — С. 297–301.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]