Участник:Kiseliov/Метод регуляризации Тихонова: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показано 27 промежуточных версий этого же участника)
Строка 22: Строка 22:
  
 
== Математическое описание алгоритма ==
 
== Математическое описание алгоритма ==
Пусть нам даны спектры искажённого изображения <math>\bar u_{ij}</math> и ядра <math>\bar K_{ij}</math> - это квадратные матрицы, состоящие из комплексных чисел, размера <math>N\timesN</math>. Тогда спектр исходного изображения в позиции <math>(i, j)</math> находится по формуле:
+
Пусть нам даны спектры искажённого изображения <math>\bar u_{ij}</math> и ядра <math>\bar K_{ij}</math> - это квадратные матрицы, состоящие из комплексных чисел, размера <math>N\times N</math>. Тогда спектр исходного изображения в позиции <math>(i, j)</math> находится по формуле:
  
 
<math>\bar z_{ij} =\frac{(\bar K_{ij})^*\bar u_{ij} (\vartriangle x)^2}{\bigl|\bar K_{ij}\bigl|^2(\vartriangle x)^4+\alpha \bigl(\bigl(\frac{\pi}{2R}\bigl)^2(i^2+j^2)\bigl)^r}</math><br>
 
<math>\bar z_{ij} =\frac{(\bar K_{ij})^*\bar u_{ij} (\vartriangle x)^2}{\bigl|\bar K_{ij}\bigl|^2(\vartriangle x)^4+\alpha \bigl(\bigl(\frac{\pi}{2R}\bigl)^2(i^2+j^2)\bigl)^r}</math><br>
Строка 28: Строка 28:
 
Здесь <math>\vartriangle x=\frac{4R}{N}</math>, где R - один из параметров ядра (известная, наперёд заданная величина).<math>(\bar K_{ij})^*</math> - сопряжение спектра ядра. <math>\alpha</math> и <math>r</math> - параметры метода.
 
Здесь <math>\vartriangle x=\frac{4R}{N}</math>, где R - один из параметров ядра (известная, наперёд заданная величина).<math>(\bar K_{ij})^*</math> - сопряжение спектра ядра. <math>\alpha</math> и <math>r</math> - параметры метода.
  
В  реализации будем считать <math>\alpha = 10000</math> и <math>r = 1</math> (наиболее подходящие параметры в нашем случае).
+
В  реализации будем считать <math>\alpha = 0.5</math> и <math>r = 1</math> (наиболее подходящие параметры в нашем случае).
  
 
== Вычислительное ядро алгоритма ==
 
== Вычислительное ядро алгоритма ==
В данном методе вычислительное ядро совпадает с алгоритмом, так как в позиции <math>(i, j)</math> вычисления производятся независимо по формуле, указанной пункте 1.2.
+
В данном методе вычислительным ядром являются вычисления по формуле из пункта 1.2 в позиции <math>(i, j)</math>, поскольку формула применяется к позиции <math>(i, j)</math> независимо от других точек.
  
 
Вычислительная сложность формулы - 3 умножения в числителе, 9 умножений, 2 сложения в знаменателе и одно общее деление (возведение в степень не считаем, так как положили <math>r = 1</math>). Итого, 15 операций. Формулы применяется независимо в каждой точке <math>(i, j)</math> <math>\Rightarrow</math> вычислительная сложность ядра алгоритма - 15 операций.
 
Вычислительная сложность формулы - 3 умножения в числителе, 9 умножений, 2 сложения в знаменателе и одно общее деление (возведение в степень не считаем, так как положили <math>r = 1</math>). Итого, 15 операций. Формулы применяется независимо в каждой точке <math>(i, j)</math> <math>\Rightarrow</math> вычислительная сложность ядра алгоритма - 15 операций.
Строка 39: Строка 39:
  
 
== Схема реализации последовательного алгоритма ==
 
== Схема реализации последовательного алгоритма ==
 +
Реализация указанного алгоритма на псевдокоде в стиле C:
 +
 +
<syntaxhighlight lang="C">
 +
tikhonov_regularization(N, *image, *kernel, R, alpha, r, *restored)
 +
{
 +
    // ДельтаX в формуле (логично посчитать эти значения единожды)
 +
    delta_x = (4 * R) / N;
 +
    delta_x_pow2 = delta_x * delta_x;
 +
    delta_x_pow4 = delta_x_pow2 * delta_x_pow2;
 +
 +
    // Это выражение - константа, и его логично посчитать заранее. CONSTANT_PI - число пи
 +
    constant = (CONSTANT_PI / (2 * R)) * (CONSTANT_PI / (2 * R));
 +
 +
    // Цикл по позициям (i, j)
 +
    for (i = 0; i < N; i++)
 +
        for (j = 0; j < N; j++)
 +
        {
 +
            // Предполагается, что матрицы хранятся в виде 1D массива, поэтому нужно вычислить правильный индекс
 +
            index = i * N + j;
 +
 +
            // Мы сдвигаем центр системы координат в центр изображения
 +
            m = i - N / 2;
 +
            n = N / 2 - j;
 +
           
 +
            // Вычисление формулы. Предполагается, что pow возводит число в степень, conj возвращает сопряжённое комплексное число,
 +
            // .real() и .imag() - получение действительной или мнимой частей комплексного числа соответственно.
 +
            upper_part = conj(kernel[index]) * image[index] * delta_x_pow2;
 +
            kernel_squared_module = kernel[index].real() * kernel[index].real() + kernel[index].imag() * kernel[index].imag();
 +
            down_part1 = kernel_squared_module * delta_x_pow4;
 +
            m_pow2 = m * m;
 +
            n_pow2 = n * n;
 +
            down_part2 = alpha * pow(constant * (m_pow2 + n_pow2), r);
 +
            down_part = down_part1 + down_part2;
 +
            restored[index] = upper_part / down_part;
 +
        }
 +
}
 +
</syntaxhighlight>
  
 
== Последовательная сложность алгоритма ==
 
== Последовательная сложность алгоритма ==
Строка 44: Строка 81:
  
 
== Информационный граф ==
 
== Информационный граф ==
 +
[[Файл:TikhonovRegularizationGraph3.png|center|thumb|'''Рисунок 1.''' Информационный граф метода регуляризации Тихонова. Вершины означают применение формулы из пункта 2 к <math>(i, j)</math> позиции входа для получения <math>(i, j)</math> позиции выхода. Рёбра в графе отсутствуют в связи с отсутствием связей по данным.]]
  
 
== Ресурс параллелизма алгоритма ==
 
== Ресурс параллелизма алгоритма ==
Поскольку в каждой позиции <math>(i, j)</math> вычисления производятся независимо, то ресурс параллелизма будет равен размеру матриц, то есть <math>N^2</math> (здесь учтено <math>M=N</math>).
+
Поскольку в каждой позиции <math>(i, j)</math> вычисления производятся независимо, то ресурс параллелизма будет равен размеру матриц, то есть <math>N^2</math>.
  
 
== Входные и выходные данные алгоритма ==
 
== Входные и выходные данные алгоритма ==
'''Входные данные:''' Спектр искажённого изображения и спектр ядра свёртки. Обычно дано искажённое изображение (от которого берётся преобразование Фурье). Мы считаем параметры ядра известными. Оно генерируется, и затем от него берётся преобразование Фурье (получается спектр).
+
'''Входные данные:''' Две матрицы размера <math>N\times N</math> - спектр искажённого изображения и спектр ядра свёртки. Обычно дано искажённое изображение (от которого берётся преобразование Фурье). Мы считаем параметры ядра известными. Оно генерируется, и затем от него берётся преобразование Фурье (получается спектр).
 +
 
 +
'''Выходные данные:''' Матрица размера <math>N\times N</math> - спектр исходного изображения. Спектр можно преобразовать в изображение с помощью обратного преобразования Фурье.
  
'''Выходные данные:''' Спектр исходного изображения. Спектр можно преобразовать в изображение с помощью обратного преобразования Фурье.
 
 
== Свойства алгоритма ==
 
== Свойства алгоритма ==
 
= ЧАСТЬ. Программная реализация алгоритма =
 
= ЧАСТЬ. Программная реализация алгоритма =
Строка 58: Строка 97:
 
== Возможные способы и особенности параллельной реализации алгоритма ==
 
== Возможные способы и особенности параллельной реализации алгоритма ==
 
== Масштабируемость алгоритма и его реализации ==
 
== Масштабируемость алгоритма и его реализации ==
 +
Все запуски программы проводились на с/к Ломоносов 2. Для распараллеливания использовалась технология OpenMP.
 +
 +
{| class="wikitable sortable"
 +
|+ Характеристики системы
 +
|-
 +
! Характеристика || Значение
 +
|-
 +
| Пиковая производительность || 5,505 Пфлопс
 +
|-
 +
| Производительность на тесте Linpack || 2,478 Пфлопс
 +
|-
 +
| Число вычислительных узлов || 1 730
 +
|-
 +
| Основной тип процессора || Intel Haswell-EP E5-2697v3, 2.6 GHz, 14 cores
 +
|-
 +
| Тип ускорителя || NVidia Tesla K40M
 +
|-
 +
| Общее число ядер || 64 384
 +
|-
 +
| Оперативная память на узел || 64 GB
 +
|-
 +
| Проектная мощность || 10 МВт
 +
|-
 +
| Операционная система || CentOS 7
 +
|}
 +
 +
Конфигурация системы в разделе compute (в котором производились запуски программы):
 +
{| class="wikitable sortable"
 +
|+ 1487 узлов со следующими характеристиками
 +
|-
 +
! Характеристика || Значение
 +
|-
 +
| Объём памяти на GPU || 11.56 GB
 +
|-
 +
| Центральный процессор || Intel Haswell-EP E5-2697v3, 2.6 GHz, 14 cores
 +
|-
 +
| Графический ускоритель || NVidia Tesla K40M
 +
|-
 +
| ОЗУ || 64 GB
 +
|}
 +
 +
Перейдём к графику времени работы алгоритма. На графике по оси Ox, Oy откладывается длина одной стороны изображения в пикселях (т.е. это - число <math>N</math>, если всё изображение - матрица размера <math>N\times N</math>) и количество нитей OpenMP (нити запускаются параллельно, если доступно необходимое количество ресурсов). По оси Oz откладывается время работы в секундах для каждой точки <math>(x, y)</math>. Временем работы здесь является время работы функции, применяющей метод регуляризации Тихонова к изображению и ядру и заполняющей выходную матрицу. Время чтения/записи файлов, а также время выполнения прочих служебных функций не учитывалось.
 +
Количество нитей задаётся следующей последовательностью: <math>1, 4, 8, 12, 16, 20, ..., 128</math>. А длина стороны изображения в пикселях задаётся такой последовательностью: <math>128, 128+16, ..., 256, 256+32, ..., 4096</math>
 +
[[Файл:TikhonovRegularizationTimes4.png|center|'''Рисунок 2.''' Время выполнения параллельной реализации метода регуляризации Тихонова в секундах.]]
  
 
== Динамические характеристики и эффективность реализации алгоритма ==
 
== Динамические характеристики и эффективность реализации алгоритма ==
 
== Выводы для классов архитектур ==
 
== Выводы для классов архитектур ==
 
== Существующие реализации алгоритма ==
 
== Существующие реализации алгоритма ==
 +
Автору такие не известны. Поиск в открытых источниках результата не дал.
  
 
= Литература =
 
= Литература =

Текущая версия на 13:38, 3 ноября 2022

Автор: Киселёв Е. И.

1 ЧАСТЬ. Свойства и структура алгоритмов

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

Метод регуляризации Тихонова заключается в следующем:

Нам дано некоторое искажённое изображение (в нашем случае мы рассматриваем статические аберрации). Фактически, некоторое изображение искажается при помощи свёртки с так называемым ядром. То есть, мы имеем уравнение Фредгольма первого рода типа свертки вида:

[math]K \circledast z =\textstyle\int\limits_{-\infty}^{\infty}\textstyle\int\limits_{-\infty}^{\infty}K(x_1-s_1,x_2-s_2 )z(s_1,s_2 )ds_1\,ds_2 = u(x_1,x_2 ), -\infty\lt x_1,x_2\lt \infty[/math]

Здесь [math]K(x_1,x_2 )∈L_2 (\mathbb{R}^2)[/math] – аппаратная функция прибора (ядро), [math]u(x_1,x_2 )∈L_2 (\mathbb{R}^2)[/math] – искаженное изображение, а [math]z(x_1,x_2 )[/math] – искомое реконструируемое изображение.

Наша задача - восстановить исходное изображение, зная параметры ядра. Метод регуляризации Тихонова говорит о том, что решение имеет вид:

[math]z(x_1 ,x_2 ) = \frac{1}{4\pi^2}\textstyle\int\limits_{-\infty}^{\infty}\textstyle\int\limits_{-\infty}^{\infty}\frac{\bar K (-\omega_1,-\omega_2)\bar u (\omega_1,\omega_2)}{|\bar K (\omega_1,\omega_2 )|^2+\alpha M(\omega_1,\omega_2)}e^{i(\omega_1x_1 + \omega_2x_2)}d\omega_1\,d\omega_2[/math]

Здесь [math]\bar K (x_1,x_2 )[/math] – спектр ядра, [math]\bar u (x_1,x_2 )[/math] – спектр искаженного изображения, а [math]M(\omega_1,\omega_2)[/math] – заданная четная функция, обладающая следующими свойствами:

  1. [math]M(\omega_1,\omega_2)[/math] кусочно-непрерывна в любой конечной области
  2. [math]M(\omega_1,\omega_2)[/math] неотрицательна: [math]M(0,0)\ge0[/math] и [math]M(\omega_1,\omega_2)\gt 0[/math] при [math]\omega_1,\omega_2\neq0[/math]
  3. для достаточно больших [math]|\omega_1|,|\omega_2| \Rightarrow M(\omega_1,\omega_2)\ge C\gt 0[/math]
  4. для [math]\forall \alpha \gt 0 \Rightarrow \frac{\bar K (-\omega_1,-\omega_2)}{|\bar K (\omega_1,\omega_2 )|^2+\alpha M(\omega_1,\omega_2)}∈L_2 (\mathbb{R}^2)[/math]

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

Пусть нам даны спектры искажённого изображения [math]\bar u_{ij}[/math] и ядра [math]\bar K_{ij}[/math] - это квадратные матрицы, состоящие из комплексных чисел, размера [math]N\times N[/math]. Тогда спектр исходного изображения в позиции [math](i, j)[/math] находится по формуле:

[math]\bar z_{ij} =\frac{(\bar K_{ij})^*\bar u_{ij} (\vartriangle x)^2}{\bigl|\bar K_{ij}\bigl|^2(\vartriangle x)^4+\alpha \bigl(\bigl(\frac{\pi}{2R}\bigl)^2(i^2+j^2)\bigl)^r}[/math]

Здесь [math]\vartriangle x=\frac{4R}{N}[/math], где R - один из параметров ядра (известная, наперёд заданная величина).[math](\bar K_{ij})^*[/math] - сопряжение спектра ядра. [math]\alpha[/math] и [math]r[/math] - параметры метода.

В реализации будем считать [math]\alpha = 0.5[/math] и [math]r = 1[/math] (наиболее подходящие параметры в нашем случае).

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

В данном методе вычислительным ядром являются вычисления по формуле из пункта 1.2 в позиции [math](i, j)[/math], поскольку формула применяется к позиции [math](i, j)[/math] независимо от других точек.

Вычислительная сложность формулы - 3 умножения в числителе, 9 умножений, 2 сложения в знаменателе и одно общее деление (возведение в степень не считаем, так как положили [math]r = 1[/math]). Итого, 15 операций. Формулы применяется независимо в каждой точке [math](i, j)[/math] [math]\Rightarrow[/math] вычислительная сложность ядра алгоритма - 15 операций.

1.4 Макроструктура алгоритма

Макроструктура алгоритма представляет собой проход по всем точкам [math][/math](i, j)[math][/math] и вычисление [math][/math](i, j)[math][/math] позиции результирующей матрицы по формуле, указанной в пункте 1.2.

1.5 Схема реализации последовательного алгоритма

Реализация указанного алгоритма на псевдокоде в стиле C:

tikhonov_regularization(N, *image, *kernel, R, alpha, r, *restored)
{
    // ДельтаX в формуле (логично посчитать эти значения единожды)
    delta_x = (4 * R) / N;
    delta_x_pow2 = delta_x * delta_x;
    delta_x_pow4 = delta_x_pow2 * delta_x_pow2;

    // Это выражение - константа, и его логично посчитать заранее. CONSTANT_PI - число пи
    constant = (CONSTANT_PI / (2 * R)) * (CONSTANT_PI / (2 * R));

    // Цикл по позициям (i, j)
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
        {
            // Предполагается, что матрицы хранятся в виде 1D массива, поэтому нужно вычислить правильный индекс
            index = i * N + j;

            // Мы сдвигаем центр системы координат в центр изображения
            m = i - N / 2;
            n = N / 2 - j;
            
            // Вычисление формулы. Предполагается, что pow возводит число в степень, conj возвращает сопряжённое комплексное число,
            // .real() и .imag() - получение действительной или мнимой частей комплексного числа соответственно.
            upper_part = conj(kernel[index]) * image[index] * delta_x_pow2;
            kernel_squared_module = kernel[index].real() * kernel[index].real() + kernel[index].imag() * kernel[index].imag();
            down_part1 = kernel_squared_module * delta_x_pow4;
            m_pow2 = m * m;
            n_pow2 = n * n;
            down_part2 = alpha * pow(constant * (m_pow2 + n_pow2), r);
            down_part = down_part1 + down_part2;
            restored[index] = upper_part / down_part;
        }
}

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

В пункте 1.3 мы получили, что число операций, необходимых для вычисления формулы, равно [math]15[/math]. Всего формула будет применена [math]N^2[/math] раз. Итого, в последовательном варианте сложность составит [math]15 * N^2[/math] операций.

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

Рисунок 1. Информационный граф метода регуляризации Тихонова. Вершины означают применение формулы из пункта 2 к [math](i, j)[/math] позиции входа для получения [math](i, j)[/math] позиции выхода. Рёбра в графе отсутствуют в связи с отсутствием связей по данным.

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

Поскольку в каждой позиции [math](i, j)[/math] вычисления производятся независимо, то ресурс параллелизма будет равен размеру матриц, то есть [math]N^2[/math].

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

Входные данные: Две матрицы размера [math]N\times N[/math] - спектр искажённого изображения и спектр ядра свёртки. Обычно дано искажённое изображение (от которого берётся преобразование Фурье). Мы считаем параметры ядра известными. Оно генерируется, и затем от него берётся преобразование Фурье (получается спектр).

Выходные данные: Матрица размера [math]N\times N[/math] - спектр исходного изображения. Спектр можно преобразовать в изображение с помощью обратного преобразования Фурье.

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

2 ЧАСТЬ. Программная реализация алгоритма

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

2.2 Локальность данных и вычислений

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

2.4 Масштабируемость алгоритма и его реализации

Все запуски программы проводились на с/к Ломоносов 2. Для распараллеливания использовалась технология OpenMP.

Характеристики системы
Характеристика Значение
Пиковая производительность 5,505 Пфлопс
Производительность на тесте Linpack 2,478 Пфлопс
Число вычислительных узлов 1 730
Основной тип процессора Intel Haswell-EP E5-2697v3, 2.6 GHz, 14 cores
Тип ускорителя NVidia Tesla K40M
Общее число ядер 64 384
Оперативная память на узел 64 GB
Проектная мощность 10 МВт
Операционная система CentOS 7

Конфигурация системы в разделе compute (в котором производились запуски программы):

1487 узлов со следующими характеристиками
Характеристика Значение
Объём памяти на GPU 11.56 GB
Центральный процессор Intel Haswell-EP E5-2697v3, 2.6 GHz, 14 cores
Графический ускоритель NVidia Tesla K40M
ОЗУ 64 GB

Перейдём к графику времени работы алгоритма. На графике по оси Ox, Oy откладывается длина одной стороны изображения в пикселях (т.е. это - число [math]N[/math], если всё изображение - матрица размера [math]N\times N[/math]) и количество нитей OpenMP (нити запускаются параллельно, если доступно необходимое количество ресурсов). По оси Oz откладывается время работы в секундах для каждой точки [math](x, y)[/math]. Временем работы здесь является время работы функции, применяющей метод регуляризации Тихонова к изображению и ядру и заполняющей выходную матрицу. Время чтения/записи файлов, а также время выполнения прочих служебных функций не учитывалось. Количество нитей задаётся следующей последовательностью: [math]1, 4, 8, 12, 16, 20, ..., 128[/math]. А длина стороны изображения в пикселях задаётся такой последовательностью: [math]128, 128+16, ..., 256, 256+32, ..., 4096[/math]

Рисунок 2. Время выполнения параллельной реализации метода регуляризации Тихонова в секундах.

2.5 Динамические характеристики и эффективность реализации алгоритма

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

2.7 Существующие реализации алгоритма

Автору такие не известны. Поиск в открытых источниках результата не дал.

3 Литература

  • Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. – М.: Наука, 1979.
  • Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Численные методы решения некорректных задач. – М.: Наука, 1990.
  • Гудмен Дж. Введение в фурье‐оптику. – М.: Мир, 1970. 364 с.