Участник:Морозов Егор Павлвочи/Метод решения недетерминированных обратных задач идентификации входных параметров с 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). Граф имеет ярко выраженную ярусно-параллельную структуру, что указывает на наличие значительного ресурса параллелизма, в основном на уровне обработки независимых точек заданной сетки.

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

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).