Алговики - Вклад [ru]

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

2025-12-07T22:42:04Z

Морозов Егор Павлвочи: Добавлены результаты прогонов и параллельное описание алгоритма

{{algorithm
| name = Метод решения недетерминированных обратных задач идентификации входных параметров с BS регуляризацией
| serial_complexity = <math>O(nm+m^2)</math>
| pf_height = <math>5</math>
| pf_width = <math>O(n)</math>
| input_data = <math>2m+n+1</math>
| output_data = <math>4n+1</math>
}}

Основные авторы описания: [[Участник:Морозов Егор Павлвочи|Е.П.Морозов]]

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

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

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

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

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

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

В математической литературе длительное время существовала точка зрения, согласно которой всякая математическая задача должна удовлетворять этим требованиям <ref>Курант Р., Вентцель Т. Д. Уравнения с частными производными. – Москва.– : Мир, 1964. – Т. 2. – С. 832.</ref>. Задачи, не удовлетворяющие перечисленным требованиям, называются некорректно поставленными. Следует отметить, что определение некорректно поставленных задач относится к данной паре метрических пространств (F, U), так как в других метриках та же задача может быть корректно поставленной. В рассматриваемом классе задач нарушается устойчивость и однозначность решения. Согласно А.Н. Тихонову, для таких задач существуют квазирешения. Однако чтобы перейти к решению, а поняв его структуру, — к вычислению обратной задачи, необходимо восстановить приблизительную зависимость между входными и выходными параметрами.

Для этого рассмотрим B-сплайны и NURBS. Их полное геометрическое описание и вывод основных свойств были представлены в «NURBS book» <ref>Piegl L., Tiller W. The NURBS book. – Springer Science & Business Media, 2012.</ref>. Помимо их основного применения — создания и сглаживания поверхностей и сложных форм на основе реперных точек, — мы можем применить их для аппроксимации исходных данных и получения системы алгебраических уравнений.

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

Данные системы функциональных уравнений и связанные с ними обратные задачи можно решать различными способами. Пример решения подобных задач для уравнений математической физики представлен в книге <ref>Самарский А. А., Вабищевич П. Н. Численные методы решения обратных задач математической физики. – 2004.</ref>, где они называются обратными задачами идентификации правой части. Однако представленный алгоритм ориентирован не на применение конечно-разностных схем, не на регуляризацию по Тихонову или иные традиционные методы, а на использование нового регуляризатора, основанного на B-сплайнах. При аппроксимации наших данных мы вводили B-сплайны как многочлены порядка m, обозначим их как <math>Q_m(u)</math>, где u - наш новый параметр. Для более точного решения нам необходимо, чтобы данный аппроксимативный многочлен удовлетворял условиям альтернанса Чебышева. Тем самым <math>Q_m(u)</math> будет полиномом наилучшего приближения. Тем самым мы можем за основу взять условия альтернанса
\begin{equation}
F(\vec{x_i})-Q_m(u)=\alpha(-1)^i||{F- Q_m (u)}||,
\end{equation}
где <math>i=1\ldots m+1</math>, <math>\alpha=\pm1</math> одновременно для всех i.
Для обобщения решения будем рассматривать Гильбертовы пространства и в качестве нормы введём интегральную гильбертову норму. А в качестве <math>F(\vec{x_i})</math> рассмотрим наше известное значение для наблюдения, идентификацию начальных параметров для которого мы хотим получить, обозначим его как <math>\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>\vec{\theta} = (\theta_1, \theta_2, \ldots, \theta_m)</math>, где m - количество наблюдений, которые были получены при начальных параметрах <math>\vec{x_i} = (x_{1i}, x_{2i}, \ldots, x_{ni})</math>, где n - количество исходных параметров для рассматриваемой системы, а также дано целевое значение наблюдения, параметры которого мы хотим узнать, обозначим его <math>\Theta, i = 0 \ldots m</math>.

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

[[file:Bspline_regulation_1.png|thumb|center|500px|Рисунок 1. Схема аппроксимативных связей начальных параметров.]]

При этом в виду особенностей B-сплайнов, о которых было упомянуто ранее мы получим, что они зависят от определённых значений параметра u, следовательно узнав при каком значении параметра U будет выполнено наблюдение <math>\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>k = \overline{1,n}</math>, n - количество исходных параметров системы, <math>x_k^*</math> вычисленные параметры системы, при которых возможно наблюдение <math>\Theta</math>

Как писалось в пункте 1.1 для введения регуляризатора нам необходимо, чтобы выполнялся альтернанс Чебышёва, то есть
чтобы многочлен <math>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> F(\vec{x_i}) </math>, необходимо и достаточно существования на <math>\left[a,\ \ b\right]</math>, где <math>a=\min{Q_m\left(u\right)}</math>, <math>b=\max{Q_m\left(u\right)}</math> по крайней мере m + 2 точек <math>x_1<\ldots<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> 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>\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>\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>\theta_i\left(u\right)</math> – i-е значение наблюдения, <math>B_{mi}(u)</math>, <math>\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>B(x,y)</math> - Бета функция, <math>\Gamma(x)</math> - Гамма функция

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

Вычисление Бета-функции в области натуральных чисел, которым и являются <math>i</math> и <math>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>(x, y)</math> из CSV файла
* Разделение на два вектора: <math>\vec x</math> и <math> \vec y</math>
2. Создание сетки параметра u
* Генерация равномерной сетки от 0.01 до 0.99 с шагом 0.001
3. Вычисление полиномов Бернштейна (функция bernstein_polynomial)
* Для каждого u вычисляется сумма:<math> \sum x_i × C_{m-1}^i × u^i × (1-u)^{m-1-i}</math>
* Используется эффективный алгоритм с накоплением
4. Вычисление метрики (функция calculate_metric)
* Вычисление основного члена: <math>(Bernstein(u) - f)^2 + f</math>
* Вычисление регуляризации через бета- и гамма-функции
* Использование кэша факториалов для ускорения
5. Поиск оптимума
* Поиск минимального значения метрики по сетке <math>u</math>
* Вывод оптимального значения <math>x</math>
6. Сохранение результатов
* Запись аппроксимации в output_aprox.dat
* Запись значений метрики в output_metric.dat

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

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

2. Вычисление комбинаторных коэффициентов
* Для каждого элемента сетки <math>u ∈ U</math> и каждого индекса <math>i ∈ [0, m-1]</math>:
** Вычисление биномиального коэффициента <math>C(m, m-i)</math>
** Вычисление степенных членов <math>u^i</math> и <math>(1-u)^{(m-i)}</math>
** Формирование взвешенной суммы элементов вектора <math>\vec x</math>
3. Решение прямой задачи аппроксимации
* Для каждого параметра <math>u ∈ U</math>:
** Инициализация накопительной суммы <math>S = 0</math>
** Для каждого <math>i ∈ [0, m-1]</math>:
*** Вычисление комбинаторного веса: <math>w = C(m, m-i) × u^i × (1-u)^{(m-i)}</math>
***Добавление взвешенного элемента: <math>S = S + w × x[i]</math>
**Сохранение результата <math>F(x, u) = S</math>
4. Вычисление регуляризационного слагаемого
* Инициализация суммы регуляризации <math>R = 0</math>
* Для каждого <math>i ∈ [1, m-1]</math>:
** Инициализация внутренней суммы <math>S_{reg} = 0</math>
** Для каждого <math>j ∈ [0, m-i-1]</math>:
*** Вычисление знака: <math>sign = (-1)^j</math>
*** Вычисление бета-функций: <math>B1 = β(m+i, j+1), B2 = β(m, i)</math>
*** Вычисление гамма-функций: <math>Γ1 = Γ(j+1), Γ2 = Γ(j+1)</math>
*** Формирование слагаемого: <math>term = \frac{sign × B1 × B2}{(Γ1 × Γ2 × (j+i+1))}</math>
*** <math>S_{reg} = S_{reg} + term</math>
** Вычисление весового множителя: <math>weight = \frac{m × (1+i)^2 }{(i × Γ(i)^2)}</math>
** <math>R = R + S_{reg} × y[i] × weight</math>
5. Формирование целевого функционала с регуляризацией
* Для каждого параметра <math>u ∈ U</math>:
** Вычисление невязки: <math>residual = F(y, u) - f</math>
** Вычисление квадратичного члена: <math>fidelity_{term} = residual^2</math>
** Формирование полного функционала: <math>J(u) = fidelity_{term} + f - R</math>
6. Минимизация функционала и определение оптимального решения
* Инициализация минимального значения <math>J_{min} = \infty </math> и индекса <math>min_{index} = 0</math>
* Для каждого <math>u ∈ U</math>:
** Если <math>J(u) < J_{min}</math>:
*** Обновление <math>J_{min} = J(u)</math>
*** Сохранение <math>min_{index} = index(u)</math>
* Определение оптимального решения: <math>x_{opt} = F(x, U[min_{index}])</math>
7. Вывод результатов
* Отображение оптимального значения <math>x_{opt}</math>
* Завершение работы алгоритма

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

1. Вычисление <math> B_{mi}(u) </math>:
* Инициализация, содержит <math>n </math> операций. Сложность - <math>O(n)</math>
* Возведение в степень, содержит <math>n×O(\log m)</math> операций. Сложность - <math>O(n \log m)</math>
* Основной цикл, содержит <math>n × m</math> операций. Сложность - <math>O(n×m)</math>
Итого получили сложность <math>O(n×m)</math>

2. Вычисление метрики:
* Кэш факториалов. Сложность - <math>O(m)</math>
* Регуляризация. Сложность - <math>O(m^2)</math>
* Конечная метрика. Сложность - <math>O(n)</math>
Итого получили сложность <math>O(m^2 + n)</math>

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

=== Информационный граф ===
Информационный граф алгоритма представляет собой направленный ациклический граф (DAG), отражающий поток данных и управляющие зависимости между основными вычислительными модулями (Рис. 2). Граф имеет ярко выраженную ярусно-параллельную структуру, что указывает на наличие значительного ресурса параллелизма, в основном на уровне обработки независимых точек заданной сетки.
[[file:Info_graph_BS_REGULATION.png|thumb|center|500px|Рисунок 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>U</math>. Основные вычислительные ядра алгоритма обладают следующими характеристиками:

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

*Функция metric: Аналогично, расчёт значения метрики в точке<math> u_k </math>зависит только от сглаженного значения в этой же точке и исходных данных, но не от значений метрики в других точках. Это позволяет независимо и параллельно вычислять метрику для разных точек сетки.

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

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

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

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

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

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

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

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

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

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

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

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

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

*Инициализация и финализация GMP-объектов — некоторые операции с числами произвольной точности требуют последовательного доступа к управляющим структурам.
=== Входные и выходные данные алгоритма ===
'''Входные данные''' 
Алгоритм принимает на вход набор экспериментальных данных, представленный в виде пар значений <math>(x_i, y_i)</math>, где <math>i = 1, 2, ..., m</math>. Данные считываются из текстового файла в формате CSV (Comma-Separated Values), где каждая строка содержит два числа, разделенных запятой или пробелом: координату <math>x_i </math> и соответствующее значение <math>y_i</math>. Формат файла допускает использование различных разделителей (запятые, пробелы, переводы строк) благодаря применению регулярных выражений при парсинге.

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

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

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

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

Значение параметра <math>u</math> в диапазоне <math>[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>O(m × n)</math>, где <math>m</math> — количество экспериментальных точек, <math>n</math> — количество значений параметра <math>U</math>

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

5. Вычислительная сложность
*Временная сложность: <math>O(m^2 × n)</math> в наихудшем случае, где <math>m</math> — размер входных данных, <math>n</math> — количество точек параметра <math>U</math>
*Пространственная сложность: <math>O(m + n)</math> для хранения входных данных и промежуточных результатов
*Сложность параллельной части: эффективно распараллеливается как по данным (разные значения <math>U</math>), так и по операциям (внутренние циклы)
== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===
В базовом варианте алгоритм реализуется через последовательное вычисление значений метрики качества для каждого параметра <math>u</math>. Основные вычислительные ядра алгоритма представлены следующими функциями:

'''Функция вычисления базовой суммы (func)''' 
<source lang="cpp">
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;
}
</source>

'''Функция вычисления метрики (metric)''' 
<source lang="cpp">
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;
}
</source>

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

'''Алгоритм быстрого возведения в степень''' 
Для реализации возведения в целую степень используется алгоритм быстрого возведения, основанный на двоичном разложении показателя степени:
<source lang="cpp">
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);
}
</source>
Этот алгоритм имеет сложность <math>O(\log n)</math> вместо <math>O(n)</math> при наивном умножении. Важно отметить, что функция работает с типом mpf_t из библиотеки GMP, что обеспечивает высокую точность вычислений даже для больших степеней.

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

'''1. Кэширование факториалов:''' 
<source lang="cpp">
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;
}
</source>

'''2. Вычисление гамма-функции для целых аргументов через факториал:''' 
<source lang="cpp">
void gamma_cached(mpf_t result, const vector<MPFWrapper> &cache, int n) {
// Γ(n) = (n-1)!
mpf_set(result, cache[n - 1].value);
}
</source>

'''3. Вычисление бета-функции через гамма-функции:''' 
<source lang="cpp">
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);
}
</source>

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

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

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

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

Для интеграции типов GMP с контейнерами STL реализованы специальные обертки:
<source lang="cpp">
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); }
// ... операторы присваивания и приведения типов
};
</source>
Эти обертки обеспечивают корректное управление памятью GMP (инициализация, копирование, освобождение) в контексте C++ RAII (Resource Acquisition Is Initialization).

=== Возможные способы и особенности параллельной реализации алгоритма ===
Как видно из структуры алгоритма, вариантов распараллеливания существует несколько. Алгоритм демонстрирует высокую степень параллелизма, поскольку вычисления для разных значений параметра <math>U</math> являются полностью независимыми. Это позволяет эффективно распределять работу между вычислительными ресурсами без необходимости сложной синхронизации.

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

'''1. Межпроцессный параллелизм (MPI)''' 
Основной массив вычислений разбивается на блоки по значениям параметра <math>U</math>, которые распределяются между MPI-процессами. Каждый процесс работает со своей порцией данных, что минимизирует необходимость в обменах информацией во время вычислений.
<source lang="cpp">
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);
</source>
Такое блочное распределение обеспечивает хорошую балансировку нагрузки, поскольку максимальная разница в количестве элементов между процессами не превышает 1. После завершения вычислений результаты собираются на главном процессе с помощью операции MPI_Gatherv.

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

<source lang="cpp">
#pragma omp parallel for schedule(dynamic)
for (int idx = 0; idx < static_cast<int>(U.size()); ++idx) {
// Вычисления для каждого значения U
}
</source>
Динамическое планирование позволяет адаптироваться к возможной неоднородности вычислений для разных значений параметра <math>U</math>, что улучшает балансировку нагрузки между потоками внутри процесса.

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

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

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

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

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

3. Минимизация критических секций - синхронизация между потоками сведена к минимуму и используется только там, где это действительно необходимо (например, при накоплении суммы).
=== Результаты прогонов ===
Для рассмотрения работы алгоритма для разных случаев вычислений лучше всего подходит решение задачи пролонгации значений. Так алгоритм позволяет пролонгировать значение например 515 на другие подробные сетки и при этом решение обратной задачи никак не изменится, но при этом сильно растёт вычислительная нагрузка и время расчёта. Рассмотрим результат прогонов одной и той же задачи сначала с пролонгацией или точностью, как указано на графике, сначала на <math>10^3</math> точек, потом на <math>10^4</math>, <math>10^5</math> и на <math>10^6</math> отразим зависимость времени вычисления от количества процессоров - Рис. 3
[[file:Execution_time_log_plot.png|thumb|center|500px|Рисунок 3. Зависимость времени выполнения задачи о пролонгации 515 точек на различные точности от числа процессоров в логарифмическом масштабе]]

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

[[file:Speedup_plot_BS_REGULATION.png|thumb|center|500px|Рисунок 4. Зависимость ускорения выполнения задачи о пролонгации 515 точек на различные точности от числа процессоров]]

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

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

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

== Литература ==

<references />

Файл:Execution time log plot.png

2025-12-07T22:21:01Z

Морозов Егор Павлвочи:

Файл:Speedup plot BS REGULATION.png

2025-12-07T22:19:12Z

Морозов Егор Павлвочи:

Файл:Info graph BS REGULATION.png

2025-12-07T21:57:12Z

Морозов Егор Павлвочи:

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

2025-12-06T21:39:39Z

Морозов Егор Павлвочи: /* Особенности реализации последовательного алгоритма */

{{algorithm
| name = Метод решения недетерминированных обратных задач идентификации входных параметров с BS регуляризацией
| serial_complexity = <math>O(nm+m^2)</math>
| pf_height = <math>5</math>
| pf_width = <math>O(n)</math>
| input_data = <math>2m+n+1</math>
| output_data = <math>4n+1</math>
}}

Основные авторы описания: [[Участник:Морозов Егор Павлвочи|Е.П.Морозов]]

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

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

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

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

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

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

В математической литературе длительное время существовала точка зрения, согласно которой всякая математическая задача должна удовлетворять этим требованиям <ref>Курант Р., Вентцель Т. Д. Уравнения с частными производными. – Москва.– : Мир, 1964. – Т. 2. – С. 832.</ref>. Задачи, не удовлетворяющие перечисленным требованиям, называются некорректно поставленными. Следует отметить, что определение некорректно поставленных задач относится к данной паре метрических пространств (F, U), так как в других метриках та же задача может быть корректно поставленной. В рассматриваемом классе задач нарушается устойчивость и однозначность решения. Согласно А.Н. Тихонову, для таких задач существуют квазирешения. Однако чтобы перейти к решению, а поняв его структуру, — к вычислению обратной задачи, необходимо восстановить приблизительную зависимость между входными и выходными параметрами.

Для этого рассмотрим B-сплайны и NURBS. Их полное геометрическое описание и вывод основных свойств были представлены в «NURBS book» <ref>Piegl L., Tiller W. The NURBS book. – Springer Science & Business Media, 2012.</ref>. Помимо их основного применения — создания и сглаживания поверхностей и сложных форм на основе реперных точек, — мы можем применить их для аппроксимации исходных данных и получения системы алгебраических уравнений.

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

Данные системы функциональных уравнений и связанные с ними обратные задачи можно решать различными способами. Пример решения подобных задач для уравнений математической физики представлен в книге <ref>Самарский А. А., Вабищевич П. Н. Численные методы решения обратных задач математической физики. – 2004.</ref>, где они называются обратными задачами идентификации правой части. Однако представленный алгоритм ориентирован не на применение конечно-разностных схем, не на регуляризацию по Тихонову или иные традиционные методы, а на использование нового регуляризатора, основанного на B-сплайнах. При аппроксимации наших данных мы вводили B-сплайны как многочлены порядка m, обозначим их как <math>Q_m(u)</math>, где u - наш новый параметр. Для более точного решения нам необходимо, чтобы данный аппроксимативный многочлен удовлетворял условиям альтернанса Чебышева. Тем самым <math>Q_m(u)</math> будет полиномом наилучшего приближения. Тем самым мы можем за основу взять условия альтернанса
\begin{equation}
F(\vec{x_i})-Q_m(u)=\alpha(-1)^i||{F- Q_m (u)}||,
\end{equation}
где <math>i=1\ldots m+1</math>, <math>\alpha=\pm1</math> одновременно для всех i.
Для обобщения решения будем рассматривать Гильбертовы пространства и в качестве нормы введём интегральную гильбертову норму. А в качестве <math>F(\vec{x_i})</math> рассмотрим наше известное значение для наблюдения, идентификацию начальных параметров для которого мы хотим получить, обозначим его как <math>\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>\vec{\theta} = (\theta_1, \theta_2, \ldots, \theta_m)</math>, где m - количество наблюдений, которые были получены при начальных параметрах <math>\vec{x_i} = (x_{1i}, x_{2i}, \ldots, x_{ni})</math>, где n - количество исходных параметров для рассматриваемой системы, а также дано целевое значение наблюдения, параметры которого мы хотим узнать, обозначим его <math>\Theta, i = 0 \ldots m</math>.

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

[[file:Bspline_regulation_1.png|thumb|center|500px|Рисунок 1. Схема аппроксимативных связей начальных параметров.]]

При этом в виду особенностей B-сплайнов, о которых было упомянуто ранее мы получим, что они зависят от определённых значений параметра u, следовательно узнав при каком значении параметра U будет выполнено наблюдение <math>\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>k = \overline{1,n}</math>, n - количество исходных параметров системы, <math>x_k^*</math> вычисленные параметры системы, при которых возможно наблюдение <math>\Theta</math>

Как писалось в пункте 1.1 для введения регуляризатора нам необходимо, чтобы выполнялся альтернанс Чебышёва, то есть
чтобы многочлен <math>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> F(\vec{x_i}) </math>, необходимо и достаточно существования на <math>\left[a,\ \ b\right]</math>, где <math>a=\min{Q_m\left(u\right)}</math>, <math>b=\max{Q_m\left(u\right)}</math> по крайней мере m + 2 точек <math>x_1<\ldots<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> 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>\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>\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>\theta_i\left(u\right)</math> – i-е значение наблюдения, <math>B_{mi}(u)</math>, <math>\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>B(x,y)</math> - Бета функция, <math>\Gamma(x)</math> - Гамма функция

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

Вычисление Бета-функции в области натуральных чисел, которым и являются <math>i</math> и <math>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>(x, y)</math> из CSV файла
* Разделение на два вектора: <math>\vec x</math> и <math> \vec y</math>
2. Создание сетки параметра u
* Генерация равномерной сетки от 0.01 до 0.99 с шагом 0.001
3. Вычисление полиномов Бернштейна (функция bernstein_polynomial)
* Для каждого u вычисляется сумма:<math> \sum x_i × C_{m-1}^i × u^i × (1-u)^{m-1-i}</math>
* Используется эффективный алгоритм с накоплением
4. Вычисление метрики (функция calculate_metric)
* Вычисление основного члена: <math>(Bernstein(u) - f)^2 + f</math>
* Вычисление регуляризации через бета- и гамма-функции
* Использование кэша факториалов для ускорения
5. Поиск оптимума
* Поиск минимального значения метрики по сетке <math>u</math>
* Вывод оптимального значения <math>x</math>
6. Сохранение результатов
* Запись аппроксимации в output_aprox.dat
* Запись значений метрики в output_metric.dat

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

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

2. Вычисление комбинаторных коэффициентов
* Для каждого элемента сетки <math>u ∈ U</math> и каждого индекса <math>i ∈ [0, m-1]</math>:
** Вычисление биномиального коэффициента <math>C(m, m-i)</math>
** Вычисление степенных членов <math>u^i</math> и <math>(1-u)^{(m-i)}</math>
** Формирование взвешенной суммы элементов вектора <math>\vec x</math>
3. Решение прямой задачи аппроксимации
* Для каждого параметра <math>u ∈ U</math>:
** Инициализация накопительной суммы <math>S = 0</math>
** Для каждого <math>i ∈ [0, m-1]</math>:
*** Вычисление комбинаторного веса: <math>w = C(m, m-i) × u^i × (1-u)^{(m-i)}</math>
***Добавление взвешенного элемента: <math>S = S + w × x[i]</math>
**Сохранение результата <math>F(x, u) = S</math>
4. Вычисление регуляризационного слагаемого
* Инициализация суммы регуляризации <math>R = 0</math>
* Для каждого <math>i ∈ [1, m-1]</math>:
** Инициализация внутренней суммы <math>S_{reg} = 0</math>
** Для каждого <math>j ∈ [0, m-i-1]</math>:
*** Вычисление знака: <math>sign = (-1)^j</math>
*** Вычисление бета-функций: <math>B1 = β(m+i, j+1), B2 = β(m, i)</math>
*** Вычисление гамма-функций: <math>Γ1 = Γ(j+1), Γ2 = Γ(j+1)</math>
*** Формирование слагаемого: <math>term = \frac{sign × B1 × B2}{(Γ1 × Γ2 × (j+i+1))}</math>
*** <math>S_{reg} = S_{reg} + term</math>
** Вычисление весового множителя: <math>weight = \frac{m × (1+i)^2 }{(i × Γ(i)^2)}</math>
** <math>R = R + S_{reg} × y[i] × weight</math>
5. Формирование целевого функционала с регуляризацией
* Для каждого параметра <math>u ∈ U</math>:
** Вычисление невязки: <math>residual = F(y, u) - f</math>
** Вычисление квадратичного члена: <math>fidelity_{term} = residual^2</math>
** Формирование полного функционала: <math>J(u) = fidelity_{term} + f - R</math>
6. Минимизация функционала и определение оптимального решения
* Инициализация минимального значения <math>J_{min} = \infty </math> и индекса <math>min_{index} = 0</math>
* Для каждого <math>u ∈ U</math>:
** Если <math>J(u) < J_{min}</math>:
*** Обновление <math>J_{min} = J(u)</math>
*** Сохранение <math>min_{index} = index(u)</math>
* Определение оптимального решения: <math>x_{opt} = F(x, U[min_{index}])</math>
7. Вывод результатов
* Отображение оптимального значения <math>x_{opt}</math>
* Завершение работы алгоритма

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

1. Вычисление <math> B_{mi}(u) </math>:
* Инициализация, содержит <math>n </math> операций. Сложность - <math>O(n)</math>
* Возведение в степень, содержит <math>n×O(\log m)</math> операций. Сложность - <math>O(n \log m)</math>
* Основной цикл, содержит <math>n × m</math> операций. Сложность - <math>O(n×m)</math>
Итого получили сложность <math>O(n×m)</math>

2. Вычисление метрики:
* Кэш факториалов. Сложность - <math>O(m)</math>
* Регуляризация. Сложность - <math>O(m^2)</math>
* Конечная метрика. Сложность - <math>O(n)</math>
Итого получили сложность <math>O(m^2 + n)</math>

Таким образом финальная сложность алгоритма будет составлять:<math>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>U</math>. Основные вычислительные ядра алгоритма обладают следующими характеристиками:

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

*Функция metric: Аналогично, расчёт значения метрики в точке<math> u_k </math>зависит только от сглаженного значения в этой же точке и исходных данных, но не от значений метрики в других точках. Это позволяет независимо и параллельно вычислять метрику для разных точек сетки.

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

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

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

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

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

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

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

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

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

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

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

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

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

*Инициализация и финализация GMP-объектов — некоторые операции с числами произвольной точности требуют последовательного доступа к управляющим структурам.
=== Входные и выходные данные алгоритма ===
'''Входные данные''' 
Алгоритм принимает на вход набор экспериментальных данных, представленный в виде пар значений <math>(x_i, y_i)</math>, где <math>i = 1, 2, ..., m</math>. Данные считываются из текстового файла в формате CSV (Comma-Separated Values), где каждая строка содержит два числа, разделенных запятой или пробелом: координату <math>x_i </math> и соответствующее значение <math>y_i</math>. Формат файла допускает использование различных разделителей (запятые, пробелы, переводы строк) благодаря применению регулярных выражений при парсинге.

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

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

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

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

Значение параметра <math>u</math> в диапазоне <math>[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>O(m × n)</math>, где <math>m</math> — количество экспериментальных точек, <math>n</math> — количество значений параметра <math>U</math>

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

5. Вычислительная сложность
*Временная сложность: <math>O(m^2 × n)</math> в наихудшем случае, где <math>m</math> — размер входных данных, <math>n</math> — количество точек параметра <math>U</math>
*Пространственная сложность: <math>O(m + n)</math> для хранения входных данных и промежуточных результатов
*Сложность параллельной части: эффективно распараллеливается как по данным (разные значения <math>U</math>), так и по операциям (внутренние циклы)
== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===
В базовом варианте алгоритм реализуется через последовательное вычисление значений метрики качества для каждого параметра <math>u</math>. Основные вычислительные ядра алгоритма представлены следующими функциями:

'''Функция вычисления базовой суммы (func)''' 
<source lang="cpp">
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;
}
</source>

'''Функция вычисления метрики (metric)''' 
<source lang="cpp">
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;
}
</source>

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

'''Алгоритм быстрого возведения в степень''' 
Для реализации возведения в целую степень используется алгоритм быстрого возведения, основанный на двоичном разложении показателя степени:
<source lang="cpp">
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);
}
</source>
Этот алгоритм имеет сложность <math>O(\log n)</math> вместо <math>O(n)</math> при наивном умножении. Важно отметить, что функция работает с типом mpf_t из библиотеки GMP, что обеспечивает высокую точность вычислений даже для больших степеней.

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

'''1. Кэширование факториалов:''' 
<source lang="cpp">
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;
}
</source>

'''2. Вычисление гамма-функции для целых аргументов через факториал:''' 
<source lang="cpp">
void gamma_cached(mpf_t result, const vector<MPFWrapper> &cache, int n) {
// Γ(n) = (n-1)!
mpf_set(result, cache[n - 1].value);
}
</source>

'''3. Вычисление бета-функции через гамма-функции:''' 
<source lang="cpp">
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);
}
</source>

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

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

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

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

Для интеграции типов GMP с контейнерами STL реализованы специальные обертки:
<source lang="cpp">
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); }
// ... операторы присваивания и приведения типов
};
</source>
Эти обертки обеспечивают корректное управление памятью GMP (инициализация, копирование, освобождение) в контексте C++ RAII (Resource Acquisition Is Initialization).

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

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

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

== Литература ==

<references />

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

2025-12-06T21:38:48Z

Морозов Егор Павлвочи: /* Особенности реализации последовательного алгоритма */

{{algorithm
| name = Метод решения недетерминированных обратных задач идентификации входных параметров с BS регуляризацией
| serial_complexity = <math>O(nm+m^2)</math>
| pf_height = <math>5</math>
| pf_width = <math>O(n)</math>
| input_data = <math>2m+n+1</math>
| output_data = <math>4n+1</math>
}}

Основные авторы описания: [[Участник:Морозов Егор Павлвочи|Е.П.Морозов]]

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

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

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

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

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

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

В математической литературе длительное время существовала точка зрения, согласно которой всякая математическая задача должна удовлетворять этим требованиям <ref>Курант Р., Вентцель Т. Д. Уравнения с частными производными. – Москва.– : Мир, 1964. – Т. 2. – С. 832.</ref>. Задачи, не удовлетворяющие перечисленным требованиям, называются некорректно поставленными. Следует отметить, что определение некорректно поставленных задач относится к данной паре метрических пространств (F, U), так как в других метриках та же задача может быть корректно поставленной. В рассматриваемом классе задач нарушается устойчивость и однозначность решения. Согласно А.Н. Тихонову, для таких задач существуют квазирешения. Однако чтобы перейти к решению, а поняв его структуру, — к вычислению обратной задачи, необходимо восстановить приблизительную зависимость между входными и выходными параметрами.

Для этого рассмотрим B-сплайны и NURBS. Их полное геометрическое описание и вывод основных свойств были представлены в «NURBS book» <ref>Piegl L., Tiller W. The NURBS book. – Springer Science & Business Media, 2012.</ref>. Помимо их основного применения — создания и сглаживания поверхностей и сложных форм на основе реперных точек, — мы можем применить их для аппроксимации исходных данных и получения системы алгебраических уравнений.

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

Данные системы функциональных уравнений и связанные с ними обратные задачи можно решать различными способами. Пример решения подобных задач для уравнений математической физики представлен в книге <ref>Самарский А. А., Вабищевич П. Н. Численные методы решения обратных задач математической физики. – 2004.</ref>, где они называются обратными задачами идентификации правой части. Однако представленный алгоритм ориентирован не на применение конечно-разностных схем, не на регуляризацию по Тихонову или иные традиционные методы, а на использование нового регуляризатора, основанного на B-сплайнах. При аппроксимации наших данных мы вводили B-сплайны как многочлены порядка m, обозначим их как <math>Q_m(u)</math>, где u - наш новый параметр. Для более точного решения нам необходимо, чтобы данный аппроксимативный многочлен удовлетворял условиям альтернанса Чебышева. Тем самым <math>Q_m(u)</math> будет полиномом наилучшего приближения. Тем самым мы можем за основу взять условия альтернанса
\begin{equation}
F(\vec{x_i})-Q_m(u)=\alpha(-1)^i||{F- Q_m (u)}||,
\end{equation}
где <math>i=1\ldots m+1</math>, <math>\alpha=\pm1</math> одновременно для всех i.
Для обобщения решения будем рассматривать Гильбертовы пространства и в качестве нормы введём интегральную гильбертову норму. А в качестве <math>F(\vec{x_i})</math> рассмотрим наше известное значение для наблюдения, идентификацию начальных параметров для которого мы хотим получить, обозначим его как <math>\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>\vec{\theta} = (\theta_1, \theta_2, \ldots, \theta_m)</math>, где m - количество наблюдений, которые были получены при начальных параметрах <math>\vec{x_i} = (x_{1i}, x_{2i}, \ldots, x_{ni})</math>, где n - количество исходных параметров для рассматриваемой системы, а также дано целевое значение наблюдения, параметры которого мы хотим узнать, обозначим его <math>\Theta, i = 0 \ldots m</math>.

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

[[file:Bspline_regulation_1.png|thumb|center|500px|Рисунок 1. Схема аппроксимативных связей начальных параметров.]]

При этом в виду особенностей B-сплайнов, о которых было упомянуто ранее мы получим, что они зависят от определённых значений параметра u, следовательно узнав при каком значении параметра U будет выполнено наблюдение <math>\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>k = \overline{1,n}</math>, n - количество исходных параметров системы, <math>x_k^*</math> вычисленные параметры системы, при которых возможно наблюдение <math>\Theta</math>

Как писалось в пункте 1.1 для введения регуляризатора нам необходимо, чтобы выполнялся альтернанс Чебышёва, то есть
чтобы многочлен <math>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> F(\vec{x_i}) </math>, необходимо и достаточно существования на <math>\left[a,\ \ b\right]</math>, где <math>a=\min{Q_m\left(u\right)}</math>, <math>b=\max{Q_m\left(u\right)}</math> по крайней мере m + 2 точек <math>x_1<\ldots<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> 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>\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>\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>\theta_i\left(u\right)</math> – i-е значение наблюдения, <math>B_{mi}(u)</math>, <math>\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>B(x,y)</math> - Бета функция, <math>\Gamma(x)</math> - Гамма функция

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

Вычисление Бета-функции в области натуральных чисел, которым и являются <math>i</math> и <math>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>(x, y)</math> из CSV файла
* Разделение на два вектора: <math>\vec x</math> и <math> \vec y</math>
2. Создание сетки параметра u
* Генерация равномерной сетки от 0.01 до 0.99 с шагом 0.001
3. Вычисление полиномов Бернштейна (функция bernstein_polynomial)
* Для каждого u вычисляется сумма:<math> \sum x_i × C_{m-1}^i × u^i × (1-u)^{m-1-i}</math>
* Используется эффективный алгоритм с накоплением
4. Вычисление метрики (функция calculate_metric)
* Вычисление основного члена: <math>(Bernstein(u) - f)^2 + f</math>
* Вычисление регуляризации через бета- и гамма-функции
* Использование кэша факториалов для ускорения
5. Поиск оптимума
* Поиск минимального значения метрики по сетке <math>u</math>
* Вывод оптимального значения <math>x</math>
6. Сохранение результатов
* Запись аппроксимации в output_aprox.dat
* Запись значений метрики в output_metric.dat

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

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

2. Вычисление комбинаторных коэффициентов
* Для каждого элемента сетки <math>u ∈ U</math> и каждого индекса <math>i ∈ [0, m-1]</math>:
** Вычисление биномиального коэффициента <math>C(m, m-i)</math>
** Вычисление степенных членов <math>u^i</math> и <math>(1-u)^{(m-i)}</math>
** Формирование взвешенной суммы элементов вектора <math>\vec x</math>
3. Решение прямой задачи аппроксимации
* Для каждого параметра <math>u ∈ U</math>:
** Инициализация накопительной суммы <math>S = 0</math>
** Для каждого <math>i ∈ [0, m-1]</math>:
*** Вычисление комбинаторного веса: <math>w = C(m, m-i) × u^i × (1-u)^{(m-i)}</math>
***Добавление взвешенного элемента: <math>S = S + w × x[i]</math>
**Сохранение результата <math>F(x, u) = S</math>
4. Вычисление регуляризационного слагаемого
* Инициализация суммы регуляризации <math>R = 0</math>
* Для каждого <math>i ∈ [1, m-1]</math>:
** Инициализация внутренней суммы <math>S_{reg} = 0</math>
** Для каждого <math>j ∈ [0, m-i-1]</math>:
*** Вычисление знака: <math>sign = (-1)^j</math>
*** Вычисление бета-функций: <math>B1 = β(m+i, j+1), B2 = β(m, i)</math>
*** Вычисление гамма-функций: <math>Γ1 = Γ(j+1), Γ2 = Γ(j+1)</math>
*** Формирование слагаемого: <math>term = \frac{sign × B1 × B2}{(Γ1 × Γ2 × (j+i+1))}</math>
*** <math>S_{reg} = S_{reg} + term</math>
** Вычисление весового множителя: <math>weight = \frac{m × (1+i)^2 }{(i × Γ(i)^2)}</math>
** <math>R = R + S_{reg} × y[i] × weight</math>
5. Формирование целевого функционала с регуляризацией
* Для каждого параметра <math>u ∈ U</math>:
** Вычисление невязки: <math>residual = F(y, u) - f</math>
** Вычисление квадратичного члена: <math>fidelity_{term} = residual^2</math>
** Формирование полного функционала: <math>J(u) = fidelity_{term} + f - R</math>
6. Минимизация функционала и определение оптимального решения
* Инициализация минимального значения <math>J_{min} = \infty </math> и индекса <math>min_{index} = 0</math>
* Для каждого <math>u ∈ U</math>:
** Если <math>J(u) < J_{min}</math>:
*** Обновление <math>J_{min} = J(u)</math>
*** Сохранение <math>min_{index} = index(u)</math>
* Определение оптимального решения: <math>x_{opt} = F(x, U[min_{index}])</math>
7. Вывод результатов
* Отображение оптимального значения <math>x_{opt}</math>
* Завершение работы алгоритма

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

1. Вычисление <math> B_{mi}(u) </math>:
* Инициализация, содержит <math>n </math> операций. Сложность - <math>O(n)</math>
* Возведение в степень, содержит <math>n×O(\log m)</math> операций. Сложность - <math>O(n \log m)</math>
* Основной цикл, содержит <math>n × m</math> операций. Сложность - <math>O(n×m)</math>
Итого получили сложность <math>O(n×m)</math>

2. Вычисление метрики:
* Кэш факториалов. Сложность - <math>O(m)</math>
* Регуляризация. Сложность - <math>O(m^2)</math>
* Конечная метрика. Сложность - <math>O(n)</math>
Итого получили сложность <math>O(m^2 + n)</math>

Таким образом финальная сложность алгоритма будет составлять:<math>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>U</math>. Основные вычислительные ядра алгоритма обладают следующими характеристиками:

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

*Функция metric: Аналогично, расчёт значения метрики в точке<math> u_k </math>зависит только от сглаженного значения в этой же точке и исходных данных, но не от значений метрики в других точках. Это позволяет независимо и параллельно вычислять метрику для разных точек сетки.

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

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

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

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

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

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

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

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

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

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

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

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

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

*Инициализация и финализация GMP-объектов — некоторые операции с числами произвольной точности требуют последовательного доступа к управляющим структурам.
=== Входные и выходные данные алгоритма ===
'''Входные данные''' 
Алгоритм принимает на вход набор экспериментальных данных, представленный в виде пар значений <math>(x_i, y_i)</math>, где <math>i = 1, 2, ..., m</math>. Данные считываются из текстового файла в формате CSV (Comma-Separated Values), где каждая строка содержит два числа, разделенных запятой или пробелом: координату <math>x_i </math> и соответствующее значение <math>y_i</math>. Формат файла допускает использование различных разделителей (запятые, пробелы, переводы строк) благодаря применению регулярных выражений при парсинге.

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

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

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

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

Значение параметра <math>u</math> в диапазоне <math>[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>O(m × n)</math>, где <math>m</math> — количество экспериментальных точек, <math>n</math> — количество значений параметра <math>U</math>

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

5. Вычислительная сложность
*Временная сложность: <math>O(m^2 × n)</math> в наихудшем случае, где <math>m</math> — размер входных данных, <math>n</math> — количество точек параметра <math>U</math>
*Пространственная сложность: <math>O(m + n)</math> для хранения входных данных и промежуточных результатов
*Сложность параллельной части: эффективно распараллеливается как по данным (разные значения <math>U</math>), так и по операциям (внутренние циклы)
== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===
В базовом варианте алгоритм реализуется через последовательное вычисление значений метрики качества для каждого параметра <math>u</math>. Основные вычислительные ядра алгоритма представлены следующими функциями:

'''Функция вычисления базовой суммы (func)''' 
<source lang="cpp">
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;
}
</source>

'''Функция вычисления метрики (metric)''' 
<source lang="cpp">
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;
}
</source>

Помимо основных оптимизированных вычислительных алгоритмов, особое внимание уделялось "подсобным" операциями таким как:
'''Алгоритм быстрого возведения в степень''' 
Для реализации возведения в целую степень используется алгоритм быстрого возведения, основанный на двоичном разложении показателя степени:
<source lang="cpp">
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);
}
</source>
Этот алгоритм имеет сложность <math>O(\log n)</math> вместо <math>O(n)</math> при наивном умножении. Важно отметить, что функция работает с типом mpf_t из библиотеки GMP, что обеспечивает высокую точность вычислений даже для больших степеней.

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

'''1. Кэширование факториалов:''' 
<source lang="cpp">
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;
}
</source>

'''2. Вычисление гамма-функции для целых аргументов через факториал:''' 
<source lang="cpp">
void gamma_cached(mpf_t result, const vector<MPFWrapper> &cache, int n) {
// Γ(n) = (n-1)!
mpf_set(result, cache[n - 1].value);
}
</source>

'''3. Вычисление бета-функции через гамма-функции:''' 
<source lang="cpp">
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);
}
</source>

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

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

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

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

Для интеграции типов GMP с контейнерами STL реализованы специальные обертки:
<source lang="cpp">
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); }
// ... операторы присваивания и приведения типов
};
</source>
Эти обертки обеспечивают корректное управление памятью GMP (инициализация, копирование, освобождение) в контексте C++ RAII (Resource Acquisition Is Initialization).

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

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

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

== Литература ==

<references />

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

2025-12-06T21:38:25Z

Морозов Егор Павлвочи: /* Особенности реализации последовательного алгоритма */

{{algorithm
| name = Метод решения недетерминированных обратных задач идентификации входных параметров с BS регуляризацией
| serial_complexity = <math>O(nm+m^2)</math>
| pf_height = <math>5</math>
| pf_width = <math>O(n)</math>
| input_data = <math>2m+n+1</math>
| output_data = <math>4n+1</math>
}}

Основные авторы описания: [[Участник:Морозов Егор Павлвочи|Е.П.Морозов]]

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

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

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

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

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

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

В математической литературе длительное время существовала точка зрения, согласно которой всякая математическая задача должна удовлетворять этим требованиям <ref>Курант Р., Вентцель Т. Д. Уравнения с частными производными. – Москва.– : Мир, 1964. – Т. 2. – С. 832.</ref>. Задачи, не удовлетворяющие перечисленным требованиям, называются некорректно поставленными. Следует отметить, что определение некорректно поставленных задач относится к данной паре метрических пространств (F, U), так как в других метриках та же задача может быть корректно поставленной. В рассматриваемом классе задач нарушается устойчивость и однозначность решения. Согласно А.Н. Тихонову, для таких задач существуют квазирешения. Однако чтобы перейти к решению, а поняв его структуру, — к вычислению обратной задачи, необходимо восстановить приблизительную зависимость между входными и выходными параметрами.

Для этого рассмотрим B-сплайны и NURBS. Их полное геометрическое описание и вывод основных свойств были представлены в «NURBS book» <ref>Piegl L., Tiller W. The NURBS book. – Springer Science & Business Media, 2012.</ref>. Помимо их основного применения — создания и сглаживания поверхностей и сложных форм на основе реперных точек, — мы можем применить их для аппроксимации исходных данных и получения системы алгебраических уравнений.

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

Данные системы функциональных уравнений и связанные с ними обратные задачи можно решать различными способами. Пример решения подобных задач для уравнений математической физики представлен в книге <ref>Самарский А. А., Вабищевич П. Н. Численные методы решения обратных задач математической физики. – 2004.</ref>, где они называются обратными задачами идентификации правой части. Однако представленный алгоритм ориентирован не на применение конечно-разностных схем, не на регуляризацию по Тихонову или иные традиционные методы, а на использование нового регуляризатора, основанного на B-сплайнах. При аппроксимации наших данных мы вводили B-сплайны как многочлены порядка m, обозначим их как <math>Q_m(u)</math>, где u - наш новый параметр. Для более точного решения нам необходимо, чтобы данный аппроксимативный многочлен удовлетворял условиям альтернанса Чебышева. Тем самым <math>Q_m(u)</math> будет полиномом наилучшего приближения. Тем самым мы можем за основу взять условия альтернанса
\begin{equation}
F(\vec{x_i})-Q_m(u)=\alpha(-1)^i||{F- Q_m (u)}||,
\end{equation}
где <math>i=1\ldots m+1</math>, <math>\alpha=\pm1</math> одновременно для всех i.
Для обобщения решения будем рассматривать Гильбертовы пространства и в качестве нормы введём интегральную гильбертову норму. А в качестве <math>F(\vec{x_i})</math> рассмотрим наше известное значение для наблюдения, идентификацию начальных параметров для которого мы хотим получить, обозначим его как <math>\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>\vec{\theta} = (\theta_1, \theta_2, \ldots, \theta_m)</math>, где m - количество наблюдений, которые были получены при начальных параметрах <math>\vec{x_i} = (x_{1i}, x_{2i}, \ldots, x_{ni})</math>, где n - количество исходных параметров для рассматриваемой системы, а также дано целевое значение наблюдения, параметры которого мы хотим узнать, обозначим его <math>\Theta, i = 0 \ldots m</math>.

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

[[file:Bspline_regulation_1.png|thumb|center|500px|Рисунок 1. Схема аппроксимативных связей начальных параметров.]]

При этом в виду особенностей B-сплайнов, о которых было упомянуто ранее мы получим, что они зависят от определённых значений параметра u, следовательно узнав при каком значении параметра U будет выполнено наблюдение <math>\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>k = \overline{1,n}</math>, n - количество исходных параметров системы, <math>x_k^*</math> вычисленные параметры системы, при которых возможно наблюдение <math>\Theta</math>

Как писалось в пункте 1.1 для введения регуляризатора нам необходимо, чтобы выполнялся альтернанс Чебышёва, то есть
чтобы многочлен <math>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> F(\vec{x_i}) </math>, необходимо и достаточно существования на <math>\left[a,\ \ b\right]</math>, где <math>a=\min{Q_m\left(u\right)}</math>, <math>b=\max{Q_m\left(u\right)}</math> по крайней мере m + 2 точек <math>x_1<\ldots<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> 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>\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>\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>\theta_i\left(u\right)</math> – i-е значение наблюдения, <math>B_{mi}(u)</math>, <math>\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>B(x,y)</math> - Бета функция, <math>\Gamma(x)</math> - Гамма функция

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

Вычисление Бета-функции в области натуральных чисел, которым и являются <math>i</math> и <math>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>(x, y)</math> из CSV файла
* Разделение на два вектора: <math>\vec x</math> и <math> \vec y</math>
2. Создание сетки параметра u
* Генерация равномерной сетки от 0.01 до 0.99 с шагом 0.001
3. Вычисление полиномов Бернштейна (функция bernstein_polynomial)
* Для каждого u вычисляется сумма:<math> \sum x_i × C_{m-1}^i × u^i × (1-u)^{m-1-i}</math>
* Используется эффективный алгоритм с накоплением
4. Вычисление метрики (функция calculate_metric)
* Вычисление основного члена: <math>(Bernstein(u) - f)^2 + f</math>
* Вычисление регуляризации через бета- и гамма-функции
* Использование кэша факториалов для ускорения
5. Поиск оптимума
* Поиск минимального значения метрики по сетке <math>u</math>
* Вывод оптимального значения <math>x</math>
6. Сохранение результатов
* Запись аппроксимации в output_aprox.dat
* Запись значений метрики в output_metric.dat

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

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

2. Вычисление комбинаторных коэффициентов
* Для каждого элемента сетки <math>u ∈ U</math> и каждого индекса <math>i ∈ [0, m-1]</math>:
** Вычисление биномиального коэффициента <math>C(m, m-i)</math>
** Вычисление степенных членов <math>u^i</math> и <math>(1-u)^{(m-i)}</math>
** Формирование взвешенной суммы элементов вектора <math>\vec x</math>
3. Решение прямой задачи аппроксимации
* Для каждого параметра <math>u ∈ U</math>:
** Инициализация накопительной суммы <math>S = 0</math>
** Для каждого <math>i ∈ [0, m-1]</math>:
*** Вычисление комбинаторного веса: <math>w = C(m, m-i) × u^i × (1-u)^{(m-i)}</math>
***Добавление взвешенного элемента: <math>S = S + w × x[i]</math>
**Сохранение результата <math>F(x, u) = S</math>
4. Вычисление регуляризационного слагаемого
* Инициализация суммы регуляризации <math>R = 0</math>
* Для каждого <math>i ∈ [1, m-1]</math>:
** Инициализация внутренней суммы <math>S_{reg} = 0</math>
** Для каждого <math>j ∈ [0, m-i-1]</math>:
*** Вычисление знака: <math>sign = (-1)^j</math>
*** Вычисление бета-функций: <math>B1 = β(m+i, j+1), B2 = β(m, i)</math>
*** Вычисление гамма-функций: <math>Γ1 = Γ(j+1), Γ2 = Γ(j+1)</math>
*** Формирование слагаемого: <math>term = \frac{sign × B1 × B2}{(Γ1 × Γ2 × (j+i+1))}</math>
*** <math>S_{reg} = S_{reg} + term</math>
** Вычисление весового множителя: <math>weight = \frac{m × (1+i)^2 }{(i × Γ(i)^2)}</math>
** <math>R = R + S_{reg} × y[i] × weight</math>
5. Формирование целевого функционала с регуляризацией
* Для каждого параметра <math>u ∈ U</math>:
** Вычисление невязки: <math>residual = F(y, u) - f</math>
** Вычисление квадратичного члена: <math>fidelity_{term} = residual^2</math>
** Формирование полного функционала: <math>J(u) = fidelity_{term} + f - R</math>
6. Минимизация функционала и определение оптимального решения
* Инициализация минимального значения <math>J_{min} = \infty </math> и индекса <math>min_{index} = 0</math>
* Для каждого <math>u ∈ U</math>:
** Если <math>J(u) < J_{min}</math>:
*** Обновление <math>J_{min} = J(u)</math>
*** Сохранение <math>min_{index} = index(u)</math>
* Определение оптимального решения: <math>x_{opt} = F(x, U[min_{index}])</math>
7. Вывод результатов
* Отображение оптимального значения <math>x_{opt}</math>
* Завершение работы алгоритма

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

1. Вычисление <math> B_{mi}(u) </math>:
* Инициализация, содержит <math>n </math> операций. Сложность - <math>O(n)</math>
* Возведение в степень, содержит <math>n×O(\log m)</math> операций. Сложность - <math>O(n \log m)</math>
* Основной цикл, содержит <math>n × m</math> операций. Сложность - <math>O(n×m)</math>
Итого получили сложность <math>O(n×m)</math>

2. Вычисление метрики:
* Кэш факториалов. Сложность - <math>O(m)</math>
* Регуляризация. Сложность - <math>O(m^2)</math>
* Конечная метрика. Сложность - <math>O(n)</math>
Итого получили сложность <math>O(m^2 + n)</math>

Таким образом финальная сложность алгоритма будет составлять:<math>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>U</math>. Основные вычислительные ядра алгоритма обладают следующими характеристиками:

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

*Функция metric: Аналогично, расчёт значения метрики в точке<math> u_k </math>зависит только от сглаженного значения в этой же точке и исходных данных, но не от значений метрики в других точках. Это позволяет независимо и параллельно вычислять метрику для разных точек сетки.

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

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

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

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

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

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

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

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

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

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

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

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

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

*Инициализация и финализация GMP-объектов — некоторые операции с числами произвольной точности требуют последовательного доступа к управляющим структурам.
=== Входные и выходные данные алгоритма ===
'''Входные данные''' 
Алгоритм принимает на вход набор экспериментальных данных, представленный в виде пар значений <math>(x_i, y_i)</math>, где <math>i = 1, 2, ..., m</math>. Данные считываются из текстового файла в формате CSV (Comma-Separated Values), где каждая строка содержит два числа, разделенных запятой или пробелом: координату <math>x_i </math> и соответствующее значение <math>y_i</math>. Формат файла допускает использование различных разделителей (запятые, пробелы, переводы строк) благодаря применению регулярных выражений при парсинге.

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

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

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

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

Значение параметра <math>u</math> в диапазоне <math>[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>O(m × n)</math>, где <math>m</math> — количество экспериментальных точек, <math>n</math> — количество значений параметра <math>U</math>

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

5. Вычислительная сложность
*Временная сложность: <math>O(m^2 × n)</math> в наихудшем случае, где <math>m</math> — размер входных данных, <math>n</math> — количество точек параметра <math>U</math>
*Пространственная сложность: <math>O(m + n)</math> для хранения входных данных и промежуточных результатов
*Сложность параллельной части: эффективно распараллеливается как по данным (разные значения <math>U</math>), так и по операциям (внутренние циклы)
== Программная реализация алгоритма ==

=== Особенности реализации последовательного алгоритма ===
В базовом варианте алгоритм реализуется через последовательное вычисление значений метрики качества для каждого параметра <math>u</math>. Основные вычислительные ядра алгоритма представлены следующими функциями:
'''Функция вычисления базовой суммы (func)''' 
<source lang="cpp">
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;
}
</source>

'''Функция вычисления метрики (metric)''' 
<source lang="cpp">
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;
}
</source>

Помимо основных оптимизированных вычислительных алгоритмов, особое внимание уделялось "подсобным" операциями таким как:
'''Алгоритм быстрого возведения в степень''' 
Для реализации возведения в целую степень используется алгоритм быстрого возведения, основанный на двоичном разложении показателя степени:
<source lang="cpp">
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);
}
</source>
Этот алгоритм имеет сложность <math>O(\log n)</math> вместо <math>O(n)</math> при наивном умножении. Важно отметить, что функция работает с типом mpf_t из библиотеки GMP, что обеспечивает высокую точность вычислений даже для больших степеней.

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

'''1. Кэширование факториалов:''' 
<source lang="cpp">
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;
}
</source>

'''2. Вычисление гамма-функции для целых аргументов через факториал:''' 
<source lang="cpp">
void gamma_cached(mpf_t result, const vector<MPFWrapper> &cache, int n) {
// Γ(n) = (n-1)!
mpf_set(result, cache[n - 1].value);
}
</source>

'''3. Вычисление бета-функции через гамма-функции:''' 
<source lang="cpp">
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);
}
</source>

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

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

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

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

Для интеграции типов GMP с контейнерами STL реализованы специальные обертки:
<source lang="cpp">
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); }
// ... операторы присваивания и приведения типов
};
</source>
Эти обертки обеспечивают корректное управление памятью GMP (инициализация, копирование, освобождение) в контексте C++ RAII (Resource Acquisition Is Initialization).

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

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

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

== Литература ==

<references />

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

2025-12-06T21:20:33Z

Морозов Егор Павлвочи: Добавление некоторых глав

{{algorithm
| name = Метод решения недетерминированных обратных задач идентификации входных параметров с BS регуляризацией
| serial_complexity = <math>O(nm+m^2)</math>
| pf_height = <math>5</math>
| pf_width = <math>O(n)</math>
| input_data = <math>2m+n+1</math>
| output_data = <math>4n+1</math>
}}

Основные авторы описания: [[Участник:Морозов Егор Павлвочи|Е.П.Морозов]]

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

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

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

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

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

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

В математической литературе длительное время существовала точка зрения, согласно которой всякая математическая задача должна удовлетворять этим требованиям <ref>Курант Р., Вентцель Т. Д. Уравнения с частными производными. – Москва.– : Мир, 1964. – Т. 2. – С. 832.</ref>. Задачи, не удовлетворяющие перечисленным требованиям, называются некорректно поставленными. Следует отметить, что определение некорректно поставленных задач относится к данной паре метрических пространств (F, U), так как в других метриках та же задача может быть корректно поставленной. В рассматриваемом классе задач нарушается устойчивость и однозначность решения. Согласно А.Н. Тихонову, для таких задач существуют квазирешения. Однако чтобы перейти к решению, а поняв его структуру, — к вычислению обратной задачи, необходимо восстановить приблизительную зависимость между входными и выходными параметрами.

Для этого рассмотрим B-сплайны и NURBS. Их полное геометрическое описание и вывод основных свойств были представлены в «NURBS book» <ref>Piegl L., Tiller W. The NURBS book. – Springer Science & Business Media, 2012.</ref>. Помимо их основного применения — создания и сглаживания поверхностей и сложных форм на основе реперных точек, — мы можем применить их для аппроксимации исходных данных и получения системы алгебраических уравнений.

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

Данные системы функциональных уравнений и связанные с ними обратные задачи можно решать различными способами. Пример решения подобных задач для уравнений математической физики представлен в книге <ref>Самарский А. А., Вабищевич П. Н. Численные методы решения обратных задач математической физики. – 2004.</ref>, где они называются обратными задачами идентификации правой части. Однако представленный алгоритм ориентирован не на применение конечно-разностных схем, не на регуляризацию по Тихонову или иные традиционные методы, а на использование нового регуляризатора, основанного на B-сплайнах. При аппроксимации наших данных мы вводили B-сплайны как многочлены порядка m, обозначим их как <math>Q_m(u)</math>, где u - наш новый параметр. Для более точного решения нам необходимо, чтобы данный аппроксимативный многочлен удовлетворял условиям альтернанса Чебышева. Тем самым <math>Q_m(u)</math> будет полиномом наилучшего приближения. Тем самым мы можем за основу взять условия альтернанса
\begin{equation}
F(\vec{x_i})-Q_m(u)=\alpha(-1)^i||{F- Q_m (u)}||,
\end{equation}
где <math>i=1\ldots m+1</math>, <math>\alpha=\pm1</math> одновременно для всех i.
Для обобщения решения будем рассматривать Гильбертовы пространства и в качестве нормы введём интегральную гильбертову норму. А в качестве <math>F(\vec{x_i})</math> рассмотрим наше известное значение для наблюдения, идентификацию начальных параметров для которого мы хотим получить, обозначим его как <math>\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>\vec{\theta} = (\theta_1, \theta_2, \ldots, \theta_m)</math>, где m - количество наблюдений, которые были получены при начальных параметрах <math>\vec{x_i} = (x_{1i}, x_{2i}, \ldots, x_{ni})</math>, где n - количество исходных параметров для рассматриваемой системы, а также дано целевое значение наблюдения, параметры которого мы хотим узнать, обозначим его <math>\Theta, i = 0 \ldots m</math>.

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

[[file:Bspline_regulation_1.png|thumb|center|500px|Рисунок 1. Схема аппроксимативных связей начальных параметров.]]

При этом в виду особенностей B-сплайнов, о которых было упомянуто ранее мы получим, что они зависят от определённых значений параметра u, следовательно узнав при каком значении параметра U будет выполнено наблюдение <math>\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>k = \overline{1,n}</math>, n - количество исходных параметров системы, <math>x_k^*</math> вычисленные параметры системы, при которых возможно наблюдение <math>\Theta</math>

Как писалось в пункте 1.1 для введения регуляризатора нам необходимо, чтобы выполнялся альтернанс Чебышёва, то есть
чтобы многочлен <math>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> F(\vec{x_i}) </math>, необходимо и достаточно существования на <math>\left[a,\ \ b\right]</math>, где <math>a=\min{Q_m\left(u\right)}</math>, <math>b=\max{Q_m\left(u\right)}</math> по крайней мере m + 2 точек <math>x_1<\ldots<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> 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>\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>\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>\theta_i\left(u\right)</math> – i-е значение наблюдения, <math>B_{mi}(u)</math>, <math>\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>B(x,y)</math> - Бета функция, <math>\Gamma(x)</math> - Гамма функция

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

Вычисление Бета-функции в области натуральных чисел, которым и являются <math>i</math> и <math>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>(x, y)</math> из CSV файла
* Разделение на два вектора: <math>\vec x</math> и <math> \vec y</math>
2. Создание сетки параметра u
* Генерация равномерной сетки от 0.01 до 0.99 с шагом 0.001
3. Вычисление полиномов Бернштейна (функция bernstein_polynomial)
* Для каждого u вычисляется сумма:<math> \sum x_i × C_{m-1}^i × u^i × (1-u)^{m-1-i}</math>
* Используется эффективный алгоритм с накоплением
4. Вычисление метрики (функция calculate_metric)
* Вычисление основного члена: <math>(Bernstein(u) - f)^2 + f</math>
* Вычисление регуляризации через бета- и гамма-функции
* Использование кэша факториалов для ускорения
5. Поиск оптимума
* Поиск минимального значения метрики по сетке <math>u</math>
* Вывод оптимального значения <math>x</math>
6. Сохранение результатов
* Запись аппроксимации в output_aprox.dat
* Запись значений метрики в output_metric.dat

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

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

2. Вычисление комбинаторных коэффициентов
* Для каждого элемента сетки <math>u ∈ U</math> и каждого индекса <math>i ∈ [0, m-1]</math>:
** Вычисление биномиального коэффициента <math>C(m, m-i)</math>
** Вычисление степенных членов <math>u^i</math> и <math>(1-u)^{(m-i)}</math>
** Формирование взвешенной суммы элементов вектора <math>\vec x</math>
3. Решение прямой задачи аппроксимации
* Для каждого параметра <math>u ∈ U</math>:
** Инициализация накопительной суммы <math>S = 0</math>
** Для каждого <math>i ∈ [0, m-1]</math>:
*** Вычисление комбинаторного веса: <math>w = C(m, m-i) × u^i × (1-u)^{(m-i)}</math>
***Добавление взвешенного элемента: <math>S = S + w × x[i]</math>
**Сохранение результата <math>F(x, u) = S</math>
4. Вычисление регуляризационного слагаемого
* Инициализация суммы регуляризации <math>R = 0</math>
* Для каждого <math>i ∈ [1, m-1]</math>:
** Инициализация внутренней суммы <math>S_{reg} = 0</math>
** Для каждого <math>j ∈ [0, m-i-1]</math>:
*** Вычисление знака: <math>sign = (-1)^j</math>
*** Вычисление бета-функций: <math>B1 = β(m+i, j+1), B2 = β(m, i)</math>
*** Вычисление гамма-функций: <math>Γ1 = Γ(j+1), Γ2 = Γ(j+1)</math>
*** Формирование слагаемого: <math>term = \frac{sign × B1 × B2}{(Γ1 × Γ2 × (j+i+1))}</math>
*** <math>S_{reg} = S_{reg} + term</math>
** Вычисление весового множителя: <math>weight = \frac{m × (1+i)^2 }{(i × Γ(i)^2)}</math>
** <math>R = R + S_{reg} × y[i] × weight</math>
5. Формирование целевого функционала с регуляризацией
* Для каждого параметра <math>u ∈ U</math>:
** Вычисление невязки: <math>residual = F(y, u) - f</math>
** Вычисление квадратичного члена: <math>fidelity_{term} = residual^2</math>
** Формирование полного функционала: <math>J(u) = fidelity_{term} + f - R</math>
6. Минимизация функционала и определение оптимального решения
* Инициализация минимального значения <math>J_{min} = \infty </math> и индекса <math>min_{index} = 0</math>
* Для каждого <math>u ∈ U</math>:
** Если <math>J(u) < J_{min}</math>:
*** Обновление <math>J_{min} = J(u)</math>
*** Сохранение <math>min_{index} = index(u)</math>
* Определение оптимального решения: <math>x_{opt} = F(x, U[min_{index}])</math>
7. Вывод результатов
* Отображение оптимального значения <math>x_{opt}</math>
* Завершение работы алгоритма

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

1. Вычисление <math> B_{mi}(u) </math>:
* Инициализация, содержит <math>n </math> операций. Сложность - <math>O(n)</math>
* Возведение в степень, содержит <math>n×O(\log m)</math> операций. Сложность - <math>O(n \log m)</math>
* Основной цикл, содержит <math>n × m</math> операций. Сложность - <math>O(n×m)</math>
Итого получили сложность <math>O(n×m)</math>

2. Вычисление метрики:
* Кэш факториалов. Сложность - <math>O(m)</math>
* Регуляризация. Сложность - <math>O(m^2)</math>
* Конечная метрика. Сложность - <math>O(n)</math>
Итого получили сложность <math>O(m^2 + n)</math>

Таким образом финальная сложность алгоритма будет составлять:<math>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>U</math>. Основные вычислительные ядра алгоритма обладают следующими характеристиками:

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

*Функция metric: Аналогично, расчёт значения метрики в точке<math> u_k </math>зависит только от сглаженного значения в этой же точке и исходных данных, но не от значений метрики в других точках. Это позволяет независимо и параллельно вычислять метрику для разных точек сетки.

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

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

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

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

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

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

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

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

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

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

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

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

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

*Инициализация и финализация GMP-объектов — некоторые операции с числами произвольной точности требуют последовательного доступа к управляющим структурам.
=== Входные и выходные данные алгоритма ===
'''Входные данные''' 
Алгоритм принимает на вход набор экспериментальных данных, представленный в виде пар значений <math>(x_i, y_i)</math>, где <math>i = 1, 2, ..., m</math>. Данные считываются из текстового файла в формате CSV (Comma-Separated Values), где каждая строка содержит два числа, разделенных запятой или пробелом: координату <math>x_i </math> и соответствующее значение <math>y_i</math>. Формат файла допускает использование различных разделителей (запятые, пробелы, переводы строк) благодаря применению регулярных выражений при парсинге.

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

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

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

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

Значение параметра <math>u</math> в диапазоне <math>[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>O(m × n)</math>, где <math>m</math> — количество экспериментальных точек, <math>n</math> — количество значений параметра <math>U</math>

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

5. Вычислительная сложность
*Временная сложность: <math>O(m^2 × n)</math> в наихудшем случае, где <math>m</math> — размер входных данных, <math>n</math> — количество точек параметра <math>U</math>
*Пространственная сложность: <math>O(m + n)</math> для хранения входных данных и промежуточных результатов
*Сложность параллельной части: эффективно распараллеливается как по данным (разные значения <math>U</math>), так и по операциям (внутренние циклы)
== Программная реализация алгоритма ==

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

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

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

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

== Литература ==

<references />

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

2025-12-06T20:34:51Z

Морозов Егор Павлвочи:

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

2025-12-06T20:14:39Z

Морозов Егор Павлвочи:

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

2025-12-06T20:07:24Z

Морозов Егор Павлвочи: /* Вычислительное ядро алгоритма */

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

2025-12-06T20:05:13Z

Морозов Егор Павлвочи: /* Схема реализации последовательного алгоритма */

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

2025-12-02T07:13:02Z

Морозов Егор Павлвочи:

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

2025-11-20T16:03:52Z

Морозов Егор Павлвочи: /* Вычислительное ядро алгоритма */

{{algorithm
| name = Метод решения недетерминированных обратных задач идентификации входных параметров с BS регуляризацией
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n (n + 1)}{2}</math>
| output_data = <math>\frac{n (n + 1)}{2}</math>
}}

Основные авторы описания: [[Участник:Морозов Егор Павлвочи|Е.П.Морозов]]

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

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

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

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

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

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

В математической литературе длительное время существовала точка зрения, согласно которой всякая математическая задача должна удовлетворять этим требованиям <ref>Курант Р., Вентцель Т. Д. Уравнения с частными производными. – Москва.– : Мир, 1964. – Т. 2. – С. 832.</ref>. Задачи, не удовлетворяющие перечисленным требованиям, называются некорректно поставленными. Следует отметить, что определение некорректно поставленных задач относится к данной паре метрических пространств (F, U), так как в других метриках та же задача может быть корректно поставленной. В рассматриваемом классе задач нарушается устойчивость и однозначность решения. Согласно А.Н. Тихонову, для таких задач существуют квазирешения. Однако чтобы перейти к решению, а поняв его структуру, — к вычислению обратной задачи, необходимо восстановить приблизительную зависимость между входными и выходными параметрами.

Для этого рассмотрим B-сплайны и NURBS. Их полное геометрическое описание и вывод основных свойств были представлены в «NURBS book» <ref>Piegl L., Tiller W. The NURBS book. – Springer Science & Business Media, 2012.</ref>. Помимо их основного применения — создания и сглаживания поверхностей и сложных форм на основе реперных точек, — мы можем применить их для аппроксимации исходных данных и получения системы алгебраических уравнений.

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

Данные системы функциональных уравнений и связанные с ними обратные задачи можно решать различными способами. Пример решения подобных задач для уравнений математической физики представлен в книге <ref>Самарский А. А., Вабищевич П. Н. Численные методы решения обратных задач математической физики. – 2004.</ref>, где они называются обратными задачами идентификации правой части. Однако представленный алгоритм ориентирован не на применение конечно-разностных схем, не на регуляризацию по Тихонову или иные традиционные методы, а на использование нового регуляризатора, основанного на B-сплайнах. При аппроксимации наших данных мы вводили B-сплайны как многочлены порядка m, обозначим их как <math>Q_m(u)</math>, где u - наш новый параметр. Для более точного решения нам необходимо, чтобы данный аппроксимативный многочлен удовлетворял условиям альтернанса Чебышева. Тем самым <math>Q_m(u)</math> будет полиномом наилучшего приближения. Тем самым мы можем за основу взять условия альтернанса
\begin{equation}
F(\vec{x_i})-Q_m(u)=\alpha(-1)^i||{F- Q_m (u)}||,
\end{equation}
где <math>i=1\ldots m+1</math>, <math>\alpha=\pm1</math> одновременно для всех i.
Для обобщения решения будем рассматривать Гильбертовы пространства и в качестве нормы введём интегральную гильбертову норму. А в качестве <math>F(\vec{x_i})</math> рассмотрим наше известное значение для наблюдения, идентификацию начальных параметров для которого мы хотим получить, обозначим его как <math>\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>\vec{\theta} = (\theta_1, \theta_2, \ldots, \theta_m)</math>, где m - количество наблюдений, которые были получены при начальных параметрах <math>\vec{x_i} = (x_{1i}, x_{2i}, \ldots, x_{ni})</math>, где n - количество исходных параметров для рассматриваемой системы, а также дано целевое значение наблюдения, параметры которого мы хотим узнать, обозначим его <math>\Theta, i = 0 \ldots m</math>.

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

[[file:Bspline_regulation_1.png|thumb|center|500px|Рисунок 1. Схема аппроксимативных связей начальных параметров.]]

При этом в виду особенностей B-сплайнов, о которых было упомянуто ранее мы получим, что они зависят от определённых значений параметра u, следовательно узнав при каком значении параметра U будет выполнено наблюдение <math>\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>k = \overline{1,n}</math>, n - количество исходных параметров системы, <math>x_k^*</math> вычисленные параметры системы, при которых возможно наблюдение <math>\Theta</math>

Как писалось в пункте 1.1 для введения регуляризатора нам необходимо, чтобы выполнялся альтернанс Чебышёва, то есть
чтобы многочлен <math>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> F(\vec{x_i}) </math>, необходимо и достаточно существования на <math>\left[a,\ \ b\right]</math>, где <math>a=\min{Q_m\left(u\right)}</math>, <math>b=\max{Q_m\left(u\right)}</math> по крайней мере m + 2 точек <math>x_1<\ldots<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> 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>\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>\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>\theta_i\left(u\right)</math> – i-е значение наблюдения, <math>B_{mi}(u)</math>, <math>\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>B(x,y)</math> - Бета функция, <math>\Gamma(x)</math> - Гамма функция

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

Вычисление Бета-функции в области натуральных чисел, которым и являются <math>i<\math> и <math>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}
Таким образом мы свели наш алгоритм к вычислению множества сочетай, факториалов и возведений в степень, что и будет вычислительным ядром нашего алгоритма.

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

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

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

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

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

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

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

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

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

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

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

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

== Литература ==

<references />

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

2025-11-18T07:51:10Z

Морозов Егор Павлвочи: Внесено вычислительное ядро алгоритма

{{algorithm
| name = Метод решения недетерминированных обратных задач идентификации входных параметров с BS регуляризацией
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n (n + 1)}{2}</math>
| output_data = <math>\frac{n (n + 1)}{2}</math>
}}

Основные авторы описания: [[Участник:Морозов Егор Павлвочи|Е.П.Морозов]]

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

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

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

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

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

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

В математической литературе длительное время существовала точка зрения, согласно которой всякая математическая задача должна удовлетворять этим требованиям <ref>Курант Р., Вентцель Т. Д. Уравнения с частными производными. – Москва.– : Мир, 1964. – Т. 2. – С. 832.</ref>. Задачи, не удовлетворяющие перечисленным требованиям, называются некорректно поставленными. Следует отметить, что определение некорректно поставленных задач относится к данной паре метрических пространств (F, U), так как в других метриках та же задача может быть корректно поставленной. В рассматриваемом классе задач нарушается устойчивость и однозначность решения. Согласно А.Н. Тихонову, для таких задач существуют квазирешения. Однако чтобы перейти к решению, а поняв его структуру, — к вычислению обратной задачи, необходимо восстановить приблизительную зависимость между входными и выходными параметрами.

Для этого рассмотрим B-сплайны и NURBS. Их полное геометрическое описание и вывод основных свойств были представлены в «NURBS book» <ref>Piegl L., Tiller W. The NURBS book. – Springer Science & Business Media, 2012.</ref>. Помимо их основного применения — создания и сглаживания поверхностей и сложных форм на основе реперных точек, — мы можем применить их для аппроксимации исходных данных и получения системы алгебраических уравнений.

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

Данные системы функциональных уравнений и связанные с ними обратные задачи можно решать различными способами. Пример решения подобных задач для уравнений математической физики представлен в книге <ref>Самарский А. А., Вабищевич П. Н. Численные методы решения обратных задач математической физики. – 2004.</ref>, где они называются обратными задачами идентификации правой части. Однако представленный алгоритм ориентирован не на применение конечно-разностных схем, не на регуляризацию по Тихонову или иные традиционные методы, а на использование нового регуляризатора, основанного на B-сплайнах. При аппроксимации наших данных мы вводили B-сплайны как многочлены порядка m, обозначим их как <math>Q_m(u)</math>, где u - наш новый параметр. Для более точного решения нам необходимо, чтобы данный аппроксимативный многочлен удовлетворял условиям альтернанса Чебышева. Тем самым <math>Q_m(u)</math> будет полиномом наилучшего приближения. Тем самым мы можем за основу взять условия альтернанса
\begin{equation}
F(\vec{x_i})-Q_m(u)=\alpha(-1)^i||{F- Q_m (u)}||,
\end{equation}
где <math>i=1\ldots m+1</math>, <math>\alpha=\pm1</math> одновременно для всех i.
Для обобщения решения будем рассматривать Гильбертовы пространства и в качестве нормы введём интегральную гильбертову норму. А в качестве <math>F(\vec{x_i})</math> рассмотрим наше известное значение для наблюдения, идентификацию начальных параметров для которого мы хотим получить, обозначим его как <math>\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>\vec{\theta} = (\theta_1, \theta_2, \ldots, \theta_m)</math>, где m - количество наблюдений, которые были получены при начальных параметрах <math>\vec{x_i} = (x_{1i}, x_{2i}, \ldots, x_{ni})</math>, где n - количество исходных параметров для рассматриваемой системы, а также дано целевое значение наблюдения, параметры которого мы хотим узнать, обозначим его <math>\Theta, i = 0 \ldots m</math>.

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

[[file:Bspline_regulation_1.png|thumb|center|500px|Рисунок 1. Схема аппроксимативных связей начальных параметров.]]

При этом в виду особенностей B-сплайнов, о которых было упомянуто ранее мы получим, что они зависят от определённых значений параметра u, следовательно узнав при каком значении параметра U будет выполнено наблюдение <math>\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>k = \overline{1,n}</math>, n - количество исходных параметров системы, <math>x_k^*</math> вычисленные параметры системы, при которых возможно наблюдение <math>\Theta</math>

Как писалось в пункте 1.1 для введения регуляризатора нам необходимо, чтобы выполнялся альтернанс Чебышёва, то есть
чтобы многочлен <math>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> F(\vec{x_i}) </math>, необходимо и достаточно существования на <math>\left[a,\ \ b\right]</math>, где <math>a=\min{Q_m\left(u\right)}</math>, <math>b=\max{Q_m\left(u\right)}</math> по крайней мере m + 2 точек <math>x_1<\ldots<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> 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>\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>\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>\theta_i\left(u\right)</math> – i-е значение наблюдения, <math>B_{mi}(u)</math>, <math>\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>B(x,y)</math> - Бета функция, <math>\Gamma(x)</math> - Гамма функция

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

Вычисление Бета-функции в области натуральных чисел, которым и являются <math>i<\math> и <math>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}
Таким образом мы свели наш алгоритм к вычислению множества сочетай, факториалов и возведений в степень, что и будет вычислительным ядром нашего алгоритма.
=== Макроструктура алгоритма ===

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

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

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

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

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

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

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

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

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

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

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

== Литература ==

<references />

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

2025-10-28T20:58:21Z

Морозов Егор Павлвочи: Новая страница: «{{algorithm | name = Метод решения недетерминированных обратных задач идентификации ...»

{{algorithm
| name = Метод решения недетерминированных обратных задач идентификации входных параметров с BS регуляризацией
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n (n + 1)}{2}</math>
| output_data = <math>\frac{n (n + 1)}{2}</math>
}}

Основные авторы описания: [[Участник:Морозов Егор Павлвочи|Е.П.Морозов]]

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

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

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

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

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

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

В математической литературе длительное время существовала точка зрения, согласно которой всякая математическая задача должна удовлетворять этим требованиям <ref>Курант Р., Вентцель Т. Д. Уравнения с частными производными. – Москва.– : Мир, 1964. – Т. 2. – С. 832.</ref>. Задачи, не удовлетворяющие перечисленным требованиям, называются некорректно поставленными. Следует отметить, что определение некорректно поставленных задач относится к данной паре метрических пространств (F, U), так как в других метриках та же задача может быть корректно поставленной. В рассматриваемом классе задач нарушается устойчивость и однозначность решения. Согласно А.Н. Тихонову, для таких задач существуют квазирешения. Однако чтобы перейти к решению, а поняв его структуру, — к вычислению обратной задачи, необходимо восстановить приблизительную зависимость между входными и выходными параметрами.

Для этого рассмотрим B-сплайны и NURBS. Их полное геометрическое описание и вывод основных свойств были представлены в «NURBS book» <ref>Piegl L., Tiller W. The NURBS book. – Springer Science & Business Media, 2012.</ref>. Помимо их основного применения — создания и сглаживания поверхностей и сложных форм на основе реперных точек, — мы можем применить их для аппроксимации исходных данных и получения системы алгебраических уравнений.

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

Данные системы функциональных уравнений и связанные с ними обратные задачи можно решать различными способами. Пример решения подобных задач для уравнений математической физики представлен в книге <ref>Самарский А. А., Вабищевич П. Н. Численные методы решения обратных задач математической физики. – 2004.</ref>, где они называются обратными задачами идентификации правой части. Однако представленный алгоритм ориентирован не на применение конечно-разностных схем, не на регуляризацию по Тихонову или иные традиционные методы, а на использование нового регуляризатора, основанного на B-сплайнах. При аппроксимации наших данных мы вводили B-сплайны как многочлены порядка m, обозначим их как <math>Q_m(u)</math>, где u - наш новый параметр. Для более точного решения нам необходимо, чтобы данный аппроксимативный многочлен удовлетворял условиям альтернанса Чебышева. Тем самым <math>Q_m(u)</math> будет полиномом наилучшего приближения. Тем самым мы можем за основу взять условия альтернанса
\begin{equation}
F(\vec{x_i})-Q_m(u)=\alpha(-1)^i||{F- Q_m (u)}||,
\end{equation}
где <math>i=1\ldots m+1</math>, <math>\alpha=\pm1</math> одновременно для всех i.
Для обобщения решения будем рассматривать Гильбертовы пространства и в качестве нормы введём интегральную гильбертову норму. А в качестве <math>F(\vec{x_i})</math> рассмотрим наше известное значение для наблюдения, идентификацию начальных параметров для которого мы хотим получить, обозначим его как <math>\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>\vec{\theta} = (\theta_1, \theta_2, \ldots, \theta_m)</math>, где m - количество наблюдений, которые были получены при начальных параметрах <math>\vec{x_i} = (x_{1i}, x_{2i}, \ldots, x_{ni})</math>, где n - количество исходных параметров для рассматриваемой системы, а также дано целевое значение наблюдения, параметры которого мы хотим узнать, обозначим его <math>\Theta, i = 0 \ldots m</math>.

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

[[file:Bspline_regulation_1.png|thumb|center|500px|Рисунок 1. Схема аппроксимативных связей начальных параметров.]]

При этом в виду особенностей B-сплайнов, о которых было упомянуто ранее мы получим, что они зависят от определённых значений параметра u, следовательно узнав при каком значении параметра U будет выполнено наблюдение <math>\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>k = \overline{1,n}</math>, n - количество исходных параметров системы, <math>x_k^*</math> вычисленные параметры системы, при которых возможно наблюдение <math>\Theta</math>

Как писалось в пункте 1.1 для введения регуляризатора нам необходимо, чтобы выполнялся альтернанс Чебышёва, то есть
чтобы многочлен <math>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> F(\vec{x_i}) </math>, необходимо и достаточно существования на <math>\left[a,\ \ b\right]</math>, где <math>a=\min{Q_m\left(u\right)}</math>, <math>b=\max{Q_m\left(u\right)}</math> по крайней мере m + 2 точек <math>x_1<\ldots<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> 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>\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>\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>\theta_i\left(u\right)</math> – i-е значение наблюдения, <math>B_{mi}(u)</math>, <math>\delta_m</math>– расчётные коэффициенты, которые вычисляются согласно формулам

\begin{equation}
B_{mi}(u) = \sum_{i = 0}^{m}\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>B(x,y)</math> - Бета функция, <math>\Gamma(x)</math> - Гамма функция

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

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

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

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

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

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

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

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

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

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

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

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

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

== Литература ==

<references />

Файл:Bspline regulation 1.png

2025-10-28T20:29:24Z

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

== Краткое описание ==
Схема аппроксимативных связей начальных параметров для метода решения недетерминированных обратных задач идентификации входных параметров с BS регуляризацией