CC BY
8Работа первого автора выполнена при частичной финансовой поддержке гранта РФФИ Л" 15 1)1 09242 А, работа второго автора — гранта РФФИ № 15-01-01132. Работа третьего автора выполнена при частичной финансовой поддержке Фонда Саймонса.
СПИСОК ЛИТЕРАТУРЫ
1. Deligne P. La conjecture de Weil. I // Publications Mathématiques de l'IHES. 1974. 43. 273-307.
2. van den Bogaart T., Edixhoven B. Algebraic stacks whose number of points over finite fields is a polynomial // Number fields and function fields — two parallel worlds. Progr. Math. Vol. 239. Boston: Birkhauser, 2005. 39-49.
3. Deligne P., Mumford D. The irreducibility of the space of curves of given genus // Publications Mathématiques de PIHÉS. 1969. 36. 75-109.
4. Etingof P., Henriques A., Kamnitzer J., Rains E. The cohomology ring of the real locus of the moduli space of stable curves of genus 0 with marked points // arXiv:math. AT/0507514, 2005.
5. Getzler E. Operands and moduli spaces of genus 0 Riemann surfaces // The Moduli Space of Curves. Progr. Math. Vol. 129. Boston: Birkhauser, 1995. 199-230. _
6. Bergstrom J., Tommasi O. The rational cohomology of ~Ml // Math. Ann. 2007. 338, N 1. 207-239.
7. van der Geer G. Counting curves over finite fields // Finite Fields and Appl. 2015. 32. 207-232.
8. Keel S. Intersection theory of moduli space of stable n-pointed curves of genus zero // Trans. Amer. Math. Soc. 1992. 330, N 2. 545-574.
Поступила в редакцию 26.10.2016
УДК 519.6
МЕТОД СОГЛАСОВАНИЯ ИЗМЕРЕНИЙ ОПОРНОЙ ФУНКЦИИ ВЫПУКЛОГО ТЕЛА В МЕТРИКЕ L^
И. А. Палачёв1
Предложен новый алгоритм решения задачи согласования измерений опорной функции выпуклого тела в метрике L00, позволяющий получить решение за квадратичное время по числу измерений без использования линейного программирования. Также доказана оценка скорости сходимости, которая имеет место при довольно слабых условиях на входные данные, что дает возможность применить метод к более широкому классу задач, чем это было прежде. Разработанный алгоритм обладает большей и гарантированной стабильностью и предсказуемостью, нежели прочие алгоритмы, существовавшие и применявшиеся для согласования измерений опорной функции. Приведены детали реализации алгоритма и результаты его тестирования.
Ключевые слова: опорная функция, восстановление геометрических тел, теневой контур, преобразование двойственности.
A new algorithm is proposed for estimation of convex body support function measurements in Loo metric, which allows us to obtain the solution in quadratic time (with respect to the number of measurements) not using linear programming. The rate of convergence is proved to be stable for quite weak conditions on input data. This fact makes the algorithm robust for a wider class of problems than it was previously. The implemented algorithm is stable and predictable unlike other existing support function estimation algorithms. Implementation details and testing results are presented.
Key words: support function, geometric body reconstruction, shadow contour, duality-transformation.
1. Введение. Настоящее исследование возникло из практической задачи восстановления трехмерного выпуклого тела по измерениям его опорной функции. Данная задача рассматривалась в
1 Палачёв Илья Александрович — соискатель каф. вычислительной математики мех.-мат. ф-та МГУ, e-mail: palachev .ilyaQyandex .ru.
связи с разнообразными приложениями: в компьютерной томографии [1], при обработке данных лазерного радара [2], в магнитно-резонансной визуализации [3, 4], при построении трехмерных моделей драгоценных камней в ювелирной промышленности [5]. Все эти прикладные задачи так или иначе сводятся к решению так называемой задачи согласования измерений опорной функции выпуклого тела (разъяснение этих понятий см. ниже).
Во всех работах, посвященных решению задачи согласования, проблема сводилась к решению задачи квадратичного (для метрики £2) или линейного (для метрик Ь\ и Ь^) программирования следующего вида:
\\к-к°\\х -тип, С/г ^ 0,
где Н° — исходные значения измерений опорной функции, Л, — неизвестные (переменные задачи), которые требуется найти, а С — линейная матрица.
Основным камнем преткновения при решении данной задачи для измерений опорных функций тел в М3 являлся тот факт, что матрица С содержит в себе квадратичное (0(Ж2), где N — размерность вектора Ь) число строк.
Поэтому главным направлением исследования данной задачи в ряде работ было упрощение условий согласованности, что позволяло записать задачу в более компактном виде и, следовательно, решить ее за гораздо меньшее время. Так, в работе [6] была найдена возможность записи задачи программирования в более удобных терминах точек касания (до того времени само нахождение условий было нетривиальной задачей).
Похожую цель преследовала и статья [5], основной результат которой дал возможность на большом классе входных данных избавиться примерно от 80% условий (последняя оценка, впрочем, была выведена экспериментально). Тем не менее это позволило лишь сдвинуть границу применимости метода, но так и не перевело его на уровень, достаточный, чтобы использовать его при любых входных данных.
Следующей сложностью на пути решения задачи стала довольно типичная для вычислительных методов ситуация: не известен алгоритм, который бы гарантированно работал для любых входных данных. Так, в мире существуют как коммерческие (лишь некоторые из них доступны по академической лицензии), так и свободные программные пакеты для решения задач квадратичного или линейного программирования. Коммерческие пакеты, часто не раскрывающие особенности своей реализации (или оглашающие их лишь частично), не всегда могут удовлетворить пользователя, желающего иметь гарантированное поведение программного пакета.
Свободные программные пакеты, такие, как, например, 1рс^ [7], который, как показано в статье [5], довольно эффективно вычисляет решение данной задачи, имеют открытый исходный код и основаны на результатах научных статей [8, 9], доказывающих свойства реализованных в них алгоритмов. Тем не менее метод внутренней точки, являющийся основой реализации этого пакета, при всех его достоинствах имеет такой недостаток, что может впасть в режим восстановления, при котором существующие оценки скорости сходимости довольно скудны [9]. Иначе говоря, алгоритм при наличии ряда неблагоприятных обстоятельств может работать сколь угодно долго (если не ограничить число итераций) или может вообще не сойтись.
В связи с этим возникает вопрос: можно ли переформулировать задачу согласования опорных чисел таким образом, чтобы алгоритм решения имел гарантированную приемлемую скорость сходимости при любых входных данных?
Настоящая работа призвана представить такой алгоритм для задачи в метрике Ь^. Метод позволяет свести задачу к алгоритму двоичного поиска на отрезке действительных чисел, на каждом шаге которого применяются алгоритмы, имеющие гарантированную вычислительную сложность (построение выпуклой оболочки, задача обнаружений пересечения отрезка и многогранника, линейные преобразования). Более того, как будет показано далее, алгоритм при любых входных данных имеет гарантированную (квадратичную) скорость сходимости, какую метод внутренней точки позволял иметь лишь при выполнении большого ряда ограничений на входные данные.
Также представлены результаты реализации и тестирования алгоритма, подтверждающие теоретически доказанные оценки и демонстирующие приемлемую производительность на ряде примеров, полученных из эмпирических данных для восстановления тела по его теневым контурам.
2. Основные понятия. Приведем далее формализацию задачи согласования измерений опорной функции выпуклого тела.
Определение 1. Опорной функцией выпуклого тела К С Кга называется функция Нк '■ —>■
М+:
Ьк(и) = 8ир(ж, и),
хек
где и — произвольный вектор в Мга, (ж, и) — евклидово скалярное произведение этого вектора на некоторую точку х € К.
В перечисленных выше работах рассматривается следующая задача. Пусть имеется некоторое реальное физическое тело К и производятся эмпирические измерения его опорной функции Л-1, Л-2, ■ ■ ■, Нт (опорные числа) на конечном наборе единичных векторов и\, 11,2, ■ ■ ■, ит (опорные направления):
Ы = Нк{щ)+£1, г = 1,2,... ,т,
где £г — погрешности измерения опорных чисел
Определение 2. Опорные числа ¡1\, ¡12, ■ ■ ■, Ьт по направлениям щ,и2, ■ ■ ■, ит называются согласованными, если существует такое выпуклое тело К* С Мга, что Нк*{щ) = Ы-
Поскольку не всякий набор действительных чисел является согласованным набором опорных чисел некоторого выпуклого тела, естественным образом возникает задача согласования опорных чисел: по набору опорных чисел Н®, Щ,..., и опорным направлениям щ, ь,2, ■ ■ ■, ит строится такой согласованный набор опорных чисел Н\, Щ,..., ¡1^, который реализует минимум функционала
1=\\Н-Н°\\Х,
где X — одна из метрик в Мт (например, ¿2, Ь\, Ь^).
В работе [6] задача согласования опорных чисел сводится к следующей задаче математического программирования:
\\{{Хг,иг)}1^1 - Ь°\\х т^ (Хг,щ) ^ (Ху,щ), 1 ^ % ф ] ^ ГП,
где х% € К3 — точки касания тела с опорной плоскостью, имеющей нормаль щ.
3. Быстрый алгоритм согласования. Основным результатом работы [5] является следующая Теорема 1. Рассмотрим следующую систему ограничений относительно переменных точек Хг € М3, г = 1,..., п, и переменного вещественного числа е > 0;
(Хг,щ) ^ (Хг,щ), | (ж*, щ) ~ Ы | ^ £, (1)
где щ, г = 1,..., п, — набор фиксированных опорных направлений, нормированных в евклидовой метрике (т.е. \\щ\\2 = = 1 ,...,п); Нг,г = 1 ,...,п, — набор соответствующих им фиксированных опорных чисел.
Пусть система ограничений (1) разрешима для некоторого е°. Введем, многогранник К по следующей формуле:
т
К = (~){(х,щ) ^ы+е0}. (2)
г=1
Пусть 5 — полярное преобразование двойственности. Тогда если отрезок, соединяющий точки 5({(х,и¿) = Л-г — е0}) и 5({(х,ь^) = Ъз — е°}); пересекает, 5(К), то условия (х^щ) ^ и,
(ху,Пу) ^ (Жг,«^) избыточны в системе ограничений (1).
Рассмотрим задачу согласования опорных чисел в метрике Ь^:
е ->• (Хг,щ) ^ (Хг,щ), \(Хг,щ) - Ы\ ^ е. (3)
По сути данная задача линейного программирования состоит в нахождении такого минимального е, при котором система ограничений (1) является разрешимой.
Теорема 2. Система ограничений (1) разрешима тогда и только тогда, когда ни одна из точек 5({(х,щ) = ^ — е}) не лежит во внутренности многогранника 5(К), где К — многогранник, определяемый соотношением (2).
Доказательство. Точка 5({(х,щ) = /г» — е}) не лежит во внутренности многогранника 5(К) тогда и только тогда, когда плоскость {(х,щ) = /г» — е} пересекает многогранник К (по одному из свойств преобразования двойственности). В данном случае можно положить Хг = Ьк{щ), и тогда
(Хг, щ) ^Ы-е, {Хг,щ) ^ Щ + £, ] = 1, . . . , 171.
Построим таким образом набор точек (х\,... ,хт). По определению опорной функции условия (хг,щ) ^ (хг,и^ выполняются автоматически. □
Полученный результат позволяет организовать процесс нахождения решения задачи линейного программирования (3) с помощью следующего алгоритма.
1. По исходным опорным данным щ,Нг,г = 1,... ,т, построить тело К0 как пересечение полупространств:
т
к° = (~){(иг,х) ^ ы}.
г=1
2. Пусть А\,..., А3 — вершины тела К0. Вычислить опорные числа тела К0 по направлениям щ, г = 1,..., т, следующим образом:
= тах (щ,А
3. Вычислить ¿оо-расстояние от полученных опорных чисел до исходных:
е° = тах —
г=1 ,...,т
Положить = 0, = ео.
4. Если |в_|_ — < Ае, то остановить алгоритм, в противном случае положить е = и перейти к следующему шагу.
5. Вычислить точки ¿({(ж,щ) = Нг + е}), г = 1,... ,гп, и построить их динамическую выпуклую оболочку И.
6. Вычислить точки ¿({(ж,щ) = Нг — е}), г = 1,... ,гп, и проверить, лежит ли хотя бы одна из них во внутренности И. Если хотя бы одна лежит, то изменить величину = е, иначе изменить величину = е. Перейти к шагу 4.
На шаге 4 величина Ае фиксирована, т.е. является параметром алгоритма, определяющим точность вычисления минимального допустимого е. Количество итераций, необходимое для завершения алгоритма, очевидно, равно д^. На каждом шаге алгоритма требуется 2т раз вычислять двойственные образы плоскостей (за константное число операций), один раз строить динамическую выпуклую оболочку т точек (что требует 0(т2) времени при использовании алгоритма построения трехмерной триангуляции Делоне [10]) и т раз проверять принадлежность точки построенному многограннику (что требует 0(1(щт) времени для каждой из т точек по отдельности [10]). Следовательно, трудоемкость работы алгоритма составляет
Поскольку в реальных задачах величины ео и Ае различаются не больше чем на несколько порядков, то можно считать, что представленный алгоритм имеет квадратичную сложность.
Также представляется целесообразным проверить гипотезу, что при соблюдении разумного предположения о равномерной распределенности опорных направлений щ,г = 1 ,...,п, на сфере возможно построение триангуляции Делоне за линейное время [11], что позволило бы говорить о линейно-логарифмической сложности алгоритма.
4. Реализация и тестирование алгоритма. Алгоритм был реализован на основе библиотеки алгоритмов вычислительной геометрии СОАЬ1. Пересечение полупространств на шаге 1 алгоритма и выпуклая оболочка на шаге 5 вычислялись на основе алгоритма построения трехмерной триангуляции Делоне [10]. Выбор в пользу данного алгоритма был сделан, поскольку он позволяет производить быструю многократную проверку принадлежности точки выпуклой оболочке [10].
1 Computational Geometry Algorithms Library, http://www.cgal.org.
Тестирование было произведено на опорных числах, полученных эмпирическим путем на основе измерений реального выпуклого тела (драгоценного камня). В метрике L^ алгоритм работает приблизительно в 300 раз быстрее, чем существующий алгоритм решения задачи линейного программирования, основанный на методе внутренней точки [7]. К примеру, решение задачи на 2912 опорных числах было получено за 0,98 с, в то время как предыдущий алгоритм позволял получить решение лишь за 300 с [5] (на одних и тех же вычислительных ресурсах).
Безусловно, основываясь лишь на тестировании алгоритма, нельзя судить о превосходстве или недостатке предлагаемого метода по отношению к предыдущему. Тем не менее определенные выводы возможно сделать на основе теоретической оценки скорости работы алгоритма, представленного в работе [5]. По теореме 4.7 статьи [8] фильтрационный метод внутренней точки с линейным поиском при выполнении ряда условий сходится сверхлинейно к минимуму задачи (3). Поскольку задача (3) содержит квадратичное число ограничений (см. [6, 5]), на каждом шаге метода выполняются вычисления общей трудоемкостью не менее 0(т2), т.е. такой же, какую дает предлагаемый метод. Между тем на практике оказалось, что по меньшей мере одно условие теоремы 4.7 из [8] не выполняется в практических примерах задачи (3), а именно алгоритм не должен входить в фазу восстановления (когда производится минимизация суммарного отклонения ограничений). На примере, рассмотренном в [5], алгоритм внутренней точки входит в фазу восстановления. Следовательно, в данном случае мы имеем лишь сходимость, гарантируемую теоремой 2 статьи [9], в которой не приводятся оценки скорости сходимости.
Таким образом, предлагаемый метод позволяет за квадратичное время гарантированно получить решение задачи (3), в отличие от метода внутренней точки, для скорости сходимости которого применительно к задаче (3) содержательных оценок пока нет.
СПИСОК ЛИТЕРАТУРЫ
1. Prince J.L., Willsky A.S. Reconstructing convex sets from support line measurements // IEEE Trans. Pattern Anal, and Machine Intel. 1990. 12, N 4. 377-389.
2. Lele A.S., Kulkarni S.R., Willsky A.S. Convex-polygon estimation from support-line measurements and applications to target reconstruction from laser-radar data //J. Opt. Soc. Amer. A. 1992. 9, N 10. 1693-1714.
3. Gregor J., Rannou F.R. Least-squares framework for projection MR I reconstruction // Proc. SPIE. 2001. 4322. 888-898.
4. Gregor J., Rannou F.R. Three-dimensional support function estimation and application for projection magnetic resonance imaging // Int. J. Imaging Systems and Technol. 2002. 12, N 1. 43-50.
5. Палачёв И.А. Метод исключения избыточных ограничений в задаче восстановления тела по измерениям его опорной функции // Вычисл. методы и програм. 2015. 16. 348-359.
6. Gardner R.J., Kiderlen М. A new algorithm for 3D reconstruction from support functions // IEEE Trans. Pattern Anal, and Machine Intel. 2009. 31, N 3. 556-562.
7. Wachter A., Biegler L.T. On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming // Math. Program. 2006. 106, N 1. 25-57.
8. Wachter A., Biegler L.T. Line search filter methods for nonlinear programming: Local convergence // SIAM J. Optim. 2005. 16, N 1. 32-48.
9. Wachter A., Biegler L. T. Line search filter methods for nonlinear programming: Motivation and global convergence // SIAM J. Optim. 2005. 16, N 1. 1-31.
10. Devillers O. The Delaunay hierarchy // Int. J. Found. Comput. Sci. 2002. 13. 163-180.
11. Attali D., Boissonnat J.-D. A linear bound on the complexity of the Delaunay triangulation of points on polyhedral surfaces // Discr. and Сотр. Geom. 2004. 31, N 3. 369-384.
Поступила в редакцию 19.09.2016
