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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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