CC BY
30ОБ АЛГОРИТМАХ ОТЫСКАНИЯ НАИМЕНЬШЕГО РАССТОЯНИЯ МЕЖДУ СИСТЕМАМИ ТОЧЕК
В ПРОСТРАНСТВЕ
И, А. Блатов1, Е, В, Китаева2
1 Поволжский государственный университет телекоммуникация и информатики, 443010, Самара 2 Самарский национальный исследовательский университет имени академика С.П. Королева, 443086, Самара
УДК 517.968
Б01: 10.24411/9999-016А-2019-10009
Рассматривается задача определения минимума сумм квадратов расстояний между двумя системами точек в пространстве. Необходимость решения такой задачи возникает в компьютерной химии и при распознавании образов. Грубый алгоритм имеет факториальную сложность. Предлагается алгоритм, позволяющий найти решение с заданной точностью е за 0(те3/е3/2) арифметических операций.
Ключевые слова: распознавание образов, совмещение объектов, венгерский алгоритм, алгоритм Кабша.
Введение
Пусть в пространстве заданы две системы точек М и Ж, каждая го которых содержит п точек. Рассматривается задача отыскания такого положения этих систем в пространстве и паросочетания их вершин, для которых сумма квадратов попарных расстояний р(М, N) между точками Ми N была бы минимальной. В случае фиксированного паросочетания данная задача решается за линейное (0(п)) количество операций алгоритмом Кабша [1, 2]. Грубый поиск оптимального паросочетания требует перебора п! различных положений. В случае фиксированного положения в пространстве систем и известен венгерский алгоритм [3] отыскания оптимального паросочетания за 0(п3) операций. Однако неизвестен [4] алгоритм общего назначения, позволяющий за полиномиальное время найти оптимальное паросочетапие и соответствующее р(М, N). Различные варианты алгоритмов, сочетающих метод Кабша и венгерский метод, рассматривались в [5-7] (см. также библиографию в [6,7]). Однако это либо эвристические алгоритмы, не гарантирующие отыскание глобального минимума, либо алгоритмы, использующие взаимозависимость и специальные свойства сравниваемых структур, либо экспоненциально или факториально сложные алгоритмы. Авторы [7]) утверждают, что предложенный ими алгоритм (¡о-РКИ.\П)!ЯТ имеет полиномиальную сложность, но не приводят точных оценок вычислительной сложности. В настоящей статье для решения данной задачи рассматриваются три алгоритма общего назначения, основанные на алгоритмах Кабша и венгерском алгоритме, позволяющие найти решение поставленной задачи с заданной точностью 0(п3 /е3/2) за арифметических операций. Проведен сравнительный экспериментальный анализ этих алгоритмов. Первый (грубый) алгоритм представляет собой сочетание венгерского алгоритма с перебором по узлам е-сети в пространстве параметров изометрий в К3. Второй алгоритм (комбинированный) отличается от первого последовательным выполнением венгерского алгоритма и алгоритма Кабша вместо одного лишь венгерского алгоритма. Третий алгоритм (итерационный) сочетает в себе венгерский алгоритм, алгоритм Кабша и итерационный метод покоординатного спуска в пространстве параметров изометрий в К3. Проведен сравнительный экспериментальный анализ этих алгоритмов. Показана высокая эффективность второго алгоритма по сравнению с первым и третьим.
Расчет проводился на семействах атомных решеток из баз данных МНИЦТМ (Самара).
!ЯВ.\ 978-5-901548-42-4
1 Постановка задачи и алгоритмы
Пусть М = {Mi(xi,yi, Zi), 1 < i < п}, N = {Ni(xi,yi, Zi), 1 < i < n} — две совокупности точек (радиус-векторов) в R3 ,заданных координатами относительно прямоугольной декартовой системы координат OXYZ. Через I обозначим единичную матрицу третьего порядка. Если w = (w1,w2,w3) вектор в R3, то || w || =
■\Jw\ + и}"2 + w¡2 обозначает его евклидову порму. через С,С1,С2 • • • будем обозначать положительные константы, не зависящие от п. Если это те вызывает недоразумений, то одним и тем же символом С могут обозначаться разные константы. Для некоторой величины f зависящей от п, запись f = О(п) означает, что f < Сп.
Рассмотрим задачу отыскания вектора v = (v\, v2, v3), ортогональной матрицы третьего порядка U = {uij} (UUт = I) и перестановки Р = (pi,p2, • • • ,Рп) первых п натуральных чисел, для которых величина Е = ^ II UMPi — Ni — V ||2 минимальна. Систему М будем называть подвижной системой, а, N —
неподвижной системой. Сделаем следующие предположения.
Al. Центры тяжести — и — ^ систем М и N совмещены сдвигом в начале координат
OXYZ.
В предположении Al рассматриваемая задача имеет вид
1 п
Е = | UMp^ — Ni min (1)
п z—' Р,и
г=1
при условиях
UUT = I, Р = (р1, ••• ,рп), pi е |1,n], pi = р.,- при i = j (2)
А2. Оси тензоров инерции систем систем М и N совмещены ортогональным преобразованием с осями декартовой системы OXYZ.
Замечание 1. В случае ограничения pi = i, i = 1, • • • ,п, задача (1)-(2) решается алгоритмом Кабша [1,2], результат, действия которого обозначим, АК(М, N). Положение UM подвижной системы, соответствующее результату действия алгоритм,а Кабша, обозначим argminAK(M,N).
Замечание 2. В случае, когда положение подвижной системы в пространстве фиксировано U = I задача (1)-(2) решается венгерским алгоритмом [3], результат, действия которого обозначим AV(M,N). Саму перенумерованную венгерским алгоритмом систему обозначим argminAV(М, N).
Перейдем к определению алгоритмов. В соответствии с [8] определим матрицу U в виде
U(ф,ф, 0) = Щ(ф,ф, в)=Р1А(ф,ф, в), i = 1,2, (3)
где
(cos0 + (1 — cos9)x2 (1 — cos9)xy — zsin0 (1 — cos 9)xz + y sin0 \ (1 — cosd)yx + zsin0 cos0+(1 — cos0)y2 (1 — cosd)yz — x sin0 I , (4)
(1 — cosd)zx — у sine (1 — cos0)zy + x sin0 cos0 + (1 — cos0) z2 J
x = cosф cos ф, у = sinф cos ф, z = sin ф, w = (x,y, z) — направляющий вектор осп поворота. Здесь ф е |—^/2,^/2], ф е |—^/2,^/2] — углы, определяющие направляющий вектор w прямой I, проходящей через начало координат, в е [0, 2^] — угол поворота точки системы М вокруг прямой I против часовой стрелки, Р1 = I, Р2 — матрица отражения [8] относительно прямой I.
Зафиксируем натуральное т. Пусть h = -к/т. Алгоритм 1 (примитивный) Al.l. Е := AV(M,N); Al.2. Для i := 1 • • • 2 выполнить Для j := 0 • • • т выполнить Для к := 0 • • •т выполнить Для I := 0 • • • т выполнить начало
ф := — § + jh, ф := — f + kh, в := 21 h; Е1 := AV (Щ(ф,ф, 6)M,N); если (Е1 < Е)то (Е1 := Е); конец.
Алгоритм 2 (комбинированный)
А1.1. Е := АК(агдтгпАУ(М, М),М); А1.2. Для г := 1 ■ ■ ■ 2 выполнить Для ] := 0 ■ ■ ■ т выполнить Для к := 0 ■ ■ ■т выполнить Для I := 0 ■ ■ ■ т выполнить начало
ф := —§ + зЬ, ф := —2 + кЬ, в := 21Н; Е1 := АК(агдтгпАУ(Щ(ф, ф, в)М, N); если (Е1 < Е) то (Е1 := Е) конец.
Третий алгоритм представляет собой комбинацию алгоритма 2 и итерационного метода покоординатного спуска. Определим сначала алгоритм покоординатного спуска. Обозначим
Е1(ф,ф, г) = шт АУ(и1(ф,ф, в)М,М),Е2(ф,в, г) = шт АУ(и1(ф,ф, в)М,М), Е3(ф,в, г) = шт АУШф,ф, в)М,М).
Фе[-ъ/2,ъ/2]
Через в = агдттЕ\(ф, ф, г), ф = агдтгпЕ2(ф, О, г), ф = агдтгпЕз(ф, в, 1,), обозначим соответствующие значения аргументов, при которых достигается минимум. Зафиксируем ее (0,1], ф € [—л/2,-л/2], ф € [—-к/2, л/2]. Рассмотрим следующий алгоритм, начало в := 0;
Р := АУ(и1(ф,ф, в)М, N) := 1 ■ ■ ■ 2 начало
Р^) := Е\(ф, ф, г); в1 := агдттЕ^ф, ф, г); Р2(г) := Е2(ф, 6и г); ф1 := агдттЕ2(ф, 0Ь г);
Рз(г) := Ез(ф1, 01, г); 61 := агдттЕз(ф, 61, г);
конец
Р4 := шт Р3(г)\
4 ie[l,2]
если |Р4 — Р| > е то (ф := фь ф := фи @ := Р := Р4; идти на 1);
Р := Р4;
конец.
Результат действия данного алгоритма (число Р ) обозначим Р := БСМ(ф,ф, е). Алгоритм 3 (итерационный) А3.1. Е := АУ(М,М) АЗ.2. Для ] := 0 ■ ■ ■ т выполнить к := 0 ■ ■ ■ т начало
ф := — § + ф := — 2 +кН; Е1 := БСМ(ф, ф, е); если (Е1 < Е) то (Е1 := Е); конец.
2 Оценка скорости сходимости алгоритмов
В данном разделе мы докажем следующее утверждение.
Теорема. Для отыскания значения Е с точностью е алгоритмами 1-8 достаточно совершить 0(п3
/е3/2)
арифметических операций.
Е
ше, чем для алгоритма 1, поэтому и для них теорема будет справедлива. Обозначим АУ(и^(ф, ф, в)М, N) = /¿(ф,ф, в), г = 1, 2. Пусть минимум Е1 функции ^(ф,ф, в) достигается при некотором паросочетанни вершин и значениях параметров ф = ф,ф = ф, в = в, г = г .Обозначим через /(ф, ф, 0) функцию, значение
которой равно AV(Щ(ф,ф, 6)M,N) для этого фиксированного паросочетання Р при всех ф, ф, 6. Тогда, очевидно, min/(ф,ф, 6) = min ^(ф,ф, 6) Но в силу своего определения (см. (3)-(4)) f есть дважды диффе-
г=1,2
ренцируемая функция с ограниченными вторыми частными производными. Поэтому в точке минимума ее первый дифференциал df (ф, ф, 6) = 0. Тогда согласно формуле Тейлора
1 3 3 с)2 f
f (ф, ф, 6) = Е + - £ ]Г ^-С-(ф, ф, 6)(xH — xi)(xj -щX (5)
i=i j=i
где обозначено x1 = ф, x2 = ф, x3 = 6, ф = ф + Х(ф — ф), ф = ф + Х(ф — ф), 6 = 6 + Х(6 — 6) для некоторого Л € [0,1]. При этом
d2f (ф,ф, 6) = ± СХС-(ф,ф, в), (в)
dxldxJ ' ' dxldxJ
k=i
где fk(ф, ф, 6) = ± || Щ(ф, ф, 6)Мрк — Nk ||2. Поскольку | (ф, ф, 0)1 < f то из (6) следует, что
Поскольку в 0(1/т) — ^^^^^^тости точки (ф, ф, 6) всегда найдется набор узлов (фj, фк, 6i), по которым идет перебор в алгоритме 1, то из (5)-(7) имеем minf^j,фк, 6i) — Е < Учитывая, что для любых (ф,ф, 6)
будет min fi(ф, ф, 6) < f(ф, ф, 6), получаем, что
г =1,2
с
min fi^j ,фк, 6t ) - Е . (8)
i=i,2 m2
Поскольку каждое значение fi(ф,ф, 6) находится за 0(п3) арифметических операций, а нахождение минимума с точностью е требует в силу (8) значения m = 0(е-1/2), то общее число операций составит 0(п3 /е 3/2). Теорема доказана.
3 Результаты численных экспериментов
Рассматривались три семейства Si,S2атомных решеток, интерпретируемых как совокупности точек в R3. Первое семейство Si состоит из 19 структур, к^дад из которых содержит 9 точек, второе S2 из 12 структур, содержащих 18 точек каждая, третье семейство S3 из 7 структур, содержащих 57 точек каждая. Вычислительный эксперимент для каждой из структур состоял в последовательном отыскании величины (1) для всевозможных пар структур семейств. Приближенное значение этой величины, найденное одним из алгоритмов 1, 2, 3, обозначим Еаррг. В алгоритме 3 минимум на каждой итерации покоординатного спуска определялся методом золотого сечения.
Результаты счета представлены в трех таблицах. В каждой из таблиц для соответствующего семейства представлена зависимость абсолютной погрешности для алгоритмов 1, 2, 3 соответственно, а также время счета на компьютере INTEL(R) Соге(ТМ) 15-2500 CPU @ 3,30GHz 3,60 GHz. При этом для
АЕ = max вычислялось с помощью алгоритма Кабша и полного перебора 9! возможных паро-
Si \ Eappr — Eexact |
сочетаний, а для S2 и S3 в качестве Eexact бралось значение алгоритма 2 для значения т = 128.
Заключение
1. Из результатов экспериментов видно, что алгоритм 2 имеет существенные преимущества перед алгоритмами 1 и 3. При сопоставимом с алгоритмом 1 времени счета, алгоритм 2, начиная с некоторого значения т, фактически, дает точное значение Е, достигая на всех парах оптимального паросочетання перед выполнения алгоритма Кабша. Остаточная погрешность связана с погрешностями численных методов алгебры, используемых в алгоритме Кабша.
Таблица 1: Семейство Si
m 4 8 16 32 64
Алгоритм 1 3.15 • 10-i 1.07 • 10-i 5.59 • 10-2 2.79 • 10-2 6.07 • 10-3
2 сек. 12 сек. 1 мин 20 сек. 9 мин. 39 сек. 1 ч. 12 мин. 42 сек.
Алгоритм 2 5.04 • 10-2 1.45 • 10-5 1.45 • 10-5 1.34 • 10-5 1.27 • 10-5
5 сек. 25 сек. 2мин. 54 сек. 21 мин. 30 сек. 2ч. 44 мин. 13 сек.
Алгоритм 3 7.18 • 10-2 5.86 • 10-2 5.86 • 10-2 5.64 • 10-2 5.64 • 10-2
1 мин. 45 сек. 5 мин. 20 сек. 18 мин. 37 сек. 1ч. 9 мин. 41 сек. 4 ч. 36 мин. 2 сек.
Таблица 2: Семейство S2
т 4 8 16 32 64
Алгоритм 1 Алгоритм 2 Алгоритм 3 3.56 • 10-i 2 сек. 8.26 • 10-i 3 сек. 8.55 • 10-2 1 мин. 10 сек. 2.40 • 10-i 14 сек. 3.11 • 10-2 21 сек. 8.55 • 10-2 3 мин. 35 сек. 1.86 • 10-i 1 мин 38 сек. 3.71 • 10-6 2мин. 20 сек. 8.55 • 10-2 12 мин. 27 сек. 2.45 • 10-2 11 мин. 45 сек. 3.38 • 10-6 17 мин. 43 сек. 8.55 • 10-2 46 мин. 4 сек. 1.38 • 10-2 1 ч. 29 мин. 22 сек. 4.22 • 10-7 2ч. 14 мин. 49 сек. 8.55 • 10-2 2 ч. 57 мин. 38 сек.
Таблица 3: Семейство S3
т 4 8 16 32 64
Алгоритм 1 7.20 • 10-i 5.02 • 10-i 1.60 • 10-i 1.48 • 10-i 1.83 • 10-2
5 сек. 1 мин. 31 сек. 3 мин 35 сек. 26 мин. 40 сек. 3 ч. 33 мин. 51 сек.
Алгоритм 2 3.80 • 10-i 1.53 • 10-i 8.10 • 10-2 1.74 • 10-5 1.39 • 10-6
6 сек. 35 сек. 4мин. 2 сек. 30 мин. 8 сек. Зч. 46 мин. 49 сек.
Алгоритм 3 1.28 • 10-i 1.28 • 10-i 1.15 • 10-2 1.15 • 10-2 1.15 • 10-2
2 мин. 28 сек. 7 мин. 28 сек. 25 мин. 25 сек. 1ч. 24 мин. 40 сек. 6 ч. 1 мин. 19 сек.
2. Алгоритм 3 демонстрирует отсутствие сходимости. Это связано с тем, что метод золотого сечения находит минимум только для унимодальных функций, а условие унимодальности не всегда выполняется. Однако применение более надежных методов отыскания одномерного минимума приведет к существенному возрастанию времени счета. Применение метода золотого сечения в итерационных алгоритмах показало хорошие результаты для отыскания максимума объема пересечения выпуклых многогранников в работе [9]. Однако для задач из настоящей статьи он оказался менее эффективен.
Список литературы
fl] Kabsch W. A solution of the best rotation to relate two sets of vectors // Acta Crystallographica. 1976. V. 32, P. 922-923.
[2] Kabsch W. A discussion of the solution for the best rotation to relate two sets of vectors // Acta Crystallographica. 1978. V. 34, P. 827-828.
[3] H.W. Kuhn. The Hungarian Method for the Assignment Problem // Naval Research Logistics. 2005. V. 52(1), P. 7-21.
[4] Ali Sadeghi, S. Alireza Ghasemi, Bastian Schaefer, Stephan Mohr, Markus A. Lill, Stefan Goedecker. Metrics for measuring distances in configuration spaces // The Journal of chemical physics. 2013. V. 139, 184118
[5] Allen, W. J.; Rizzo, R. C. Implementation Of The Hungarian Algorithm To Account For Ligand Symmetry And Similarity In Structure-Based Design //J. Chem. Inf. Model. 2014. V. 54, P. 518-529.
[6] Matthew Griffiths, Samuel P. Niblett, and David J. Wales. Optimal Alignment of Structures for Finite and Periodic Systems
J. Chem. Theory Comput. 2017. V. 13, P. 4914-4931
[7] Berhane Temelso, Joel М. ЛI л bey. Toshiro Kubota, Nana Appiah-padi, George C. Shields. ArbAlign: A Tool for Optimal Alignment of Arbitrarily Ordered Isomers Using the Kuhn-Munkres Algorithm
J. Chem. Info. Model. 2017. V. 57 (5), P. 1045-1054
[8] Szymanski, John E. (1989). Basic Mathematics for Electronic Engineers:Models and Applications. Taylor & Francis. P. 154. ISBN 0278000681.
[9] Shevchenko A.H., Blatov I. A., Kitaeva E.V., Blatov V.A. Local coordination versus overall topology in crystal structures: deriving knowledge from crystallographic databases
Cristal Growth and Design. 2017. V. 17, №2. P. 774-785
Блатов Игорь Анатольевич — д.ф.-м.н., зав. кафедрой высшей математики Поволжского государственного университета телекоммуникаций и информатики;
e-mail: blatow@mail.ru; Китаева Елена Викторовна — к.ф.-м.н., доцент кафедры дифференциальных уравнений и теории управления; Самарский национальный исследовательский университет имени академика С.П. Королева;
e-mail: el kitaeva@mail.ru.
