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

{{algorithm
| name              = Метод решения недетерминированных <br>обратных задач идентификации <br>входных параметров <br>с 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. Инициализация и ввод данных
*	Загрузка экспериментальных данных из 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>
*	Завершение работы алгоритма

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

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

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

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

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

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

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

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

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

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



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

<references />