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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

• Инициализация параметрической сетки [math]\displaystyle{ U = [0.01, 0.02, ..., 0.99] }[/math]

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

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

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

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

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

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

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

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

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

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

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

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