Регуляризация

Рассмотрим  метод регуляризации, предложенный Тихоновым для решения обратных задач на примере линейной задачи. Как известно, она ставится в форме решения плохо обусловленной системы М линейных уравнений:
(1)
с неизвестным вектором N^J, подлежащим определению (по повторяющимся индексам здесь и далее мы предполагаем суммирование). Алгоритм построения матрицы А и конкретный вид вектора правых частей уравнений f выражают физическую постановку задачи. В частности, задача томографии состоит в численном решении системы большого числа интегральных уравнений, дискретизация который приводит к  системе линейных уравнений с матрицей А и правыми частями f, выражающими результаты измерений.

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

Общеупотребительный подход к решению некорректных задач заключается в их регуляризации. При этом используется дополнительная априорная информация о решении, которая может быть как качественной, так и количественной. Например можно искать решение, максимально близкое к некоторому профилю , т.е. к некоторому вектору  N₀. Концепция регуляризации применительно к (1) сводится к замене задачи на поиск псевдорешения на задачу о минимизации следующего функционала:
, (2)
где  a - малый положительный параметр регуляризации, который необходимо подобрать определенным способом. 
Задача о минимизации функционала (2) для системы линейных алгебраических уравнений (1) эквивалентна системе, содержащей M линейных уравнений: 
, (3)
где символ E обозначает единичную матрицу размера MхM. Решая систему (3), можно получить регуляризованное решение, зависящее от a. Из (3) хорошо ясен смысл параметра a: при малых a~0 обусловленность системы (3) близка к плохой обусловленности (1), а при больших a система (3) обусловлена хорошо, но её решение далеко от решения исходной обратной задачи. А именно, чем больше параметр регуляризации, тем ближе решение к N₀. Очевидно, что на практике необходимо выбирать промежуточные a. Принцип выбора конкретного значения параметра регуляризации будет обсуждён на следующей странице при анализе результатов численного эксперимента.

Описанный регуляризационный подход сводит некорректную задачу  к условно-корректной (по Тихонову) задаче отыскания решения системы (3), которое, в силу линейности задачи, является единственным и устойчивым. Заметим, что метод регуляризации успешно применяется и для решения нелинейных некорректных обратных задач. Однако в этом случае заменить задачу минимизации функционала (2) линейной системой типа (3) невозможно. С вычислительной точки зрения, это означает, что вместо задачи решения системы алгебраических линейных уравнений придется численно решать задачу на глобальный экстремум функции, которая часто более сложна (в частности, из-за возможного наличия в нелинейном случае многих локальных минимумов).

Компьютерное моделирование.
Проиллюстрируем сказанное конкретным вычислительным примером. Зададим некоторую модель ионосферы (рис.1), которая образована наложением на регулярный высотный профиль неоднородности, расположенной в центре, и промоделируем процесс доплеровских измерений. Высота ИСЗ в компьютерном эксперименте составляла 900 км, два приёмника располагались в точках с координатами 0 км и 500 км по поверхности Земли в плоскости пролёта ИСЗ. Дискретизация интегральных уравнений проводилась, в соответствии с [1], методом кусочно-планарной аппроксимации на сетке 21х21 узел. Такая крупная (по сравнению с реальными задачами томографии ионосферы) сетка была взята нами для того, чтобы проверить выдвигаемую методику с помощью точного метода решения (5), не прибегая к приближённым итерационным алгоритмам. Поскольку число арифметических операций, требующихся компьютеру для завершения точного метода, порядка (MK)6 , то для расчёта на сетке 100х100 узлов необходимо время в 15,000 раз большее. В скобках заметим, что на выполнение программы решения системы (5) методом Гаусса с выбором главного элемента [3] занимает на Pentium-200 время порядка минуты, поэтому на практике следует, конечно, использовать итерационные методы. Эффективным может оказаться, например, метод сопряжённых градиентов, поскольку в его рамках можно организовать эффективный спуск по параметру a.
Задав реалистичный профиль N0(z), проиллюстрируем предложенный подход, рассчитав точное решение регуляризованной задачи. Для того, чтобы определить оптимальный параметр регуляризации, проведём серию расчётов с различными a, контролируя невязку. На рис.2 изображена зависимость нормы невязки 
(6) 
решения регуляризованной задачи от параметра регуляризации a. На том же рисунке изображена суммарная ошибка (т.е. норма отклонения решения от модельного распределения).
Для реконструкции можно использовать значение a, соответствующее глобальному минимуму зависимости e(a) (рис.3), либо применить т.н. принцип невязки [6], который требует выбора a, с которым невязка приблизительно равна сумме погрешностей измерений (т.е. заданий правой части) и аппроксимации (результат реконструкции изображён на рис.4). Из сравнения четырёх полей видно, что наилучшее совпадение с моделью демонстрирует последний рисунок.
В заключение заметим, что применяемые на практике методы реконструкции неявно также используют элементы регуляризации (для метода ART3 параметрами регуляризации являются число итераций и индексы релаксации, графики, подобные рис.3., приведены в [5]). Однако такую регуляризацию по-прежнему нельзя считать доказательством устойчивости решения. 
Дальнейшие компьютерные эксперименты и сопоставление с результатами других алгоритмов, видимо, позволит выделить круг задач, в которых метод регуляризации может быть наиболее эффективным по сравнению с остальными.