Участник:Морозов Егор Павлвочи/Метод решения недетерминированных обратных задач идентификации входных параметров с BS регуляризацией

Метод решения недетерминированных обратных задач идентификации входных параметров с BS регуляризацией
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ O(nm+m^2) }[/math]
Объём входных данных	[math]\displaystyle{ 2m+n+1 }[/math]
Объём выходных данных	[math]\displaystyle{ 4n+1 }[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]\displaystyle{ 5 }[/math]
Ширина ярусно-параллельной формы	[math]\displaystyle{ O(n) }[/math]

Основные авторы описания: Е.П.Морозов

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

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

Недетерминированные обратные задачи идентификации входных параметров относятся к классу задач, в которых искомое решение зависит от входных данных с некоторым разбросом и характеризуется отсутствием однозначной связи между наблюдаемыми данными и решением. Формально постановка такой задачи выглядит следующим образом. Пусть имеется некоторый вектор неизвестных величин [math]\displaystyle{ \vec{x} }[/math], который мы хотим восстановить по наблюдениям [math]\displaystyle{ \vec{\theta} }[/math]. Между ними существует связь вида: \begin{equation} \vec{\theta} = F(\vec{x}, \vec{\varepsilon}), \end{equation} где F — оператор прямого преобразования, [math]\displaystyle{ \vec{\varepsilon} }[/math] — случайная ошибка измерения или воздействия среды. Задача восстановления входных параметров [math]\displaystyle{ \vec{x} }[/math] называется обратной задачей идентификации входных параметров. Недетерминированные обратные задачи зачастую являются в том числе некорректными. Чтобы понять какие задачи лежат в данном классе, рассмотрим его антоним - корректные задачи. Согласно определению, которое было дано Тихоновым А.Н. :

Задача определения решения [math]\displaystyle{ \vec{\theta} }[/math] из пространства F по «исходным данным» и [math]\displaystyle{ \vec{x} }[/math] из пространства U называется корректно поставленной на паре метрических пространств (F, U), если удовлетворяются требования (условия):

• для всякого элемента [math]\displaystyle{ x \in U }[/math] существует решение [math]\displaystyle{ {\theta} }[/math] из пространства F;

• решение определяется однозначно;

• задача устойчива на пространствах (F, U).

В математической литературе длительное время существовала точка зрения, согласно которой всякая математическая задача должна удовлетворять этим требованиям . Задачи, не удовлетворяющие перечисленным требованиям, называются некорректно поставленными. Следует отметить, что определение некорректно поставленных задач относится к данной паре метрических пространств (F, U), так как в других метриках та же задача может быть корректно поставленной. В рассматриваемом классе задач нарушается устойчивость и однозначность решения. Согласно А.Н. Тихонову, для таких задач существуют квазирешения. Однако чтобы перейти к решению, а поняв его структуру, — к вычислению обратной задачи, необходимо восстановить приблизительную зависимость между входными и выходными параметрами.

Для этого рассмотрим B-сплайны и NURBS. Их полное геометрическое описание и вывод основных свойств были представлены в «NURBS book» . Помимо их основного применения — создания и сглаживания поверхностей и сложных форм на основе реперных точек, — мы можем применить их для аппроксимации исходных данных и получения системы алгебраических уравнений.

\begin{equation} B\vec{x} \approx \vec{\theta}, \end{equation}

Данные системы функциональных уравнений и связанные с ними обратные задачи можно решать различными способами. Пример решения подобных задач для уравнений математической физики представлен в книге , где они называются обратными задачами идентификации правой части. Однако представленный алгоритм ориентирован не на применение конечно-разностных схем, не на регуляризацию по Тихонову или иные традиционные методы, а на использование нового регуляризатора, основанного на B-сплайнах. При аппроксимации наших данных мы вводили B-сплайны как многочлены порядка m, обозначим их как [math]\displaystyle{ Q_m(u) }[/math], где u - наш новый параметр. Для более точного решения нам необходимо, чтобы данный аппроксимативный многочлен удовлетворял условиям альтернанса Чебышева. Тем самым [math]\displaystyle{ Q_m(u) }[/math] будет полиномом наилучшего приближения. Тем самым мы можем за основу взять условия альтернанса \begin{equation} F(\vec{x_i})-Q_m(u)=\alpha(-1)^i||{F- Q_m (u)}||, \end{equation} где [math]\displaystyle{ i=1\ldots m+1 }[/math], [math]\displaystyle{ \alpha=\pm1 }[/math] одновременно для всех i. Для обобщения решения будем рассматривать Гильбертовы пространства и в качестве нормы введём интегральную гильбертову норму. А в качестве [math]\displaystyle{ F(\vec{x_i}) }[/math] рассмотрим наше известное значение для наблюдения, идентификацию начальных параметров для которого мы хотим получить, обозначим его как [math]\displaystyle{ \Theta }[/math], таким образом для решения исходной обратной задачи с учётом аппроксимации B-сплайнами, нам необходимо найти параметр U, который вычисляется из выражения \begin{equation} U = \operatorname*{argmin}_u ((Q_m (u) - \Theta)^2 + \Theta - \int_{0}^{1}Q_m (u)\,du) \end{equation}

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

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

Для математического описания алгоритма рассмотрим рассуждения которые были описаны в пункте 1.1, но акцент сделаем именно не его математической части. Таким образом пусть дан определённый массив наблюдений [math]\displaystyle{ \vec{\theta} = (\theta_1, \theta_2, \ldots, \theta_m) }[/math], где m - количество наблюдений, которые были получены при начальных параметрах [math]\displaystyle{ \vec{x_i} = (x_{1i}, x_{2i}, \ldots, x_{ni}) }[/math], где n - количество исходных параметров для рассматриваемой системы, а также дано целевое значение наблюдения, параметры которого мы хотим узнать, обозначим его [math]\displaystyle{ \Theta, i = 0 \ldots m }[/math].

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

Рисунок 1. Схема аппроксимативных связей начальных параметров.

При этом в виду особенностей B-сплайнов, о которых было упомянуто ранее мы получим, что они зависят от определённых значений параметра u, следовательно узнав при каком значении параметра U будет выполнено наблюдение [math]\displaystyle{ \Theta }[/math] и благодаря связи параметров мы получим искомые значения. Или же учитывая введённые обозначения получим связь параметров как:

\begin{equation} x_k^* \approx \sum_{i = 0}^{m}\frac{U^i(1-U)^{m-i}}{iB(i,m-i+1)}x_{ki}, \end{equation} где [math]\displaystyle{ k = \overline{1,n} }[/math], n - количество исходных параметров системы, [math]\displaystyle{ x_k^* }[/math] вычисленные параметры системы, при которых возможно наблюдение [math]\displaystyle{ \Theta }[/math]

Как писалось в пункте 1.1 для введения регуляризатора нам необходимо, чтобы выполнялся альтернанс Чебышёва, то есть чтобы многочлен [math]\displaystyle{ Q_m(u) = \sum_{i = 0}^{m}\frac{U^i(1-U)^{m-i}}{iB(i,m-i+1)}\vec{x} }[/math] при i = const, степени m был многочленом наилучшего равномерного приближения непрерывной функции [math]\displaystyle{ F(\vec{x_i}) }[/math], необходимо и достаточно существования на [math]\displaystyle{ \left[a,\ \ b\right] }[/math], где [math]\displaystyle{ a=\min{Q_m\left(u\right)} }[/math], [math]\displaystyle{ b=\max{Q_m\left(u\right)} }[/math] по крайней мере m + 2 точек [math]\displaystyle{ x_1\lt \ldots\lt x_{m+1} }[/math] таких, что \begin{equation}\label{iline-shevtsova-eq1} F(\vec{x_i})-Q_m(u)=\alpha(-1)^i||{F- Q_m (u)}||, \end{equation} где [math]\displaystyle{ i=1 \ldots m+1, \alpha= \pm1 }[/math] одновременно для всех i.

Или с учётом введённых обозначений и интегральной нормы гильбертова пространства

\begin{equation} \Theta-Q_m\left(u\right)=\alpha\left(-1\right)^i\sqrt{\int_{0}^{1}{\left[\Theta-\ Q_m\left(u\right)\right]\ du}}, \end{equation}

Тогда для решения нашей исходной задачи нам необходимо ввести сетку для параметра u, в данной реализации вводится равномерная сетка с шириной клетки в [math]\displaystyle{ \frac{1}{m} }[/math] и вычисляется

\begin{equation} U = \operatorname*{argmin}_u ((Q_m (u) - \Theta)^2 + \Theta - \int_{0}^{1}Q_m (u)\,du) \end{equation}

А находя значение [math]\displaystyle{ \int_{0}^{1}Q_m (u)\,du }[/math] получим окончательное условие в для решения обратной задачи, как

\begin{equation} U = \operatorname*{argmin}_u ((\sum_{i=0}^{m}B_{mi}(u)\theta_i(u) - \Theta)^2 + \Theta - \delta_m) \end{equation} где [math]\displaystyle{ \theta_i\left(u\right) }[/math] – i-е значение наблюдения, [math]\displaystyle{ B_{mi}(u) }[/math], [math]\displaystyle{ \delta_m }[/math]– расчётные коэффициенты, которые вычисляются согласно формулам

\begin{equation} B_{mi}(u) = \frac{U^i(1-U)^{m-i}}{iB(i,m-i+1)} , \end{equation}

\begin{equation} \delta_m = \sum_{i=0}^{m}\theta_i(u)\frac{m(1+i)^2}{i\Gamma(i)^2}\sum_{j=0}^{m-i}\frac{(-1)^j B(m+i, j+1)B(m,i)}{\Gamma(j+1)^2(j+i+1)}, \end{equation}

где [math]\displaystyle{ B(x,y) }[/math] - Бета функция, [math]\displaystyle{ \Gamma(x) }[/math] - Гамма функция

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

Вычислительное ядро для данного алгоритма по сути модно поделить на две части. Первая часть - это операции связанные с введением регулятора и частичным нахождением метрики для решения обратной задачи. То есть в данном алгоритме большую часть времени превалирует операция по расчёту [math]\displaystyle{ B_{mi}(u) }[/math], которую можно рассчитать как: \begin{equation} B_{mi}(u) = \frac{U^i(1-U)^{m-i}}{iB(i,m-i+1)}. \end{equation} При вычислении данного коэффициента, основное время работы занимают операции возведения в степень [math]\displaystyle{ U^i }[/math] и [math]\displaystyle{ (1-U)^{m-i} }[/math], а также вычисление Бета-функции.

Вычисление Бета-функции в области натуральных чисел, которым и являются [math]\displaystyle{ i }[/math] и [math]\displaystyle{ m-i+1 }[/math], можно свести к вычислению количества сочетаний, то есть к: \begin{equation} \frac{m!}{i!(m-i)!}. \end{equation} Реализация данного алгоритма состоит из последовательных умножений и делений. Второй частью вычислительного ядра является является расчёт параметра: \begin{equation} \delta_m = \sum_{i=0}^{m}\theta_i(u)\frac{m(1+i)^2}{i\Gamma(i)^2}\sum_{j=0}^{m-i}\frac{(-1)^j B(m+i, j+1)B(m,i)}{\Gamma(j+1)^2(j+i+1)}, \end{equation} Ядром которого является вычисление гамма-функций, в свою очередь для натуральных чисел мы можем заменить его вычислением факториала по свойствам Гамма-функций, то есть: \begin{equation} \Gamma(n-1) = n! \end{equation} Таким образом мы свели наш алгоритм к вычислению множества сочетай, факториалов и возведений в степень, что и будет вычислительным ядром нашего алгоритма.

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

1. Инициализация и чтение данных

• Чтение пар [math]\displaystyle{ (x, y) }[/math] из CSV файла
• Разделение на два вектора: [math]\displaystyle{ \vec x }[/math] и [math]\displaystyle{ \vec y }[/math]

2. Создание сетки параметра u

• Генерация равномерной сетки от 0.01 до 0.99 с шагом 0.001

3. Вычисление полиномов Бернштейна (функция bernstein_polynomial)

• Для каждого u вычисляется сумма:[math]\displaystyle{ \sum x_i × C_{m-1}^i × u^i × (1-u)^{m-1-i} }[/math]
• Используется эффективный алгоритм с накоплением

4. Вычисление метрики (функция calculate_metric)

• Вычисление основного члена: [math]\displaystyle{ (Bernstein(u) - f)^2 + f }[/math]
• Вычисление регуляризации через бета- и гамма-функции
• Использование кэша факториалов для ускорения

5. Поиск оптимума

• Поиск минимального значения метрики по сетке [math]\displaystyle{ u }[/math]
• Вывод оптимального значения [math]\displaystyle{ x }[/math]

6. Сохранение результатов

• Запись аппроксимации в output_aprox.dat
• Запись значений метрики в output_metric.dat

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

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

1. Инициализация и ввод данных

• Загрузка экспериментальных данных из CSV-файла
• Парсинг и преобразование строковых данных в числовые значения
• Разделение данных на векторы [math]\displaystyle{ \vec x }[/math] (входные параметры) и [math]\displaystyle{ \vec y }[/math] (наблюдения)
• Инициализация параметрической сетки [math]\displaystyle{ U = [0.01, 0.02, ..., 0.99] }[/math]

2. Вычисление комбинаторных коэффициентов

• Для каждого элемента сетки [math]\displaystyle{ u ∈ U }[/math] и каждого индекса [math]\displaystyle{ i ∈ [0, m-1] }[/math]:
  • Вычисление биномиального коэффициента [math]\displaystyle{ C(m, m-i) }[/math]
  • Вычисление степенных членов [math]\displaystyle{ u^i }[/math] и [math]\displaystyle{ (1-u)^{(m-i)} }[/math]
  • Формирование взвешенной суммы элементов вектора [math]\displaystyle{ \vec x }[/math]

3. Решение прямой задачи аппроксимации

• Для каждого параметра [math]\displaystyle{ u ∈ U }[/math]:
  • Инициализация накопительной суммы [math]\displaystyle{ S = 0 }[/math]
  • Для каждого [math]\displaystyle{ i ∈ [0, m-1] }[/math]:
    • Вычисление комбинаторного веса: [math]\displaystyle{ w = C(m, m-i) × u^i × (1-u)^{(m-i)} }[/math]
    • Добавление взвешенного элемента: [math]\displaystyle{ S = S + w × x[i] }[/math]
  • Сохранение результата [math]\displaystyle{ F(x, u) = S }[/math]

4. Вычисление регуляризационного слагаемого

• Инициализация суммы регуляризации [math]\displaystyle{ R = 0 }[/math]
• Для каждого [math]\displaystyle{ i ∈ [1, m-1] }[/math]:
  • Инициализация внутренней суммы [math]\displaystyle{ S_{reg} = 0 }[/math]
  • Для каждого [math]\displaystyle{ j ∈ [0, m-i-1] }[/math]:
    • Вычисление знака: [math]\displaystyle{ sign = (-1)^j }[/math]
    • Вычисление бета-функций: [math]\displaystyle{ B1 = β(m+i, j+1), B2 = β(m, i) }[/math]
    • Вычисление гамма-функций: [math]\displaystyle{ Γ1 = Γ(j+1), Γ2 = Γ(j+1) }[/math]
    • Формирование слагаемого: [math]\displaystyle{ term = \frac{sign × B1 × B2}{(Γ1 × Γ2 × (j+i+1))} }[/math]
    • [math]\displaystyle{ S_{reg} = S_{reg} + term }[/math]
  • Вычисление весового множителя: [math]\displaystyle{ weight = \frac{m × (1+i)^2 }{(i × Γ(i)^2)} }[/math]
  • [math]\displaystyle{ R = R + S_{reg} × y[i] × weight }[/math]

5. Формирование целевого функционала с регуляризацией

• Для каждого параметра [math]\displaystyle{ u ∈ U }[/math]:
  • Вычисление невязки: [math]\displaystyle{ residual = F(y, u) - f }[/math]
  • Вычисление квадратичного члена: [math]\displaystyle{ fidelity_{term} = residual^2 }[/math]
  • Формирование полного функционала: [math]\displaystyle{ J(u) = fidelity_{term} + f - R }[/math]

6. Минимизация функционала и определение оптимального решения

• Инициализация минимального значения [math]\displaystyle{ J_{min} = \infty }[/math] и индекса [math]\displaystyle{ min_{index} = 0 }[/math]
• Для каждого [math]\displaystyle{ u ∈ U }[/math]:
  • Если [math]\displaystyle{ J(u) \lt J_{min} }[/math]:
    • Обновление [math]\displaystyle{ J_{min} = J(u) }[/math]
    • Сохранение [math]\displaystyle{ min_{index} = index(u) }[/math]
• Определение оптимального решения: [math]\displaystyle{ x_{opt} = F(x, U[min_{index}]) }[/math]

7. Вывод результатов

• Отображение оптимального значения [math]\displaystyle{ x_{opt} }[/math]
• Завершение работы алгоритма

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

Как писалось ранее основными функциями, которые входят в вычислительное ядро алгоритма - то вычисление полинома Бернштейна [math]\displaystyle{ B_{mi}(u) }[/math] и вычисление метрики, в которую основную сложность вносит параметр [math]\displaystyle{ \delta_m }[/math]. Рассмотрим их сложность вычисления:

1. Вычисление [math]\displaystyle{ B_{mi}(u) }[/math]:

• Инициализация, содержит [math]\displaystyle{ n }[/math] операций. Сложность - [math]\displaystyle{ O(n) }[/math]
• Возведение в степень, содержит [math]\displaystyle{ n×O(\log m) }[/math] операций. Сложность - [math]\displaystyle{ O(n \log m) }[/math]
• Основной цикл, содержит [math]\displaystyle{ n × m }[/math] операций. Сложность - [math]\displaystyle{ O(n×m) }[/math]

Итого получили сложность [math]\displaystyle{ O(n×m) }[/math]

2. Вычисление метрики:

• Кэш факториалов. Сложность - [math]\displaystyle{ O(m) }[/math]
• Регуляризация. Сложность - [math]\displaystyle{ O(m^2) }[/math]
• Конечная метрика. Сложность - [math]\displaystyle{ O(n) }[/math]

Итого получили сложность [math]\displaystyle{ O(m^2 + n) }[/math]

Таким образом финальная сложность алгоритма будет составлять:[math]\displaystyle{ O(n×m + m^2) }[/math]

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

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

Рисунок 2. Информационный граф алгоритма

Высокоуровневые узлы и их семантика.
Граф можно декомпозировать на семь логических узлов:

1. main — корневой узел, инициирующий выполнение всего алгоритма. Он задаёт параметры и запускает независимые процессы чтения данных и построения вычислительной сетки.

2. readData — узел ввода. Выполняет синтаксический разбор CSV-файла data.csv, используя регулярные выражения для токенизации строк, и формирует два вектора исходных данных: x (аргументы) и y (значения).

3. buildGrid — узел генерации области исследования. Создаёт равномерную одномерную сетку U в диапазоне от 0.01 до 0.99 с фиксированным шагом 0.001, определяя множество точек, в которых будет оцениваться интерполирующая функция.

4. func — вычислительное ядро сглаживания. Для каждой точки u из сетки U вычисляет взвешенную сумму элементов вектора x. Веса рассчитываются на основе заданных параметров (C_func, pow_m, gcd), что реализует непараметрическую оценку функции.

5. metric — вычислительное ядро анализа. Для каждой точки u вычисляет значение целевой метрики, используя исходные данные y, результаты работы узла func для текущего u, а также ряд специальных функций (beta_function, gamma_int, factorial). Данный узел содержит наиболее сложные и затратные операции, включая вложенные суммирования и работу с числами произвольной точности (GMP).

6. minSearch — узел редукции. Выполняет последовательный линейный поиск минимального значения в массиве metric_arr, полученном от узла metric, и определяет индекс min_idx, соответствующий оптимальной точке сетки.

7. output — терминальный узел. Формирует итоговый результат работы алгоритма, выводя значение U[min_idx] — оценку аргумента, минимизирующего целевую метрику.

Потоки данных и зависимости. Рёбра графа отображают передачу данных и управляющие зависимости:

• Узел main активирует параллельное выполнение узлов readData и buildGrid.
• Узлы func и metric являются потребителями исходных данных: они начинают работу после завершения readData и получения векторов x и y.
• Оба вычислительных ядра (func и metric) также требуют задания области исследования, поэтому зависят от узла buildGrid, поставляющего сетку U.
• Между ядрами существует строгая зависимость по данным: для вычисления метрики в точке u необходимо предварительно рассчитать значение сглаженной функции в этой же точке. Таким образом, узел metric для конкретного u ожидает результат работы узла func для того же u.
• После полного формирования массива метрик управление передаётся узлу minSearch, который, завершив поиск минимума, инициирует финальный вывод в узле output.

Ярусно-параллельная форма (ЯПФ). Вычислительный процесс естественным образом распределяется по пяти последовательным ярусам:

Ярус 0: Содержит единственный узел main, запускающий вычисления.

Ярус 1: Состоит из двух независимых узлов readData и buildGrid, которые могут выполняться параллельно. Это максимальный по ширине ярус.

Ярус 2: Включает вычислительные узлы func и metric. Между ними существует тонкая зависимость: начало вычисления метрики для точки u возможно только после получения сглаженного значения для этой точки от func. Однако, если рассматривать обработку всей сетки U, эти узлы могут работать в конвейерном режиме: пока func вычисляет значение для точки u_{k+1}, узел metric может обрабатывать предыдущую точку u_k.

Ярус 3: Представлен узлом minSearch, который начинает работу после полного завершения расчётов на предыдущем ярусе и выполняет редукционную операцию.

Ярус 4: Финальный ярус с узлом output, выводящим результат.

Таким образом, высота ЯПФ равна 5, что характеризует длину критического пути. Ширина графа (степень параллелизма) достигает максимума на первом вычислительном ярусе (равна 2) и теоретически может быть расширена до размера сетки U на втором ярусе при распараллеливании циклов по точкам. Информационный граф наглядно демонстрирует, что основная вычислительная нагрузка сосредоточена на втором ярусе, а потенциальным узким местом является последовательный поиск минимума и, в меньшей степени, этапы ввода-вывода.

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

Анализ информационного графа и структуры вычислений демонстрирует наличие значительного потенциала для параллельного исполнения алгоритма. Ресурс параллелизма проявляется на нескольких уровнях, что позволяет адаптировать реализацию к различным вычислительным архитектурам — от многоядерных процессоров до распределённых кластерных систем.

1. Параллелизм на уровне макроопераций (задач). На верхнем уровне алгоритм содержит независимые ветви исполнения, которые могут выполняться одновременно:

• Процесс чтения входных данных из файла (readData) полностью независим от построения вычислительной сетки (buildGrid). Эта независимость позволяет выполнять эти операции параллельно, что особенно актуально при работе с большими файлами данных или при необходимости динамического построения сложных сеток.

• Оба указанных процесса не требуют синхронизации и могут быть запущены одновременно сразу после инициализации параметров в узле main.

2. Параллелизм на уровне данных (data parallelism). Наиболее существенный ресурс параллелизма заключён в обработке точек вычислительной сетки [math]\displaystyle{ U }[/math]. Основные вычислительные ядра алгоритма обладают следующими характеристиками:

• Функция func: Вычисление сглаженного значения в точке u_k не зависит от вычислений в точках [math]\displaystyle{ u_m }[/math] при [math]\displaystyle{ k ≠ m }[/math]. Таким образом, обработка всех точек сетки может быть распределена между вычислительными элементами без необходимости обмена промежуточными данными.

• Функция metric: Аналогично, расчёт значения метрики в точке[math]\displaystyle{ u_k }[/math]зависит только от сглаженного значения в этой же точке и исходных данных, но не от значений метрики в других точках. Это позволяет независимо и параллельно вычислять метрику для разных точек сетки.

Теоретическая оценка: Максимальная степень параллелизма на этом уровне ограничена только размером сетки [math]\displaystyle{ N = |U| }[/math] и может достигать [math]\displaystyle{ O(N) }[/math], что для типичных параметров (сетка с шагом 0.001 на интервале [math]\displaystyle{ [0.01, 0.99) }[/math]) составляет порядка 980 параллельно выполняемых задач.

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

Взвешенное суммирование в func: Вычисление суммы [math]\displaystyle{ S(u) = \sum w_i(u) * x_i }[/math] по всем входным точкам [math]\displaystyle{ i = 1..n }[/math] представляет собой операцию редукции, которая может быть эффективно распараллелена с использованием техник параллельного суммирования (например, методом сдваивания).

Вложенные циклы в metric: Расчёт метрики содержит двойные и тройные вложенные циклы по индексам входных данных. Эти циклы могут быть распараллелены, однако требуется аккуратное обращение с операциями над числами произвольной точности (GMP), которые могут создавать contention при работе с общими ресурсами памяти.

4. Конвейерный параллелизм (pipelining). Между основными вычислительными стадиями существует возможность организации конвейерной обработки:

После вычисления сглаженного значения для первой точки сетки [math]\displaystyle{ u_1 }[/math] в модуле func, модуль metric может начать обработку этой точки, в то время как func продолжает вычисления для точки [math]\displaystyle{ u_2 }[/math].

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

5. Параллелизм на уровне инструкций (ILP). Благодаря регулярной структуре вычислений и отсутствию сложных условных переходов в основных циклах, алгоритм обладает хорошим потенциалом для использования возможностей современных процессоров:

Векторизация (SIMD): операции с плавающей точкой и целочисленной арифметикой в циклах могут быть векторизованы с использованием инструкций AVX/SSE (для частей кода, не задействующих GMP).

Суперскалярное исполнение: независимые инструкции внутри итераций циклов могут выполняться параллельно на нескольких исполнительных устройствах процессора.

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

• Чтение данных из файла — операция ввода-вывода, которая может стать узким местом, хотя и может быть перекрыта с другими вычислениями.

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

• Инициализация и финализация GMP-объектов — некоторые операции с числами произвольной точности требуют последовательного доступа к управляющим структурам.

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

Входные данные
Алгоритм принимает на вход набор экспериментальных данных, представленный в виде пар значений [math]\displaystyle{ (x_i, y_i) }[/math], где [math]\displaystyle{ i = 1, 2, ..., m }[/math]. Данные считываются из текстового файла в формате CSV (Comma-Separated Values), где каждая строка содержит два числа, разделенных запятой или пробелом: координату [math]\displaystyle{ x_i }[/math] и соответствующее значение [math]\displaystyle{ y_i }[/math]. Формат файла допускает использование различных разделителей (запятые, пробелы, переводы строк) благодаря применению регулярных выражений при парсинге.

Дополнительным входным параметром алгоритма является вещественное число [math]\displaystyle{ f }[/math], которое используется в расчете метрики качества. В представленной реализации значение [math]\displaystyle{ f }[/math] жестко задано как 20.0, но может быть легко модифицировано для конфигурации через аргументы командной строки или конфигурационный файл.

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

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

• Первое число — значение аппроксимации для координаты x
• Второе число — значение аппроксимации для координаты y

Файл output_metric.dat — содержит значения метрики качества для различных значений параметра u. Каждая строка содержит:

Значение параметра [math]\displaystyle{ u }[/math] в диапазоне [math]\displaystyle{ [0.01, 0.99] }[/math] с шагом 0.001

Соответствующее значение метрики качества с точностью 12 знаков после десятичной точки

Стандартный вывод — в консоль выводится оптимальное значение аппроксимации, соответствующее минимальному значению метрики качества. Также выводится информация о времени выполнения программы с указанием количества использованных MPI-процессов и OpenMP-нитей.

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

1. Детерминированность Алгоритм является детерминированным — при одинаковых входных данных и одинаковых параметрах параллельного выполнения (количестве MPI-процессов и OpenMP-нитей) результаты вычислений будут идентичны. Это обеспечивается:

• Отсутствием случайных элементов в вычислениях
• Детерминированными формулами для расчета значений функций
• Фиксированными диапазонами параметров (U от 0.01 до 0.99 с шагом 0.001)

2. Устойчивость Алгоритм обладает высокой численной устойчивостью благодаря:

• Использованию библиотеки GMP (GNU Multiple Precision) для вычислений с произвольной точностью
• Реализации собственных функций для операций с высокой точностью (fast_pow, gcd)
• Кэшированию значений факториалов для ускорения вычислений гамма- и бета-функций
• Обработке особых случаев (например, деления на ноль в вычислении отношения)

3. Точность

• Входные данные: представлены числами с плавающей запятой двойной точности (double)
• Промежуточные вычисления: выполняются с использованием библиотеки GMP, обеспечивающей произвольную точность
• Выходные данные: выводятся с точностью 12 знаков после десятичной точки
• Контроль точности: обеспечивается за счет использования точных рациональных вычислений для суммирования

4. Вычислительная мощность Алгоритм характеризуется высокой вычислительной сложностью:

• Основные вычисления сосредоточены в функциях func и metric
• Количество операций растет как [math]\displaystyle{ O(m × n) }[/math], где [math]\displaystyle{ m }[/math] — количество экспериментальных точек, [math]\displaystyle{ n }[/math] — количество значений параметра [math]\displaystyle{ U }[/math]

Используются вычисления комбинаторных коэффициентов через гамма- и бета-функции

5. Вычислительная сложность

• Временная сложность: [math]\displaystyle{ O(m^2 × n) }[/math] в наихудшем случае, где [math]\displaystyle{ m }[/math] — размер входных данных, [math]\displaystyle{ n }[/math] — количество точек параметра [math]\displaystyle{ U }[/math]
• Пространственная сложность: [math]\displaystyle{ O(m + n) }[/math] для хранения входных данных и промежуточных результатов
• Сложность параллельной части: эффективно распараллеливается как по данным (разные значения [math]\displaystyle{ U }[/math]), так и по операциям (внутренние циклы)

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

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

В базовом варианте алгоритм реализуется через последовательное вычисление значений метрики качества для каждого параметра [math]\displaystyle{ u }[/math]. Основные вычислительные ядра алгоритма представлены следующими функциями:

Функция вычисления базовой суммы (func)

vector<MPFWrapper> func(const vector<double> &x, const vector<double> &U) {
    size_t m = x.size();
    vector<MPFWrapper> sum_arr(U.size());

    for (int idx = 0; idx < static_cast<int>(U.size()); ++idx) {
        mpf_t sum, u_val, one_minus_u, term, ratio, temp1, temp2, scale;
        // Инициализация переменных GMP
        
        mpf_set_d(u_val, U[idx]);
        mpf_set_d(one_minus_u, 1.0);
        mpf_sub(one_minus_u, one_minus_u, u_val);
        
        // Быстрое возведение в степень (1-u)^m
        fast_pow(term, one_minus_u, static_cast<int>(m));
        
        // Вычисление отношения u/(1-u) с обработкой особого случая
        if (mpf_cmp_ui(one_minus_u, 0) == 0) {
            mpf_set_ui(ratio, 0);
        } else {
            mpf_div(ratio, u_val, one_minus_u);
        }

        for (int i = 0; i < static_cast<int>(m); i++) {
            // Добавление слагаемого: x[i] * term
            mpf_set_d(temp1, x[i]);
            mpf_mul(temp2, term, temp1);
            mpf_add(sum, sum, temp2);
            
            // Рекуррентное обновление term для следующего i:
            // term = term * (m-i)/(i+1) * u/(1-u)
            if (i < static_cast<int>(m) - 1 && mpf_cmp_ui(one_minus_u, 0) != 0) {
                mpf_set_d(temp1, static_cast<double>(m - i));
                mpf_set_d(temp2, static_cast<double>(i + 1));
                mpf_div(scale, temp1, temp2);
                mpf_mul(temp1, scale, ratio);
                mpf_mul(term, term, temp1);
            }
        }
        
        mpf_set(sum_arr[idx].value, sum);
        // Освобождение переменных GMP
    }
    return sum_arr;
}

Функция вычисления метрики (metric)

vector<MPFWrapper> metric(const vector<double> &x,
                          const vector<double> &U,
                          double f) {
    vector<MPFWrapper> sum_1 = func(x, U);
    int m = static_cast<int>(x.size());

    // Построение кэша факториалов
    const int gamma_limit = max(3 * m, 1);
    vector<MPFWrapper> factorial_cache = build_factorial_cache(gamma_limit);

    MPFWrapper sum_2_wrapper;
    mpf_set_ui(sum_2_wrapper.value, 0);

    // Вычисление регуляризирующего слагаемого
    for (int i = 1; i < m; i++) {
        MPFWrapper sum_2_regulation_wrapper;
        mpf_set_ui(sum_2_regulation_wrapper.value, 0);
        
        int sign = 1;
        mpf_t beta_m_i, gamma_i, numerator, gamma_j, denominator, temp1, temp2, temp3, coeff;
        // Инициализация переменных GMP
        
        // Вычисление B(m, i) через гамма-функции
        beta_cached(beta_m_i, m, i, factorial_cache);
        gamma_cached(gamma_i, factorial_cache, i);
        
        for (int j = 0; j < m - i; j++) {
            // Вычисление B(m+i, j+1)
            beta_cached(temp1, m + i, j + 1, factorial_cache);
            mpf_mul(numerator, temp1, beta_m_i);
            
            // Вычисление (j+1)! = gamma(j+2)
            gamma_cached(gamma_j, factorial_cache, j + 1);
            
            // Знаменатель: (j+1)!^2 * (j+i+1)
            mpf_mul(temp1, gamma_j, gamma_j);
            mpf_set_d(temp2, static_cast<double>(j + i + 1));
            mpf_mul(denominator, temp1, temp2);
            
            // Слагаемое с чередующимся знаком
            mpf_div(temp1, numerator, denominator);
            if (sign < 0) {
                mpf_neg(temp1, temp1);
            }
            mpf_add(sum_2_regulation_wrapper.value, sum_2_regulation_wrapper.value, temp1);
            sign = -sign;
        }
        
        // Коэффициент C_i = m * (1+i)^2 / (i * (i-1)!^2)
        mpf_set_d(temp1, static_cast<double>(m));
        mpf_set_d(temp2, static_cast<double>(1 + i));
        mpf_mul(temp3, temp2, temp2);
        mpf_mul(temp1, temp1, temp3);
        mpf_set_d(temp2, static_cast<double>(i));
        mpf_mul(temp3, temp2, gamma_i);
        mpf_mul(temp3, temp3, gamma_i);
        mpf_div(coeff, temp1, temp3);
        
        // Умножение на x[i] и коэффициент
        mpf_set_d(temp1, x[i]);
        mpf_mul(temp2, sum_2_regulation_wrapper.value, temp1);
        mpf_mul(temp1, temp2, coeff);
        mpf_add(sum_2_wrapper.value, sum_2_wrapper.value, temp1);
    }

    // Финальное вычисление метрики
    vector<MPFWrapper> out(U.size());
    for (int u_idx = 0; u_idx < static_cast<int>(U.size()); u_idx++) {
        mpf_t diff, temp1, temp2;
        
        // M(u) = (S(u) - f)^2 + f - R
        mpf_set_d(temp1, f);
        mpf_sub(diff, sum_1[u_idx].value, temp1);
        mpf_mul(temp1, diff, diff);
        mpf_set_d(temp2, f);
        mpf_add(temp1, temp1, temp2);
        mpf_sub(out[u_idx].value, temp1, sum_2_wrapper.value);
    }

    return out;
}

Помимо основных оптимизированных вычислительных алгоритмов, особое внимание уделялось "подсобным" операциями таким как:

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

void fast_pow(mpf_t result, mpf_t base, int exp) {
    mpf_set_ui(result, 1);
    if (exp == 0) return;
    
    mpf_t multiplier, temp;
    mpf_init_set(multiplier, base);
    mpf_init(temp);
    
    int power = exp;
    while (power > 0) {
        if (power & 1) {
            mpf_mul(temp, result, multiplier);
            mpf_set(result, temp);
        }
        power >>= 1;
        if (power) {
            mpf_mul(temp, multiplier, multiplier);
            mpf_set(multiplier, temp);
        }
    }
    
    mpf_clear(multiplier);
    mpf_clear(temp);
}

Этот алгоритм имеет сложность [math]\displaystyle{ O(\log n) }[/math] вместо [math]\displaystyle{ O(n) }[/math] при наивном умножении. Важно отметить, что функция работает с типом mpf_t из библиотеки GMP, что обеспечивает высокую точность вычислений даже для больших степеней.

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

1. Кэширование факториалов:

vector<MPFWrapper> build_factorial_cache(int limit) {
    vector<MPFWrapper> cache(limit + 1);
    mpf_set_ui(cache[0].value, 1);
    
    mpf_t temp;
    mpf_init(temp);
    for (int i = 1; i <= limit; ++i) {
        mpf_set_d(temp, static_cast<double>(i));
        mpf_mul(cache[i].value, cache[i - 1].value, temp);
    }
    mpf_clear(temp);
    return cache;
}

2. Вычисление гамма-функции для целых аргументов через факториал:

void gamma_cached(mpf_t result, const vector<MPFWrapper> &cache, int n) {
    // Γ(n) = (n-1)!
    mpf_set(result, cache[n - 1].value);
}

3. Вычисление бета-функции через гамма-функции:

void beta_cached(mpf_t result, int a, int b, const vector<MPFWrapper> &cache) {
    // B(a,b) = Γ(a)Γ(b)/Γ(a+b)
    mpf_t temp1, temp2, temp3;
    mpf_init(temp1);
    mpf_init(temp2);
    mpf_init(temp3);
    
    gamma_cached(temp1, cache, a);
    gamma_cached(temp2, cache, b);
    gamma_cached(temp3, cache, a + b);
    
    mpf_mul(temp1, temp1, temp2);
    mpf_div(result, temp1, temp3);
    
    mpf_clear(temp1);
    mpf_clear(temp2);
    mpf_clear(temp3);
}

Использование библиотеки GMP для высокой точности
Особенностью данной реализации является использование библиотеки GNU Multiple Precision (GMP) для всех критичных вычислений. Это обеспечивает:

1.Произвольную точность: В отличие от стандартных типов double (53 бита мантиссы) или long double (64 бита), GMP позволяет работать с числами произвольной точности, ограниченной только доступной памятью.

2.Контроль над округлением: При умножении двух чисел с плавающей запятой стандартных типов результат получается с удвоенной длиной мантиссы и порядка, и только при присваивании происходит округление до точности типа переменной. В GMP это поведение контролируется явно через установку точности чисел.

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

Для интеграции типов GMP с контейнерами STL реализованы специальные обертки:

struct MPFWrapper {
    mpf_t value;
    MPFWrapper() { mpf_init(value); }
    MPFWrapper(const MPFWrapper& other) {
        mpf_init(value);
        mpf_set(value, other.value);
    }
    ~MPFWrapper() { mpf_clear(value); }
    // ... операторы присваивания и приведения типов
};

Эти обертки обеспечивают корректное управление памятью GMP (инициализация, копирование, освобождение) в контексте C++ RAII (Resource Acquisition Is Initialization).

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

Как видно из структуры алгоритма, вариантов распараллеливания существует несколько. Алгоритм демонстрирует высокую степень параллелизма, поскольку вычисления для разных значений параметра [math]\displaystyle{ U }[/math] являются полностью независимыми. Это позволяет эффективно распределять работу между вычислительными ресурсами без необходимости сложной синхронизации.

В представленной реализации используются два основных подхода к параллелизации:

1. Межпроцессный параллелизм (MPI)
Основной массив вычислений разбивается на блоки по значениям параметра [math]\displaystyle{ U }[/math], которые распределяются между MPI-процессами. Каждый процесс работает со своей порцией данных, что минимизирует необходимость в обменах информацией во время вычислений.

int total_U = static_cast<int>(U.size());
int base = total_U / world_size;
int remainder = total_U % world_size;
int local_count = base + (world_rank < remainder ? 1 : 0);

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

2. Внутрипроцессный параллелизм (OpenMP)
Внутри каждого MPI-процесса используется параллелизм на уровне потоков с помощью технологии OpenMP. Особенностью реализации является применение динамического планирования:

#pragma omp parallel for schedule(dynamic)
for (int idx = 0; idx < static_cast<int>(U.size()); ++idx) {
    // Вычисления для каждого значения U
}

Динамическое планирование позволяет адаптироваться к возможной неоднородности вычислений для разных значений параметра [math]\displaystyle{ U }[/math], что улучшает балансировку нагрузки между потоками внутри процесса.

Основная сложность при параллелизации данного алгоритма связана с использованием библиотеки точной арифметики GMP (GNU Multiple Precision Arithmetic Library). Библиотека GMP не является потокобезопасной, что требует специальных мер при использовании в многопоточной среде.

В реализации эта проблема решается следующим образом:

1. Локальные переменные в каждом потоке - каждая параллельная область создает собственные экземпляры переменных GMP, что исключает конфликты при обращении к разделяемым данным.

2. Обертки с автоматическим управлением памятью - используются структуры MPFWrapper и MPZWrapper, которые обеспечивают безопасное копирование и освобождение ресурсов:

struct MPFWrapper {
    mpf_t value;
    MPFWrapper() { mpf_init(value); }
    ~MPFWrapper() { mpf_clear(value); }
    // ... методы копирования и присваивания ...
};

3. Минимизация критических секций - синхронизация между потоками сведена к минимуму и используется только там, где это действительно необходимо (например, при накоплении суммы).

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

Для рассмотрения работы алгоритма для разных случаев вычислений лучше всего подходит решение задачи пролонгации значений. Так алгоритм позволяет пролонгировать значение например 515 на другие подробные сетки и при этом решение обратной задачи никак не изменится, но при этом сильно растёт вычислительная нагрузка и время расчёта. Рассмотрим результат прогонов одной и той же задачи сначала с пролонгацией или точностью, как указано на графике, сначала на [math]\displaystyle{ 10^3 }[/math] точек, потом на [math]\displaystyle{ 10^4 }[/math], [math]\displaystyle{ 10^5 }[/math] и на [math]\displaystyle{ 10^6 }[/math] отразим зависимость времени вычисления от количества процессоров - Рис. 3

Рисунок 3. Зависимость времени выполнения задачи о пролонгации 515 точек на различные точности от числа процессоров в логарифмическом масштабе

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

Рисунок 4. Зависимость ускорения выполнения задачи о пролонгации 515 точек на различные точности от числа процессоров

На основе этих данных мы можем сделать вывод, что при увеличении сложности решения задачи всё больше будет эффект или ускорение от её паралелизации и тем ближе мы к линейному ускорению. Также если обратимся к не сложным решениям с малым числом точек - при точности в [math]\displaystyle{ 10^3 }[/math] и [math]\displaystyle{ 10^4 }[/math] мы можем даже заметить увеличение времени расчёта программы при большом числе процессоров - это связанно с тем, что время вычислений на одном процессоре меньше, чем накладные расходы на выделение памяти, распределения ресурсов и т.п.

Таким образом, модно сделать вывод о том, что данный алгоритм в своей параллельной реализации имеет большой вес при работе с большим числом данных или их размерностью, а также нахождения решения с высокой точностью, то есть при погрешности выходных значений ниже [math]\displaystyle{ 10^{-4} }[/math]