МГУ имени М.В.Ломоносова
Московский государственный университет
имени М.В.Ломоносова
 

 ЧИСЛЕННЫЕ МЕТОДЫ ОБРАБОТКИ 


Experimental data
 
 


background

    Все новости.


      новое - от  01.05.2017 :

Численные методы в практике экспериментатора

           Текст заказчика:

"Восстановление профилей лазерных пучков с большим количеством разномасштабных экстремумов"

           Решается система линейных уравнений AX=B c прямоугольной матрицей A (NIZ, NDIM), NIZ< NDIM. Здесь NIZ – число полосовых измерений интенсивности пучка c несимметричным профилем, расположенного в прямоугольной области [-x, x] * [-y, y]. Число точек сетки по каждому направлению: N по оси X, M по оси Y. NDIM = N*M – число неизвестных. В файлах с данными для каждого примера содержатся: матрица системы А, записанная по столбцам, правая часть B размером NIZ, а также точное решение системы X, записанное в виде матрицы размером N*M по столбцам.
           Пример 1. 16 экстремумов. Область [-8,8]*[-8,8]. Число точек сетки по каждому направлению N=M=81. Шаг сетки hx=hy=0.2 . Размер матрицы системы: NIZ=1621, NDIM= 6561. Относительная погрешность восстановления SVD-методом .

Рис.1. Вид точного решения (16-ть экстремумов).

Скачать матрицу А - здесь (5мб, здесь и далее формат zip).
Скачать правую часть b - здесь .

Пример 2. 25 экстремумов. Область [-10,10]*[-10,10]. Число точек сетки по каждому направлению N=M=101. Шаг сетки hx=hy=0.2 . Размер матрицы системы: NIZ=2561, NDIM= 10201. Относительная погрешность восстановления SVD-методом .

Рис.2. Вид точного решения (25-ть экстремумов).

Скачать матрицу А - здесь (8.33мб).
Скачать правую часть b - здесь .

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

(конец текста заказчика)


           Ну, насколько мы поняли в лаборатории - решается интегральное 2-у мерное уравнение Фредгольма I-го рода. Известно ядро интегрального оператора по которому насчитывалась после дискретизации уравнения матрица А; правая часть b вначале модельная, потом из эксперимента. Для нас заказчик моделировал точное решение(см. рис. 1,2) и по матрице А насчитывал правые части. Таким образом, в наличии: прямоугольная недоопределенная матрица(средних размеров - 1621*6561, 2561*10201) и правая часть. Требуется построить быстрый алгоритм нахождения решения с возможностью реализации кода в приборе. Визуально было определено, что матрицы разрежены.

           У заказчика работающие алгоритм/программа в наличии уже были (использовался метод решения по SVD), но их не устраивало время решения - порядка 1,5 часов. Для практики требовалось - считанные секунды (не более минуты). Очевидно, открытым оставался и вопрос выбора границы отсечения малых сингулярных чисел. Так например, в Примере 1 (где 16 экстремумов) эта граница у них была найдена как 0.034 , а в Примере 2 - 0.047 .

           Перед решением задачи было определено, что матрица не является матрицей специальной струкуры(ленточной/теплицевой). Анализ на блочность структуры провести не удалось. Информацию о ядре заказчик не предоставлял.

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

           Найти вектор минимальной нормы минимизирующий невязку ||Ах — b||.

Такой вектор существует и единственен; известны теоремы:
Теорема 1. Рассмотрим СЛАУ Ах=b. Величина ||Ах — b|| минимальна для таких x что Ах = АА+b, а среди всех таких x наименьшую величину ||x|| имеет х0 = А+b, где А+ - обобщенная обратная(псевдообратная) матрица.
Теорема 2(Мур — Пенроуз). Для любой матрицы А существует единственная обобщенная обратная матрица А+.

Укажем способы вычисления псевдообратной матрицы для нашего случая(недооопределенной матрицы) :
1. Если матрица А имеет ранг равный количеству строк(т.е. все строки нашей матрицы линейно-независимы), то А+T(AАT)-1
Но, выяснилось, что rank A < m, таким образом, данным представлением мы воспользоваться не могли. Рассмотрим следующую возможность.
2. Стандартный путь получения псевдообратной матрицы это использование SVD: Если A = UΣVT — сингулярное разложение, то A+ = VΣ+UT.
Посмотрим, что получается в этом случае. Отметим, что решение получали без отсечения малых сингулярных чисел. Результаты в Табл.1.

          

Табл.1. Данные примера 2 - матрица 2561*10201.
   заказчик  наша лаборатория
 ||A||*  37.0389
 ||b||  257.2287
 метод решения  SVD с отсечением сингулярных
чисел меньше 0.047
 SVD без отсечения
 ||Ax||  257.0518  257.2287
 ||Ax-b||  0.7252  0.00026
 100%||Ax-b||/||b||  0.282 %   0.0001 %
 время решения  1 h. 30 min.  1 h. 2 min. 23.60 sec.
процессор 2.94GHz
*евклидовы нормы.

           Сверить/повторить результаты Табл.1 можно по решению заказчика(скачать здесь) и по нашему решению(здесь).
На рис.3 дан вид полученного решения - негладкая функция с количеством пиков более чем 25.

Рис.3. Вид полученного решения.


Уменьшение времени с 1.5 часов до 1часа удалось за счет использования операций с разреженной матрицей. Но, тем не менее, результаты оказались по временному фактору совершенно неудовлетворительны (собственно, как и у заказчика) - т.к. требовалось получение решения менее, чем за 1 мин.

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

           Выяснилось:
1. методы регуляризации(с одним параметром регуляризации a или с 3-мя: при нормах решения, первой и второй производных решения) давали более-менее приемлимое приближение 100%||Ax-b||/||b||=0.3%, но с неприемлимо большим временем - порядка 1часа.
Полученное регуляризованное решение при a=0.01 можно скачать здесь.
2. комбинация метода регуляризации с итерационным методом Крейга дал уже существенно более лучшие показатели: время порядка 1мин.
Полученное решение можно скачать здесь.

           Таким образом, стало ясно, что необходимо искать принципиально отличный подход от принципа решения обратной задачи СЛАУ с большой размерностью искомого вектора(NDIM= 10201). Таким подходом мы предложили:

           Моделирование решения двумерными пиками и далее поиск параметров пиков из решения задачи минимизации функционала невязки ||Ax-b||.

Каждый пик моделировался асимметричным модифицированным гауссианом (не классическим) - всего 5-ть параметров на каждый пик. Суммарно получается существенно меньше, чем 10201 искомых значений.
Но, при таком подходе, есть принципиальный момент: пользователю априори неизвестно количество пиков. Предполагать изначально, что их м.б. миллион и потом отсекать те пики, у которых получились малые(на уровне погрешностей) амплитуды(высоты) просто нереально по численке. Этот эффект еще называют "корытом": часть искомых параметров практически не влияют на изменения функционала, образно говоря - "гуляют по дну корыта". Фактически это или линейно-зависимые параметры модели, или параметры, которые не влияют на модель по физике. И определить их практически невозможно, тем более, при многоэкстремальности функционала. Рассматривалось также и такое предложение: увеличивать количество пиков по одному. Не прошло из-за лимита времени.

          Таким образом, вначале решалась задача определения количества пиков, а далее уже использовался алгоритм поиска глобального минимума. Алгоритм/метод решения первой задачи - определения количества пиков, - мы не приводим всилу того, что разработанный алгоритм является эвристическим, без доказательной базы(алгоритм подобен тому, который мы используем при поиске глобального минимума: переход из одного локального максимума в другой). Тем не менее, идейка сработала. И если у вас есть, например, 2-у мерная поверхность, имеющая множество пиков разной высоты и разных полуширин (по сечениям параллельно координатам), но с неизвестными координатами центров пиков - можете присылать свой массив. Алгоритм/программа выдаст координаты центров пиков, их высоты и полуширины. При необходимости, можно посчитать и параметры аппроксимирующего эллипса на полувысоте. Да, данный алгоритм в принципе должен срабатывать не только для 2-у мерной поверхности, но и для произвольной.

           Результаты счета даны в Табл.2.

          

Табл.2. Счет данных Примеров 1,2.
 x - решение  Пример 1  Пример 2
 ||Ax-b||  0.007  0.016
 100%||Ax-b||/||b||   0.005 %   0.006 %
 Ввод данных(сек.)   0.86   1.67
 Формирование портретов
матриц A, АT.
Перемножение матриц A*АT,
формирование ее портрета (сек.)
  1.72   4.35
 Поиск количества пиков (сек.)   1.73   1.58
 Глобальная минимизация (сек.)   5.58   16.43
 Общее время (сек.)
(время вывода не учтено)
  9.91   24.01
 Решение скачать   здесь   здесь

   
 Решение для Примера 1.  Решение для Примера 2.

          Заключение.
          Алгоритм/программа может использоваться если априори известно, что решение представляет собой суперпозицию пиков. При этом, моделирование решения пиками и дальнейшая глобальная минимизация позволила избежать проблемных моментов связанных с неполным рангом недоопределенной алгебраической системы, возникающей после дискретизации 2-у мерного уравнения Фредгольма I-го рода.





opening 15.11.2009    © math-lab.ru    All rights reserved.

  Яндекс цитирования
  Rambler's Top100
 
  Яндекс.Метрика
  Locations of visitors to this page