Викиконспекты - Вклад участника [ru]

Алгоритмы точного вычисления гиперобъема

2012-06-20T12:23:00Z

Joshik:

{{В разработке}}

== Постановка задачи ==
<tex>x = (x_1, x_2, ..., x_d; x_i \ge 0) \in \mathbb{R}^d</tex> - точка в <tex>d</tex>-мерном пространстве.

Точка <tex>x</tex> доминирует точку <tex>y</tex> (<tex>x \succ y</tex>), если <tex>\forall i : x_i \ge y_i, \exists j : x_j > y_j</tex>.

<tex>X = \{x^1, x^2, ..., x^n\} \subset \mathbb{R}^d</tex> - множество из <tex>n</tex> точек в <tex>d</tex>-мерном пространстве таких, что <tex>\nexists i \neq j : x_i \succ x_j</tex> - никакая точка не доминируется другой точкой из этого множества.

<tex>S(X) = \mu (\bigcup \limits_{x \in X} \{y | x \succ y\}) </tex> - гиперобъем множества <tex>X</tex>, где <tex>\mu</tex> - мера Лебега.

В частности, если <tex>X = \{x\}</tex>, то <tex>S(X) = \prod \limits_{i=1}^{d} x_i</tex>.

Задача: найти точное значение гиперобъема <tex>S(X)</tex> множества из <tex>n</tex> точек <tex>d</tex>-мерного пространства.

== Алгоритм включения-исключения (Inclusion-Exclusion Algorithm, IEA) ==

Самый простой алгоритм нахождения гиперобъема базируется на идее комбинаторной [[формула включения-исключения|формулы включения-искючения]].
Все множество <tex>X</tex> представляется в виде объединения <tex>n</tex> гиперпараллелепипедов (<tex>X^i</tex>), соответствующих отдельным точкам <tex>x^i</tex>.

После этого объем всего множества вычисляется по формуле: <center> <tex> S(X) = \sum \limits_{I \in 2^n} (-1)^{|I|+1} S(\bigcap \limits_{ j \in I} X^j) </tex> </center>

Объем пересечения гиперпараллелепипедов легко считается как произведение по каждой координате минимального значения этой координаты среди всех точек, которым соответствуют гиперпараллелепипеды.

Таким образом, в этом алгоритме перебираются все подмножества точек множества <tex>X</tex>, для каждого множества находится гиперобъем пересечения соответствующих гиперпараллелепипедов, и он прибавляется с соответствующим знаком к ответу.
Время работы этого алгоритма составляет <tex>O(n 2^n)</tex>.

== Алгоритм LebMeasure ==

Алгоритм LebMeasure<ref>Bringmann K., Friedrich. T.: An Efficient Algorithm for Computing Hypervolume Contributions. Evolutionary Computation, Vol. 18, No. 3, MIT Press, pp 383-402. (2010)</ref> обрабатывает точки множества <tex>X</tex> по очереди. Для каждой очередной точки <tex>x</tex> находится объем некоторого максимального по включению гиперпараллелепипеда, доминируемого эксклюзивно только этой точкой <tex>x</tex> и она заменяется на некоторое множество ''порожденных'' точек, которые доминируют оставшуюся область, доминировавшуюся до этого точкой <tex>x</tex>.

Например, если изначально было четыре точки в трехмерном пространстве <tex>(6, 7, 4), (9, 5, 5), (1, 9, 3), (4, 1, 9)</tex>, то точка <tex>(6, 7, 4)</tex> эксклюзивно доминирует параллелепипед с одним концов в <tex>(6, 7, 4)</tex>, а другим - в <tex>(4, 5, 3)</tex>. После добавления объема этого параллелепипеда к ответу, точка <tex>(6, 7, 4)</tex> порождает три точки: <tex>(4, 7, 4), (6, 5, 3), (6, 7, 3)</tex>. При этом точка <tex>(6, 5, 4)</tex> доминируется точкой <tex>(9, 5, 5)</tex> и сразу удаляется из множества <tex>X</tex>.

Обработка точек продолжается, пока все точки не будут обработаны. Таким образом, время работы алгоритма напрямую зависит от количества порожденных точек. Легко заметить, что таких точек не более, чем <tex>n^d</tex>, поскольку каждая координата каждой порожденной точки равна соответствующей координате некоторой точки исходного множества <tex>X</tex>.

К сожалению, эта верхняя оценка является достижимой. Например, если исходное множество <tex>X</tex> имеет вид: <tex>(1, n, n, ..., n), (2, n - 1, n - 1, ..., n - 1), ..., (n, 1, 1, ..., 1)</tex>, и точки обрабатываются именно в этом порядке, то всего будет обработано <tex>n^{d-1}</tex> точка, что показано в <ref name = "while"/>. Правда, если в этом примере точки обрабатываются в обратном порядке, то суммарное количество обработанных точек линейно зависит от <tex>n</tex> и <tex>d</tex>. Тем не менее, существуют примеры, для которых любой порядок обработки приводит к экспоненциальной зависимости числа порожденных точек от размерности пространства <tex>d</tex> и эта зависимость близка к <tex>n^d</tex> <ref name="while">
While, L., Hingston, P., Barone, L., Huband, S.: A Faster Algorithm for Calculating Hypervolume. IEEE Transaction on Evolutionary Computation, Vol.10, No. 1, pp 29–38. (2006)</ref>.

== Алгоритм Hypervolume by Slicing Objectives (HSO) ==

Под ''Objectives'' в названии данного алгоритма имеются в виду координаты пространства <tex>\mathbb{R}^d</tex>.

Если алгоритм LebMeasure по очереди рассматривает все точки, то алгоритм HSO<ref>Zhou X., Sun C., Mao N., Li W.: Generalization of HSO algortihms for computing hypervolume for mutiobjective optimization problems, IEEE Congress on Evolutionary Computation. (2007)</ref> рассматривает по очереди все координаты, сводя задачу к меньшей на единицу размерности.

Изначально все точки сортируются по первой координате. Значения этой координаты используются для ''расслоения'' (разрезания) всего множества на <tex>n</tex> частей, внутри каждой из которых при движении вдоль первой координаты форма ''разреза'' (гиперплоскости), проходящего перпендикулярно оси первой координаты, не меняется. Таким образом, для подсчета объема каждой части необходимо найти объем разреза, и умножить на длину части вдоль первой координаты. При этом получившийся разрез имеет на единицу меньшую размерность.
Заметим, что после сортировки и расслоения первая часть содержит все исходные <tex>n</tex> точек множества <tex>X</tex>, вторая - все, помимо точки с минимальной первой координатой, вдоль которой происходит расслоение, и т.д., а последняя часть содержит только одну точку с максимальной этой координатой.

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

Описанный процесс можно реализовать как в виде рекурсивной процедуры, расслаивающей множество вдоль очередной координаты и вызывающей себя рекурсивно для каждой части и следующей координаты, так и нерекурсивно, если поддерживать множество всех текущих частей и на очередной итерации разбивать их все вдоль очередной координаты.

Время работы алгоритма напрямую зависит от суммарного количества частей, на которые будет разбито исходное множество. По аналогичным предыдущему алгоритму рассуждениям легко показывается, что частей не более <tex>n^d</tex>. Тем не менее, существует более точная оценка.

{{Теорема
|statement=Суммарное количество частей, полученных алгоритмом HSO, не превышает <center><tex>f(n, d) = {n + d - 2 \choose d - 1}</tex>.</center>

|proof=
База:<center><tex>f(n, 1) = 1</tex>.</center>

Также, из алгоритма расслоения следует, что из части с <tex>n</tex> точками будет получено по одной части для каждого <tex>i = 1...n</tex>, то есть: <center><tex>f(n, d) = \sum \limits_{i = 1}^{n}f(i, d - 1)</tex>.</center>

Докажем полезные леммы:

{{Лемма
|statement=<center><tex>{x \choose y} = {x - r \choose y} + \sum \limits_{i = 1}^{r} {x - i \choose y - 1}</tex>.</center>
|proof=Индукция по <tex>r</tex>. База: <center><tex>{x \choose y} = {x - 1 \choose y} + {x - 1 \choose y - 1}</tex>.</center>
Переход: предполагаем, что утверждение верно для <tex>r = j</tex>. Раскрываем первое слагаемое в правой части предположения аналогично базе и получаем:
<center><tex>{x \choose y} = {x - j - 1 \choose y} + {x - (j + 1) \choose y - 1} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1} = {x - (j + 1) \choose y} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1}</tex>.</center>
}}

{{Лемма
|statement=<center><tex>{x \choose y} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
|proof=Подставляем в предыдущую лемму <tex>r = x - y</tex> и получаем:
<center><tex>{x \choose y} = {y \choose y} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = {x - (x - y + 1) \choose y - 1} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
}}

продолжаем доказательство теоремы по индукции:
<center><tex>f(n, d) = {n + d - 2 \choose d - 1} = \sum \limits_{i = 1}^{n} {n + d - 2 - i \choose d - 2} = \sum \limits_{k = n}^{1} {n + d + 2 - (n - k + 1) \choose d - 2} = \sum \limits_{k = 1}^{n} {k + d - 3 \choose d - 2} = \sum \limits_{k = 1}^{n}{k + (d - 1) - 2 \choose (d - 1) - 1} = \sum \limits_{k = 1}^{n}{f(k, d - 1)}</tex></center>

}}

Таким образом, на этот алгоритм получена более низкая оценка времени работы, чем на предыдущий, и это подтверждается на практике <ref name ="while"/>.

== Сведение к задаче KMP (Klee's Measure Problem) ==

Задача KMP состоит в нахождении объема объединения прямоугольных гиперпараллелепипедов в <tex>d</tex>-мерном пространстве. Как показано в описании алгоритма IEA, исходная задача легко сводится к этой, если каждой точке поставить в соответствие гиперпараллелепипед с одной вершиной в центре координат, а противоположной - в этой точке.

Существует несколько алгоритмов решения задачи KMP<ref name = "overmars">Overmars M.H., Yap C.K.: New upper bounds in Klee’s measure problem. SIAM Journal on Computing 20(6), pp 1034–1045. (1991)</ref><ref>Beume N.: S-Metric calculation by considering dominated hypervolume as Klee’s measure problem. Evolutionary Computation, 17(4), pp 477–492. (2009)</ref><ref>Beume N., Rudolph G.: Faster S-metric calculation by considering dominated hypervolume as Klee’s measure problem. In Proc. Second International Conference on Computational Intelligence (IASTED ’06), pp 233–238. (2006)</ref>, самый оптимальный из которых использует идеи сканирующей гиперплоскости, [[Статистики на отрезках. Корневая эвристика|корневой эвристики]] и КД-дерево, и позволяет решать задачу за время <tex>O(n^{d/2}\log n)</tex>. Рассмотрим его.

Для начала изложим идею сканирующей гиперплоскости. Рассмотрим всевозможные значения <tex>d</tex>-ой координаты у точек множества <tex>X</tex>.
Возьмем гиперплоскость, перпендикулярную оси <tex>d</tex>-ой координаты, будем ее двигать по этим значениям от минимального до максимального и рассматривать проекцию всех гиперпараллелепипедов на эту гиперплоскость. Заметим, что между точками проекция меняться не будет, а в каждой точке будут появляться или исчезать проекции некоторых гиперпараллелепипедов. Поэтому, если уметь эффективно пересчитывать объем проекции при каждом появлении/исчезновении и прибавлять это значение, умноженное на расстояние между последовательными значениями <tex>d</tex>-ой координаты, то можно легко посчитать суммарный гиперобъем. Как раз для эффективного пересчета используется такая структура, как КД-дерево.

Теперь изложим идею алгоритма на примере трехмерного пространства.
В трехмерном пространстве сканирующая гиперплоскость превращается в обычную плоскость. Будем считать, что она всегда располагается перпендикулярно оси <tex>OZ</tex> и двигается от минимальных значений по этой оси до максимальных. В каждый момент проекция всех параллелепипедов на эту плоскость представляет собой множество прямоугольников и необходимо уметь обрабатывать следующие три запроса:
# добавить прямоугольник
# удалить прямоугольник
# посчитать площадь объединения всех прямоугольников

{{Определение
|definition=КД-деревом называется [[Дерево поиска, наивная реализация|сбалансированное двоичное дерево]], где каждой вершине <tex>\alpha</tex> соответствует некоторая область <tex>C_{\alpha}</tex> <tex>d</tex>-мерного пространства, удовлетворяющая следующим свойствам:
# Корневой вершине <tex>root</tex> соответствует <tex>C_{root}</tex> - все пространство.
# Каждое <tex>C_{\alpha}</tex> является гиперпараллелепипедом (возможно бесконечным в некоторых направлениях).
# Для каждой вершины <tex>\alpha</tex> с детьми <tex>\beta</tex> и <tex>\gamma</tex> верно, что области <tex>C_{\beta}</tex> и <tex>C_{\gamma}</tex> не пересекаются и <tex>C_{\alpha} = C_{\beta} \cup C_{\gamma}</tex>.
}}

Будем хранить прямоугольники в КД-дереве следующим образом:
# Для каждого листа <tex>\alpha</tex> будем хранить все прямоугольники, которые пересекают внутреннюю часть области <tex>C_{\alpha}</tex>.
# Для каждой внутренней вершины <tex>\alpha</tex> будем считать число <tex>Cover_{\alpha}</tex> - количество прямоугольников, которые полностью покрывают область <tex>C_{\alpha}</tex>, но не полностью покрывают область родителя этой вершины <tex>C_{father(\alpha)}</tex>.

Также для каждой вершины будем хранить значение площади пересечения области этой вершины со всеми прямоугольниками <tex>M</tex>, которое определяется следующим образом:
# Для каждого листа <tex>M</tex> равно площади пересечения всех прямоугольников с областью этого листа.
# Для каждой внутренней вершины <tex>\alpha</tex> если <tex>Cover_{\alpha} > 0</tex>, тогда <tex>M</tex> равно площади всей области <tex>C_{\alpha}</tex>, иначе <tex>M</tex> равно сумме этих значений для двух сыновей этой вершины.

Таким образом, значение для корневой вершины <tex>M_{root}</tex> и является искомой площадью объединения всех прямоугольников.

Рассмотрим операцию добавления прямоугольника (box) в КД-дерево.

procedure Insert(box, <tex>\alpha</tex>)
if <tex>\alpha</tex> - лист then
добавить box в эту вершину
пересчитать <tex>M_{\alpha}</tex>
elseif box покрывает область <tex>C_{\alpha}</tex> then
<tex>Cover_{\alpha} = Cover_{\alpha} + 1</tex>
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
elseif box пересекает область <tex>C_{\alpha}</tex> then
Insert(box, leftson(<tex>\alpha</tex>))
Insert(box, rightson(<tex>\alpha</tex>))
if <tex>Cover_{\alpha} > 0 </tex> then
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
else
<tex>M_{\alpha} = M_{leftson(\alpha)} + M_{rightson(\alpha)}</tex>

Аналогично устроена операция удаления прямоугольника - она выполняет обратные действия.

Рассмотрим теперь как именно будет устроено разделение пространства на области КД-деревом для эффективного его использования.
Возьмем проекции всех параллелепипедов на плоскость, перпендикулярную оси <tex>OZ</tex> - множество прямоугольников <tex>V</tex>. Рассмотрим все <tex>x</tex>-координаты границ прямоугольников и разобьем ось <tex>OX</tex> на <tex>2\sqrt{n}</tex> промежутков так, чтобы в каждом из них было не более <tex>\sqrt{n}</tex> границ. Получившиеся <tex>2\sqrt{n}</tex> полос разобьем горизонтальными линиями на итоговые области по отдельности следующим образом. Для каждой полосы у всех прямоугольников, имеющих вертикальную границу в этой полосе, возьмем все горизонтальные границы (их будет не более <tex>2\sqrt{n}</tex>), а среди всех горизонтальных границ, пересекающих область, возьмем каждую <tex>\sqrt{n}</tex>-ую. Таким образом, в итоге получим не более <tex>4\sqrt{n}</tex> областей в каждой вертикальной полосе.

Полученные области обладают следующими свойствами:
# Всего <tex>O(n)</tex> областей.
# Каждый прямоугольник частично покрывает не более <tex>O(\sqrt{n})</tex> областей.
# Никакая вершина прямоугольника не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти утверждения несложно следуют из алгоритма разбиения на области.

Каждая из этих областей будет одним из листов КД-дерева. Для построения всего дерева сначала объединим все области для одной вертикальной полосы, сливая соседние области в одну вершину дерева. После этого аналогично будем сливать соседние вертикальные области и в итоге получим корневую вершину, представляющую всю плоскость.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n)</tex> вершин.
# Каждый прямоугольник хранится в <tex>O(\sqrt{n})</tex> листьях.
# Каждый прямоугольник влияет на <tex>O(\sqrt{n} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина прямоугольника не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти свойства следуют из свойств разбиения на области и логарифмической высоты дерева.

Наконец, получаем итоговый алгоритм. Идем сканирующей плоскостью, перпендикулярной оси <tex>OZ</tex>, добавляя или удаляя прямоугольники из КД-дерева. При этом каждый раз пересчитываем площадь объединения всех текущих прямоугольников и, складывая эти величины, умноженные на соответствующую разницу <tex>z</tex>-координат, получаем общий объем.

Каждая операция добавления или удаления прямоугольников требует изменения <tex>O(\sqrt{n} \log n)</tex> полей <tex>Cover</tex>, а также операций добавления или удаления прямоугольника в <tex>O(\sqrt{n})</tex> листьев, каждая из которых требует <tex>O(\log n)</tex> времени. Итого, на одну операцию добавления или удаления требуется <tex>O(\sqrt{n} \log n)</tex> времени, суммарно таких операций <tex>O(n)</tex>, и получается, что весь алгоритм работает за <tex>O(n \sqrt{n} \log n)</tex>. При этом на хранение всей информации требуется <tex>O(n\sqrt{n}) </tex> памяти - каждый лист может хранить <tex>O(\sqrt{n})</tex> прямоугольников.

Все идеи этого алгоритма обобщаются с трехмерного пространства на <tex>d</tex>-мерное следующим образом.

Построение разбиения на области в <tex>d</tex>-мерном пространстве вместо двухмерного ведется сначала аналогичным образом, только после того, как в двухмерном пространстве получаются итоговые прямоугольные области, здесь продолжаются аналогичные разбиения вдоль третьей координаты - берутся все границы третьей координаты гиперпараллелепипедов, имеющие границу по первой или второй координате строго внутри области (их получается не более <tex>O(\sqrt{n})</tex>), а также каждую <tex>O(\sqrt{n})</tex>-ую границу по третьей координате гиперепараллелепипедов, пересекающих эту область, и так далее.

Это разбиение на области имеет следующие аналогичные свойства:
# Всего <tex>O(n^{d/2})</tex> областей.
# Каждый гиперпараллелепипед частично покрывает не более <tex>O(n^{(d-1)/2})</tex> областей.
# Никакая вершина гиперпараллелепипеда не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Общее дерево строится аналогичным образом и оно в итоге состоит из <tex>d</tex> уровней, каждый из которых имеет высоту <tex>O(\log n)</tex> в дереве - самый верхний уровень содержит разбиение по первой координате, следующий - по второй и так далее.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n^{d/2})</tex> вершин.
# Каждый гиперпараллелепипед хранится в <tex>O(n^{(d-1)/2})</tex> листьях.
# Каждый гиперпараллелепипед влияет на <tex>O(n^{(d-1)/2} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина гиперпараллелепипеда не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Операции в этом дереве выполняются аналогично трехмерному случаю и каждая операция добавления или удаления в КД-дерево размерности <tex>d</tex> занимает <tex>O(n^{(d-1)/2} \log n)</tex> времени. Выполняя сканирование гиперплоскостью вдоль одной из координат получаем КД-дерево размерности <tex>d-1</tex>, в котором необходимо провести <tex>O(n)</tex> операций, то есть суммарное время работы алгоритма равно: <tex>O(n \times n^{(d-1-1)/2} \log n) = O(n^{d/2} \log n)</tex>. При этом хранение КД-дерева требует <tex>O(n^{d/2})</tex> памяти.

== Источники ==
<references/>

Алгоритмы точного вычисления гиперобъема

2012-06-20T12:07:05Z

Joshik:

{{В разработке}}

== Постановка задачи ==
<tex>x = (x_1, x_2, ..., x_d; x_i \ge 0) \in \mathbb{R}^d</tex> - точка в <tex>d</tex>-мерном пространстве.

Точка <tex>x</tex> доминирует точку <tex>y</tex> (<tex>x \succ y</tex>), если <tex>\forall i : x_i \ge y_i, \exists j : x_j > y_j</tex>.

<tex>X = \{x^1, x^2, ..., x^n\} \subset \mathbb{R}^d</tex> - множество из <tex>n</tex> точек в <tex>d</tex>-мерном пространстве таких, что <tex>\nexists i \neq j : x_i \succ x_j</tex> - никакая точка не доминируется другой точкой из этого множества.

<tex>S(X) = \mu (\bigcup \limits_{x \in X} \{y | x \succ y\}) </tex> - гиперобъем множества <tex>X</tex>, где <tex>\mu</tex> - мера Лебега.

В частности, если <tex>X = \{x\}</tex>, то <tex>S(X) = \prod \limits_{i=1}^{d} x_i</tex>.

Задача: найти точное значение гиперобъема <tex>S(X)</tex> множества из <tex>n</tex> точек <tex>d</tex>-мерного пространства.

== Алгоритм включения-исключения (Inclusion-Exclusion Algorithm, IEA) ==

Самый простой алгоритм нахождения гиперобъема базируется на идее комбинаторной [[формула включения-исключения|формулы включения-искючения]].
Все множество <tex>X</tex> представляется в виде объединения <tex>n</tex> гиперкубов (<tex>X^i</tex>), соответствующих отдельным точкам <tex>x^i</tex>.

После этого объем всего множества вычисляется по формуле: <center> <tex> S(X) = \sum \limits_{I \in 2^n} (-1)^{|I|+1} S(\bigcap \limits_{ j \in I} X^j) </tex> </center>

Объем пересечения гиперкубов легко считается как произведение по каждой координате минимального значения этой координаты среди всех точек, которым соответствуют гиперкубы.

Таким образом, в этом алгоритме перебираются все подмножества точек множества <tex>X</tex>, для каждого множества находится гиперобъем пересечения соответствующих гиперкубов, и он прибавляется с соответствующим знаком к ответу.
Время работы этого алгоритма составляет <tex>O(n 2^n)</tex>.

== Алгоритм LebMeasure ==

Алгоритм LebMeasure<ref>Bringmann K., Friedrich. T.: An Efficient Algorithm for Computing Hypervolume Contributions. Evolutionary Computation, Vol. 18, No. 3, MIT Press, pp 383-402. (2010)</ref> обрабатывает точки множества <tex>X</tex> по очереди. Для каждой очередной точки <tex>x</tex> находится объем некоторого максимального по включению гиперкуба, доминируемого эксклюзивно только этой точкой <tex>x</tex> и она заменяется на некоторое множество ''порожденных'' точек, которые доминируют оставшуюся область, доминировавшуюся до этого точкой <tex>x</tex>.

Например, если изначально было четыре точки в трехмерном пространстве <tex>(6, 7, 4), (9, 5, 5), (1, 9, 3), (4, 1, 9)</tex>, то точка <tex>(6, 7, 4)</tex> эксклюзивно доминирует куб с одним концов в <tex>(6, 7, 4)</tex>, а другим - в <tex>(4, 5, 3)</tex>. После добавления объема этого куба к ответу, точка <tex>(6, 7, 4)</tex> порождает три точки: <tex>(4, 7, 4), (6, 5, 3), (6, 7, 3)</tex>. При этом точка <tex>(6, 5, 4)</tex> доминируется точкой <tex>(9, 5, 5)</tex> и сразу удаляется из множества <tex>X</tex>.

Обработка точек продолжается, пока все точки не будут обработаны. Таким образом, время работы алгоритма напрямую зависит от количества порожденных точек. Легко заметить, что таких точек не более, чем <tex>n^d</tex>, поскольку каждая координата каждой порожденной точки равна соответствующей координате некоторой точки исходного множества <tex>X</tex>.

К сожалению, эта верхняя оценка является достижимой. Например, если исходное множество <tex>X</tex> имеет вид: <tex>(1, n, n, ..., n), (2, n - 1, n - 1, ..., n - 1), ..., (n, 1, 1, ..., 1)</tex>, и точки обрабатываются именно в этом порядке, то всего будет обработано <tex>n^{d-1}</tex> точка, что показано в <ref name = "while"/>. Правда, если в этом примере точки обрабатываются в обратном порядке, то суммарное количество обработанных точек линейно зависит от <tex>n</tex> и <tex>d</tex>. Тем не менее, существуют примеры, для которых любой порядок обработки приводит к экспоненциальной зависимости числа порожденных точек от размерности пространства <tex>d</tex> и эта зависимость близка к <tex>n^d</tex> <ref name="while">
While, L., Hingston, P., Barone, L., Huband, S.: A Faster Algorithm for Calculating Hypervolume. IEEE Transaction on Evolutionary Computation, Vol.10, No. 1, pp 29–38 (2006)</ref>.

== Алгоритм Hypervolume by Slicing Objectives (HSO) ==

Под ''Objectives'' в названии данного алгоритма имеются в виду координаты пространства <tex>\mathbb{R}^d</tex>.

Если алгоритм LebMeasure по очереди рассматривает все точки, то алгоритм HSO<ref>Zhou X., Sun C., Mao N., Li W.: Generalization of HSO algortihms for computing hypervolume for mutiobjective optimization problems, IEEE Congress on Evolutionary Computation. (2007)</ref> рассматривает по очереди все координаты, сводя задачу к меньшей на единицу размерности.

Изначально все точки сортируются по первой координате. Значения этой координаты используются для ''расслоения'' (разрезания) всего множества на <tex>n</tex> частей, внутри каждой из которых при движении вдоль первой координаты форма ''разреза'' (гиперплоскости), проходящего перпендикулярно оси первой координаты, не меняется. Таким образом, для подсчета объема каждой части необходимо найти объем разреза, и умножить на длину части вдоль первой координаты. При этом получившийся разрез имеет на единицу меньшую размерность.
Заметим, что после сортировки и расслоения первая часть содержит все исходные <tex>n</tex> точек множества <tex>X</tex>, вторая - все, помимо точки с минимальной первой координатой, вдоль которой происходит расслоение, и т.д., а последняя часть содержит только одну точку с максимальной этой координатой.

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

Описанный процесс можно реализовать как в виде рекурсивной процедуры, расслаивающей множество вдоль очередной координаты и вызывающей себя рекурсивно для каждой части и следующей координаты, так и нерекурсивно, если поддерживать множество всех текущих частей и на очередной итерации разбивать их все вдоль очередной координаты.

Время работы алгоритма напрямую зависит от суммарного количества частей, на которые будет разбито исходное множество. По аналогичным предыдущему алгоритму рассуждениям легко показывается, что частей не более <tex>n^d</tex>. Тем не менее, существует более точная оценка.

{{Теорема
|statement=Суммарное количество частей, полученных алгоритмом HSO, не превышает <center><tex>f(n, d) = {n + d - 2 \choose d - 1}</tex>.</center>

|proof=
База:<center><tex>f(n, 1) = 1</tex>.</center>

Также, из алгоритма расслоения следует, что из части с <tex>n</tex> точками будет получено по одной части для каждого <tex>i = 1...n</tex>, то есть: <center><tex>f(n, d) = \sum \limits_{i = 1}^{n}f(i, d - 1)</tex>.</center>

Докажем полезные леммы:

{{Лемма
|statement=<center><tex>{x \choose y} = {x - r \choose y} + \sum \limits_{i = 1}^{r} {x - i \choose y - 1}</tex>.</center>
|proof=Индукция по <tex>r</tex>. База: <center><tex>{x \choose y} = {x - 1 \choose y} + {x - 1 \choose y - 1}</tex>.</center>
Переход: предполагаем, что утверждение верно для <tex>r = j</tex>. Раскрываем первое слагаемое в правой части предположения аналогично базе и получаем:
<center><tex>{x \choose y} = {x - j - 1 \choose y} + {x - (j + 1) \choose y - 1} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1} = {x - (j + 1) \choose y} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1}</tex>.</center>
}}

{{Лемма
|statement=<center><tex>{x \choose y} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
|proof=Подставляем в предыдущую лемму <tex>r = x - y</tex> и получаем:
<center><tex>{x \choose y} = {y \choose y} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = {x - (x - y + 1) \choose y - 1} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
}}

продолжаем доказательство теоремы по индукции:
<center><tex>f(n, d) = {n + d - 2 \choose d - 1} = \sum \limits_{i = 1}^{n} {n + d - 2 - i \choose d - 2} = \sum \limits_{k = n}^{1} {n + d + 2 - (n - k + 1) \choose d - 2} = \sum \limits_{k = 1}^{n} {k + d - 3 \choose d - 2} = \sum \limits_{k = 1}^{n}{k + (d - 1) - 2 \choose (d - 1) - 1} = \sum \limits_{k = 1}^{n}{f(k, d - 1)}</tex></center>

}}

Таким образом, на этот алгоритм получена более низкая оценка времени работы, чем на предыдущий, и это подтверждается на практике <ref name ="while"/>.

== Сведение к задаче KMP (Klee's Measure Problem) ==

Задача KMP состоит в нахождении объема объединения прямоугольных гиперпараллелепипедов в <tex>d</tex>-мерном пространстве. Как показано в описании алгоритма IEA, исходная задача легко сводится к этой, если каждой точке поставить в соответствие гиперкуб с одной вершиной в центре координат, а противоположной - в этой точке.

Существует несколько алгоритмов решения задачи KMP<ref name = "overmars">Overmars M.H., Yap C.K.: New upper bounds in Klee’s measure problem. SIAM Journal on Computing 20(6), pp 1034–1045. (1991)</ref><ref>Beume N.: S-Metric calculation by considering dominated hypervolume as Klee’s measure problem. Evolutionary Computation, 17(4), pp 477–492. (2009)</ref><ref>Beume N., Rudolph G.: Faster S-metric calculation by considering dominated hypervolume as Klee’s measure problem. In Proc. Second International Conference on Computational Intelligence (IASTED ’06), pp 233–238. (2006)</ref>, самый оптимальный из которых использует идеи сканирующей гиперплоскости, [[Статистики на отрезках. Корневая эвристика|корневой эвристики]] и КД-дерево, и позволяет решать задачу за время <tex>O(n^{d/2}\log n)</tex>. Рассмотрим его.

Для начала изложим идею сканирующей гиперплоскости. Рассмотрим всевозможные значения <tex>d</tex>-ой координаты у точек множества <tex>X</tex>.
Возьмем гиперплоскость, перпендикулярную оси <tex>d</tex>-ой координаты, будем ее двигать по этим значениям от минимального до максимального и рассматривать проекцию всех гиперкубов на эту гиперплоскость. Заметим, что между точками проекция меняться не будет, а в каждой точке будут появляться или исчезать проекции некоторых гиперкубов. Поэтому, если уметь эффективно пересчитывать объем проекции при каждом появлении/исчезновении и прибавлять это значение, умноженное на расстояние между последовательными значениями <tex>d</tex>-ой координаты, то можно легко посчитать суммарный гиперобъем. Как раз для эффективного пересчета используется такая структура, как КД-дерево.

Теперь изложим идею алгоритма на примере трехмерного пространства.
В трехмерном пространстве сканирующая гиперплоскость превращается в обычную плоскость. Будем считать, что она всегда располагается перпендикулярно оси <tex>OZ</tex> и двигается от минимальных значений по этой оси до максимальных. В каждый момент проекция всех параллелепипедов на эту плоскость представляет собой множество прямоугольников и необходимо уметь обрабатывать следующие три запроса:
# добавить прямоугольник
# удалить прямоугольник
# посчитать площадь объединения всех прямоугольников

{{Определение
|definition=КД-деревом называется [[Дерево поиска, наивная реализация|сбалансированное двоичное дерево]], где каждой вершине <tex>\alpha</tex> соответствует некоторая область <tex>C_{\alpha}</tex> <tex>d</tex>-мерного пространства, удовлетворяющая следующим свойствам:
# Корневой вершине <tex>root</tex> соответствует <tex>C_{root}</tex> - все пространство.
# Каждое <tex>C_{\alpha}</tex> является гиперпараллелепипедом (возможно бесконечным в некоторых направлениях).
# Для каждой вершины <tex>\alpha</tex> с детьми <tex>\beta</tex> и <tex>\gamma</tex> верно, что области <tex>C_{\beta}</tex> и <tex>C_{\gamma}</tex> не пересекаются и <tex>C_{\alpha} = C_{\beta} \cup C_{\gamma}</tex>.
}}

Будем хранить прямоугольники в КД-дереве следующим образом:
# Для каждого листа <tex>\alpha</tex> будем хранить все прямоугольники, которые пересекают внутреннюю часть области <tex>C_{\alpha}</tex>.
# Для каждой внутренней вершины <tex>\alpha</tex> будем считать число <tex>Cover_{\alpha}</tex> - количество прямоугольников, которые полностью покрывают область <tex>C_{\alpha}</tex>, но не полностью покрывают область родителя этой вершины <tex>C_{father(\alpha)}</tex>.

Также для каждой вершины будем хранить значение площади пересечения области этой вершины со всеми прямоугольниками <tex>M</tex>, которое определяется следующим образом:
# Для каждого листа <tex>M</tex> равно площади пересечения всех прямоугольников с областью этого листа.
# Для каждой внутренней вершины <tex>\alpha</tex> если <tex>Cover_{\alpha} > 0</tex>, тогда <tex>M</tex> равно площади всей области <tex>C_{\alpha}</tex>, иначе <tex>M</tex> равно сумме этих значений для двух сыновей этой вершины.

Таким образом, значение для корневой вершины <tex>M_{root}</tex> и является искомой площадью объединения всех прямоугольников.

Рассмотрим операцию добавления прямоугольника (box) в КД-дерево.

procedure Insert(box, <tex>\alpha</tex>)
if <tex>\alpha</tex> - лист then
добавить box в эту вершину
пересчитать <tex>M_{\alpha}</tex>
elseif box покрывает область <tex>C_{\alpha}</tex> then
<tex>Cover_{\alpha} = Cover_{\alpha} + 1</tex>
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
elseif box пересекает область <tex>C_{\alpha}</tex> then
Insert(box, leftson(<tex>\alpha</tex>))
Insert(box, rightson(<tex>\alpha</tex>))
if <tex>Cover_{\alpha} > 0 </tex> then
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
else
<tex>M_{\alpha} = M_{leftson(\alpha)} + M_{rightson(\alpha)}</tex>

Аналогично устроена операция удаления прямоугольника - она выполняет обратные действия.

Рассмотрим теперь как именно будет устроено разделение пространства на области КД-деревом для эффективного его использования.
Возьмем проекции всех параллелепипедов на плоскость, перпендикулярную оси <tex>OZ</tex> - множество прямоугольников <tex>V</tex>. Рассмотрим все <tex>x</tex>-координаты границ прямоугольников и разобьем ось <tex>OX</tex> на <tex>2\sqrt{n}</tex> промежутков так, чтобы в каждом из них было не более <tex>\sqrt{n}</tex> границ. Получившиеся <tex>2\sqrt{n}</tex> полос разобьем горизонтальными линиями на итоговые области по отдельности следующим образом. Для каждой полосы у всех прямоугольников, имеющих вертикальную границу в этой полосе, возьмем все горизонтальные границы (их будет не более <tex>2\sqrt{n}</tex>), а среди всех горизонтальных границ, пересекающих область, возьмем каждую <tex>\sqrt{n}</tex>-ую. Таким образом, в итоге получим не более <tex>4\sqrt{n}</tex> областей в каждой вертикальной полосе.

Полученные области обладают следующими свойствами:
# Всего <tex>O(n)</tex> областей.
# Каждый прямоугольник частично покрывает не более <tex>O(\sqrt{n})</tex> областей.
# Никакая вершина прямоугольника не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти утверждения несложно следуют из алгоритма разбиения на области.

Каждая из этих областей будет одним из листов КД-дерева. Для построения всего дерева сначала объединим все области для одной вертикальной полосы, сливая соседние области в одну вершину дерева. После этого аналогично будем сливать соседние вертикальные области и в итоге получим корневую вершину, представляющую всю плоскость.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n)</tex> вершин.
# Каждый прямоугольник хранится в <tex>O(\sqrt{n})</tex> листьях.
# Каждый прямоугольник влияет на <tex>O(\sqrt{n} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина прямоугольника не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти свойства следуют из свойств разбиения на области и логарифмической высоты дерева.

Наконец, получаем итоговый алгоритм. Идем сканирующей плоскостью, перпендикулярной оси <tex>OZ</tex>, добавляя или удаляя прямоугольники из КД-дерева. При этом каждый раз пересчитываем площадь объединения всех текущих прямоугольников и, складывая эти величины, умноженные на соответствующую разницу <tex>z</tex>-координат, получаем общий объем.

Каждая операция добавления или удаления прямоугольников требует изменения <tex>O(\sqrt{n} \log n)</tex> полей <tex>Cover</tex>, а также операций добавления или удаления прямоугольника в <tex>O(\sqrt{n})</tex> листьев, каждая из которых требует <tex>O(\log n)</tex> времени. Итого, на одну операцию добавления или удаления требуется <tex>O(\sqrt{n} \log n)</tex> времени, суммарно таких операций <tex>O(n)</tex>, и получается, что весь алгоритм работает за <tex>O(n \sqrt{n} \log n)</tex>. При этом на хранение всей информации требуется <tex>O(n\sqrt{n}) </tex> памяти - каждый лист может хранить <tex>O(\sqrt{n})</tex> прямоугольников.

Все идеи этого алгоритма обобщаются с трехмерного пространства на <tex>d</tex>-мерное следующим образом.

Построение разбиения на области в <tex>d</tex>-мерном пространстве вместо двухмерного ведется сначала аналогичным образом, только после того, как в двухмерном пространстве получаются итоговые прямоугольные области, здесь продолжаются аналогичные разбиения вдоль третьей координаты - берутся все границы третьей координаты гиперпараллелепипедов, имеющие границу по первой или второй координате строго внутри области (их получается не более <tex>O(\sqrt{n})</tex>), а также каждую <tex>O(\sqrt{n})</tex>-ую границу по третьей координате гиперепараллелепипедов, пересекающих эту область, и так далее.

Это разбиение на области имеет следующие аналогичные свойства:
# Всего <tex>O(n^{d/2})</tex> областей.
# Каждый гиперпараллелепипед частично покрывает не более <tex>O(n^{(d-1)/2})</tex> областей.
# Никакая вершина гиперпараллелепипеда не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Общее дерево строится аналогичным образом и оно в итоге состоит из <tex>d</tex> уровней, каждый из которых имеет высоту <tex>O(\log n)</tex> в дереве - самый верхний уровень содержит разбиение по первой координате, следующий - по второй и так далее.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n^{d/2})</tex> вершин.
# Каждый гиперпараллелепипед хранится в <tex>O(n^{(d-1)/2})</tex> листьях.
# Каждый гиперпараллелепипед влияет на <tex>O(n^{(d-1)/2} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина гиперпараллелепипеда не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Операции в этом дереве выполняются аналогично трехмерному случаю и каждая операция добавления или удаления в КД-дерево размерности <tex>d</tex> занимает <tex>O(n^{(d-1)/2} \log n)</tex> времени. Выполняя сканирование гиперплоскостью вдоль одной из координат получаем КД-дерево размерности <tex>d-1</tex>, в котором необходимо провести <tex>O(n)</tex> операций, то есть суммарное время работы алгоритма равно: <tex>O(n \times n^{(d-1-1)/2} \log n) = O(n^{d/2} \log n)</tex>. При этом хранение КД-дерева требует <tex>O(n^{d/2})</tex> памяти.

== Источники ==
<references/>

Алгоритмы точного вычисления гиперобъема

2012-06-20T12:01:57Z

Joshik:

{{В разработке}}

== Постановка задачи ==
<tex>x = (x_1, x_2, ..., x_d; x_i \ge 0) \in \mathbb{R}^d</tex> - точка в <tex>d</tex>-мерном пространстве.

Точка <tex>x</tex> доминирует точку <tex>y</tex> (<tex>x \succ y</tex>), если <tex>\forall i : x_i \ge y_i, \exists j : x_j > y_j</tex>.

<tex>X = \{x^1, x^2, ..., x^n\} \subset \mathbb{R}^d</tex> - множество из <tex>n</tex> точек в <tex>d</tex>-мерном пространстве таких, что <tex>\nexists i \neq j : x_i \succ x_j</tex> - никакая точка не доминируется другой точкой из этого множества.

<tex>S(X) = \mu (\bigcup \limits_{x \in X} \{y | x \succ y\}) </tex> - гиперобъем множества <tex>X</tex>, где <tex>\mu</tex> - мера Лебега.

В частности, если <tex>X = \{x\}</tex>, то <tex>S(X) = \prod \limits_{i=1}^{d} x_i</tex>.

Задача: найти точное значение гиперобъема <tex>S(X)</tex> множества из <tex>n</tex> точек <tex>d</tex>-мерного пространства.

== Алгоритм включения-исключения (Inclusion-Exclusion Algorithm, IEA) ==

Самый простой алгоритм нахождения гиперобъема базируется на идее комбинаторной [[формула включения-исключения|формулы включения-искючения]].
Все множество <tex>X</tex> представляется в виде объединения <tex>n</tex> гиперкубов (<tex>X^i</tex>), соответствующих отдельным точкам <tex>x^i</tex>.

После этого объем всего множества вычисляется по формуле: <center> <tex> S(X) = \sum \limits_{I \in 2^n} (-1)^{|I|+1} S(\bigcap \limits_{ j \in I} X^j) </tex> </center>

Объем пересечения гиперкубов легко считается как произведение по каждой координате минимального значения этой координаты среди всех точек, которым соответствуют гиперкубы.

Таким образом, в этом алгоритме перебираются все подмножества точек множества <tex>X</tex>, для каждого множества находится гиперобъем пересечения соответствующих гиперкубов, и он прибавляется с соответствующим знаком к ответу.
Время работы этого алгоритма составляет <tex>O(n 2^n)</tex>.

== Алгоритм LebMeasure ==

Алгоритм LebMeasure<ref>Bringmann K., Friedrich. T.: An Efficient Algorithm for Computing Hypervolume Contributions. Evolutionary Computation, Vol. 18, No. 3, pp 383-402, MIT Press. (2010)</ref> обрабатывает точки множества <tex>X</tex> по очереди. Для каждой очередной точки <tex>x</tex> находится объем некоторого максимального по включению гиперкуба, доминируемого эксклюзивно только этой точкой <tex>x</tex> и она заменяется на некоторое множество ''порожденных'' точек, которые доминируют оставшуюся область, доминировавшуюся до этого точкой <tex>x</tex>.

Например, если изначально было четыре точки в трехмерном пространстве <tex>(6, 7, 4), (9, 5, 5), (1, 9, 3), (4, 1, 9)</tex>, то точка <tex>(6, 7, 4)</tex> эксклюзивно доминирует куб с одним концов в <tex>(6, 7, 4)</tex>, а другим - в <tex>(4, 5, 3)</tex>. После добавления объема этого куба к ответу, точка <tex>(6, 7, 4)</tex> порождает три точки: <tex>(4, 7, 4), (6, 5, 3), (6, 7, 3)</tex>. При этом точка <tex>(6, 5, 4)</tex> доминируется точкой <tex>(9, 5, 5)</tex> и сразу удаляется из множества <tex>X</tex>.

Обработка точек продолжается, пока все точки не будут обработаны. Таким образом, время работы алгоритма напрямую зависит от количества порожденных точек. Легко заметить, что таких точек не более, чем <tex>n^d</tex>, поскольку каждая координата каждой порожденной точки равна соответствующей координате некоторой точки исходного множества <tex>X</tex>.

К сожалению, эта верхняя оценка является достижимой. Например, если исходное множество <tex>X</tex> имеет вид: <tex>(1, n, n, ..., n), (2, n - 1, n - 1, ..., n - 1), ..., (n, 1, 1, ..., 1)</tex>, и точки обрабатываются именно в этом порядке, то всего будет обработано <tex>n^{d-1}</tex> точка, что показано в <ref name = "while"/>. Правда, если в этом примере точки обрабатываются в обратном порядке, то суммарное количество обработанных точек линейно зависит от <tex>n</tex> и <tex>d</tex>. Тем не менее, существуют примеры, для которых любой порядок обработки приводит к экспоненциальной зависимости числа порожденных точек от размерности пространства <tex>d</tex> и эта зависимость близка к <tex>n^d</tex> <ref name="while">
While, L., Hingston, P., Barone, L., Huband, S.: A Faster Algorithm for Calculating Hypervolume. IEEE Transaction on Evolutionary Computation, Vol.10, No. 1, pp 29–38 (2006)</ref>.

== Алгоритм Hypervolume by Slicing Objectives (HSO) ==

Под ''Objectives'' в названии данного алгоритма имеются в виду координаты пространства <tex>\mathbb{R}^d</tex>.

Если алгоритм LebMeasure по очереди рассматривает все точки, то алгоритм HSO<ref>Zhou X., Sun C., Mao N., Li W.: Generalization of HSO algortihms for computing hypervolume for mutiobjective optimization problems, IEEE Congress on Evolutionary Computation (2007)</ref> рассматривает по очереди все координаты, сводя задачу к меньшей на единицу размерности.

Изначально все точки сортируются по первой координате. Значения этой координаты используются для ''расслоения'' (разрезания) всего множества на <tex>n</tex> частей, внутри каждой из которых при движении вдоль первой координаты форма ''разреза'' (гиперплоскости), проходящего перпендикулярно оси первой координаты, не меняется. Таким образом, для подсчета объема каждой части необходимо найти объем разреза, и умножить на длину части вдоль первой координаты. При этом получившийся разрез имеет на единицу меньшую размерность.
Заметим, что после сортировки и расслоения первая часть содержит все исходные <tex>n</tex> точек множества <tex>X</tex>, вторая - все, помимо точки с минимальной первой координатой, вдоль которой происходит расслоение, и т.д., а последняя часть содержит только одну точку с максимальной этой координатой.

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

Описанный процесс можно реализовать как в виде рекурсивной процедуры, расслаивающей множество вдоль очередной координаты и вызывающей себя рекурсивно для каждой части и следующей координаты, так и нерекурсивно, если поддерживать множество всех текущих частей и на очередной итерации разбивать их все вдоль очередной координаты.

Время работы алгоритма напрямую зависит от суммарного количества частей, на которые будет разбито исходное множество. По аналогичным предыдущему алгоритму рассуждениям легко показывается, что частей не более <tex>n^d</tex>. Тем не менее, существует более точная оценка.

{{Теорема
|statement=Суммарное количество частей, полученных алгоритмом HSO, не превышает <center><tex>f(n, d) = {n + d - 2 \choose d - 1}</tex>.</center>

|proof=
База:<center><tex>f(n, 1) = 1</tex>.</center>

Также, из алгоритма расслоения следует, что из части с <tex>n</tex> точками будет получено по одной части для каждого <tex>i = 1...n</tex>, то есть: <center><tex>f(n, d) = \sum \limits_{i = 1}^{n}f(i, d - 1)</tex>.</center>

Докажем полезные леммы:

{{Лемма
|statement=<center><tex>{x \choose y} = {x - r \choose y} + \sum \limits_{i = 1}^{r} {x - i \choose y - 1}</tex>.</center>
|proof=Индукция по <tex>r</tex>. База: <center><tex>{x \choose y} = {x - 1 \choose y} + {x - 1 \choose y - 1}</tex>.</center>
Переход: предполагаем, что утверждение верно для <tex>r = j</tex>. Раскрываем первое слагаемое в правой части предположения аналогично базе и получаем:
<center><tex>{x \choose y} = {x - j - 1 \choose y} + {x - (j + 1) \choose y - 1} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1} = {x - (j + 1) \choose y} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1}</tex>.</center>
}}

{{Лемма
|statement=<center><tex>{x \choose y} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
|proof=Подставляем в предыдущую лемму <tex>r = x - y</tex> и получаем:
<center><tex>{x \choose y} = {y \choose y} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = {x - (x - y + 1) \choose y - 1} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
}}

продолжаем доказательство теоремы по индукции:
<center><tex>f(n, d) = {n + d - 2 \choose d - 1} = \sum \limits_{i = 1}^{n} {n + d - 2 - i \choose d - 2} = \sum \limits_{k = n}^{1} {n + d + 2 - (n - k + 1) \choose d - 2} = \sum \limits_{k = 1}^{n} {k + d - 3 \choose d - 2} = \sum \limits_{k = 1}^{n}{k + (d - 1) - 2 \choose (d - 1) - 1} = \sum \limits_{k = 1}^{n}{f(k, d - 1)}</tex></center>

}}

Таким образом, на этот алгоритм получена более низкая оценка времени работы, чем на предыдущий, и это подтверждается на практике <ref name ="while"/>.

== Сведение к задаче KMP (Klee's Measure Problem) ==

Задача KMP состоит в нахождении объема объединения прямоугольных гиперпараллелепипедов в <tex>d</tex>-мерном пространстве. Как показано в описании алгоритма IEA, исходная задача легко сводится к этой, если каждой точке поставить в соответствие гиперкуб с одной вершиной в центре координат, а противоположной - в этой точке.

Существует несколько алгоритмов решения задачи KMP<ref name = "overmars">Overmars M.H., Yap C.K.: New upper bounds in Klee’s measure problem. SIAM Journal on Computing 20(6), pp 1034–1045. (1991)</ref><ref>Beume N.: S-Metric calculation by considering dominated hypervolume as Klee’s measure problem. Evolutionary Computation, 17(4), pp 477–492. (2009)</ref><ref>Beume N., Rudolph G.: Faster S-metric calculation by considering dominated hypervolume as Klee’s measure problem. In Proc. Second International Conference on Computational Intelligence (IASTED ’06), pp 233–238. (2006)
</ref>, самый оптимальный из которых использует идеи сканирующей гиперплоскости, [[Статистики на отрезках. Корневая эвристика|корневой эвристики]] и КД-дерево, и позволяет решать задачу за время <tex>O(n^{d/2}\log n)</tex>. Рассмотрим его.

Для начала изложим идею сканирующей гиперплоскости. Рассмотрим всевозможные значения <tex>d</tex>-ой координаты у точек множества <tex>X</tex>.
Возьмем гиперплоскость, перпендикулярную оси <tex>d</tex>-ой координаты, будем ее двигать по этим значениям от минимального до максимального и рассматривать проекцию всех гиперкубов на эту гиперплоскость. Заметим, что между точками проекция меняться не будет, а в каждой точке будут появляться или исчезать проекции некоторых гиперкубов. Поэтому, если уметь эффективно пересчитывать объем проекции при каждом появлении/исчезновении и прибавлять это значение, умноженное на расстояние между последовательными значениями <tex>d</tex>-ой координаты, то можно легко посчитать суммарный гиперобъем. Как раз для эффективного пересчета используется такая структура, как КД-дерево.

Теперь изложим идею алгоритма на примере трехмерного пространства.
В трехмерном пространстве сканирующая гиперплоскость превращается в обычную плоскость. Будем считать, что она всегда располагается перпендикулярно оси <tex>OZ</tex> и двигается от минимальных значений по этой оси до максимальных. В каждый момент проекция всех параллелепипедов на эту плоскость представляет собой множество прямоугольников и необходимо уметь обрабатывать следующие три запроса:
# добавить прямоугольник
# удалить прямоугольник
# посчитать площадь объединения всех прямоугольников

{{Определение
|definition=КД-деревом называется [[Дерево поиска, наивная реализация|сбалансированное двоичное дерево]], где каждой вершине <tex>\alpha</tex> соответствует некоторая область <tex>C_{\alpha}</tex> <tex>d</tex>-мерного пространства, удовлетворяющая следующим свойствам:
# Корневой вершине <tex>root</tex> соответствует <tex>C_{root}</tex> - все пространство.
# Каждое <tex>C_{\alpha}</tex> является гиперпараллелепипедом (возможно бесконечным в некоторых направлениях).
# Для каждой вершины <tex>\alpha</tex> с детьми <tex>\beta</tex> и <tex>\gamma</tex> верно, что области <tex>C_{\beta}</tex> и <tex>C_{\gamma}</tex> не пересекаются и <tex>C_{\alpha} = C_{\beta} \cup C_{\gamma}</tex>.
}}

Будем хранить прямоугольники в КД-дереве следующим образом:
# Для каждого листа <tex>\alpha</tex> будем хранить все прямоугольники, которые пересекают внутреннюю часть области <tex>C_{\alpha}</tex>.
# Для каждой внутренней вершины <tex>\alpha</tex> будем считать число <tex>Cover_{\alpha}</tex> - количество прямоугольников, которые полностью покрывают область <tex>C_{\alpha}</tex>, но не полностью покрывают область родителя этой вершины <tex>C_{father(\alpha)}</tex>.

Также для каждой вершины будем хранить значение площади пересечения области этой вершины со всеми прямоугольниками <tex>M</tex>, которое определяется следующим образом:
# Для каждого листа <tex>M</tex> равно площади пересечения всех прямоугольников с областью этого листа.
# Для каждой внутренней вершины <tex>\alpha</tex> если <tex>Cover_{\alpha} > 0</tex>, тогда <tex>M</tex> равно площади всей области <tex>C_{\alpha}</tex>, иначе <tex>M</tex> равно сумме этих значений для двух сыновей этой вершины.

Таким образом, значение для корневой вершины <tex>M_{root}</tex> и является искомой площадью объединения всех прямоугольников.

Рассмотрим операцию добавления прямоугольника (box) в КД-дерево.

procedure Insert(box, <tex>\alpha</tex>)
if <tex>\alpha</tex> - лист then
добавить box в эту вершину
пересчитать <tex>M_{\alpha}</tex>
elseif box покрывает область <tex>C_{\alpha}</tex> then
<tex>Cover_{\alpha} = Cover_{\alpha} + 1</tex>
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
elseif box пересекает область <tex>C_{\alpha}</tex> then
Insert(box, leftson(<tex>\alpha</tex>))
Insert(box, rightson(<tex>\alpha</tex>))
if <tex>Cover_{\alpha} > 0 </tex> then
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
else
<tex>M_{\alpha} = M_{leftson(\alpha)} + M_{rightson(\alpha)}</tex>

Аналогично устроена операция удаления прямоугольника - она выполняет обратные действия.

Рассмотрим теперь как именно будет устроено разделение пространства на области КД-деревом для эффективного его использования.
Возьмем проекции всех параллелепипедов на плоскость, перпендикулярную оси <tex>OZ</tex> - множество прямоугольников <tex>V</tex>. Рассмотрим все <tex>x</tex>-координаты границ прямоугольников и разобьем ось <tex>OX</tex> на <tex>2\sqrt{n}</tex> промежутков так, чтобы в каждом из них было не более <tex>\sqrt{n}</tex> границ. Получившиеся <tex>2\sqrt{n}</tex> полос разобьем горизонтальными линиями на итоговые области по отдельности следующим образом. Для каждой полосы у всех прямоугольников, имеющих вертикальную границу в этой полосе, возьмем все горизонтальные границы (их будет не более <tex>2\sqrt{n}</tex>), а среди всех горизонтальных границ, пересекающих область, возьмем каждую <tex>\sqrt{n}</tex>-ую. Таким образом, в итоге получим не более <tex>4\sqrt{n}</tex> областей в каждой вертикальной полосе.

Полученные области обладают следующими свойствами:
# Всего <tex>O(n)</tex> областей.
# Каждый прямоугольник частично покрывает не более <tex>O(\sqrt{n})</tex> областей.
# Никакая вершина прямоугольника не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти утверждения несложно следуют из алгоритма разбиения на области.

Каждая из этих областей будет одним из листов КД-дерева. Для построения всего дерева сначала объединим все области для одной вертикальной полосы, сливая соседние области в одну вершину дерева. После этого аналогично будем сливать соседние вертикальные области и в итоге получим корневую вершину, представляющую всю плоскость.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n)</tex> вершин.
# Каждый прямоугольник хранится в <tex>O(\sqrt{n})</tex> листьях.
# Каждый прямоугольник влияет на <tex>O(\sqrt{n} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина прямоугольника не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти свойства следуют из свойств разбиения на области и логарифмической высоты дерева.

Наконец, получаем итоговый алгоритм. Идем сканирующей плоскостью, перпендикулярной оси <tex>OZ</tex>, добавляя или удаляя прямоугольники из КД-дерева. При этом каждый раз пересчитываем площадь объединения всех текущих прямоугольников и, складывая эти величины, умноженные на соответствующую разницу <tex>z</tex>-координат, получаем общий объем.

Каждая операция добавления или удаления прямоугольников требует изменения <tex>O(\sqrt{n} \log n)</tex> полей <tex>Cover</tex>, а также операций добавления или удаления прямоугольника в <tex>O(\sqrt{n})</tex> листьев, каждая из которых требует <tex>O(\log n)</tex> времени. Итого, на одну операцию добавления или удаления требуется <tex>O(\sqrt{n} \log n)</tex> времени, суммарно таких операций <tex>O(n)</tex>, и получается, что весь алгоритм работает за <tex>O(n \sqrt{n} \log n)</tex>. При этом на хранение всей информации требуется <tex>O(n\sqrt{n}) </tex> памяти - каждый лист может хранить <tex>O(\sqrt{n})</tex> прямоугольников.

Все идеи этого алгоритма обобщаются с трехмерного пространства на <tex>d</tex>-мерное следующим образом.

Построение разбиения на области в <tex>d</tex>-мерном пространстве вместо двухмерного ведется сначала аналогичным образом, только после того, как в двухмерном пространстве получаются итоговые прямоугольные области, здесь продолжаются аналогичные разбиения вдоль третьей координаты - берутся все границы третьей координаты гиперпараллелепипедов, имеющие границу по первой или второй координате строго внутри области (их получается не более <tex>O(\sqrt{n})</tex>), а также каждую <tex>O(\sqrt{n})</tex>-ую границу по третьей координате гиперепараллелепипедов, пересекающих эту область, и так далее.

Это разбиение на области имеет следующие аналогичные свойства:
# Всего <tex>O(n^{d/2})</tex> областей.
# Каждый гиперпараллелепипед частично покрывает не более <tex>O(n^{(d-1)/2})</tex> областей.
# Никакая вершина гиперпараллелепипеда не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Общее дерево строится аналогичным образом и оно в итоге состоит из <tex>d</tex> уровней, каждый из которых имеет высоту <tex>O(\log n)</tex> в дереве - самый верхний уровень содержит разбиение по первой координате, следующий - по второй и так далее.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n^{d/2})</tex> вершин.
# Каждый гиперпараллелепипед хранится в <tex>O(n^{(d-1)/2})</tex> листьях.
# Каждый гиперпараллелепипед влияет на <tex>O(n^{(d-1)/2} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина гиперпараллелепипеда не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Операции в этом дереве выполняются аналогично трехмерному случаю и каждая операция добавления или удаления в КД-дерево размерности <tex>d</tex> занимает <tex>O(n^{(d-1)/2} \log n)</tex> времени. Выполняя сканирование гиперплоскостью вдоль одной из координат получаем КД-дерево размерности <tex>d-1</tex>, в котором необходимо провести <tex>O(n)</tex> операций, то есть суммарное время работы алгоритма равно: <tex>O(n \times n^{(d-1-1)/2} \log n) = O(n^{d/2} \log n)</tex>. При этом хранение КД-дерева требует <tex>O(n^{d/2})</tex> памяти.

== Источники ==
<references/>

Алгоритмы точного вычисления гиперобъема

2012-06-20T11:58:40Z

Joshik:

{{В разработке}}

== Постановка задачи ==
<tex>x = (x_1, x_2, ..., x_d; x_i \ge 0) \in \mathbb{R}^d</tex> - точка в <tex>d</tex>-мерном пространстве.

Точка <tex>x</tex> доминирует точку <tex>y</tex> (<tex>x \succ y</tex>), если <tex>\forall i : x_i \ge y_i, \exists j : x_j > y_j</tex>.

<tex>X = \{x^1, x^2, ..., x^n\} \subset \mathbb{R}^d</tex> - множество из <tex>n</tex> точек в <tex>d</tex>-мерном пространстве таких, что <tex>\nexists i \neq j : x_i \succ x_j</tex> - никакая точка не доминируется другой точкой из этого множества.

<tex>S(X) = \mu (\bigcup \limits_{x \in X} \{y | x \succ y\}) </tex> - гиперобъем множества <tex>X</tex>, где <tex>\mu</tex> - мера Лебега.

В частности, если <tex>X = \{x\}</tex>, то <tex>S(X) = \prod \limits_{i=1}^{d} x_i</tex>.

Задача: найти точное значение гиперобъема <tex>S(X)</tex> множества из <tex>n</tex> точек <tex>d</tex>-мерного пространства.

== Алгоритм включения-исключения (Inclusion-Exclusion Algorithm, IEA) ==

Самый простой алгоритм нахождения гиперобъема базируется на идее комбинаторной [[формула включения-исключения|формулы включения-искючения]].
Все множество <tex>X</tex> представляется в виде объединения <tex>n</tex> гиперкубов (<tex>X^i</tex>), соответствующих отдельным точкам <tex>x^i</tex>.

После этого объем всего множества вычисляется по формуле: <center> <tex> S(X) = \sum \limits_{I \in 2^n} (-1)^{|I|+1} S(\bigcap \limits_{ j \in I} X^j) </tex> </center>

Объем пересечения гиперкубов легко считается как произведение по каждой координате минимального значения этой координаты среди всех точек, которым соответствуют гиперкубы.

Таким образом, в этом алгоритме перебираются все подмножества точек множества <tex>X</tex>, для каждого множества находится гиперобъем пересечения соответствующих гиперкубов, и он прибавляется с соответствующим знаком к ответу.
Время работы этого алгоритма составляет <tex>O(n 2^n)</tex>.

== Алгоритм LebMeasure ==

Алгоритм LebMeasure<ref>K. Bringmann, T. Friedrich.: An Efficient Algorithm for Computing Hypervolume Contributions. Evolutionary Computation, Vol. 18, No. 3, pp 383-402, MIT Press. (2010)</ref> обрабатывает точки множества <tex>X</tex> по очереди. Для каждой очередной точки <tex>x</tex> находится объем некоторого максимального по включению гиперкуба, доминируемого эксклюзивно только этой точкой <tex>x</tex> и она заменяется на некоторое множество ''порожденных'' точек, которые доминируют оставшуюся область, доминировавшуюся до этого точкой <tex>x</tex>.

Например, если изначально было четыре точки в трехмерном пространстве <tex>(6, 7, 4), (9, 5, 5), (1, 9, 3), (4, 1, 9)</tex>, то точка <tex>(6, 7, 4)</tex> эксклюзивно доминирует куб с одним концов в <tex>(6, 7, 4)</tex>, а другим - в <tex>(4, 5, 3)</tex>. После добавления объема этого куба к ответу, точка <tex>(6, 7, 4)</tex> порождает три точки: <tex>(4, 7, 4), (6, 5, 3), (6, 7, 3)</tex>. При этом точка <tex>(6, 5, 4)</tex> доминируется точкой <tex>(9, 5, 5)</tex> и сразу удаляется из множества <tex>X</tex>.

Обработка точек продолжается, пока все точки не будут обработаны. Таким образом, время работы алгоритма напрямую зависит от количества порожденных точек. Легко заметить, что таких точек не более, чем <tex>n^d</tex>, поскольку каждая координата каждой порожденной точки равна соответствующей координате некоторой точки исходного множества <tex>X</tex>.

К сожалению, эта верхняя оценка является достижимой. Например, если исходное множество <tex>X</tex> имеет вид: <tex>(1, n, n, ..., n), (2, n - 1, n - 1, ..., n - 1), ..., (n, 1, 1, ..., 1)</tex>, и точки обрабатываются именно в этом порядке, то всего будет обработано <tex>n^{d-1}</tex> точка, что показано в <ref name = "while"/>. Правда, если в этом примере точки обрабатываются в обратном порядке, то суммарное количество обработанных точек линейно зависит от <tex>n</tex> и <tex>d</tex>. Тем не менее, существуют примеры, для которых любой порядок обработки приводит к экспоненциальной зависимости числа порожденных точек от размерности пространства <tex>d</tex> и эта зависимость близка к <tex>n^d</tex> <ref name="while">
While, L., Hingston, P., Barone, L., Huband, S.: A Faster Algorithm for Calculating Hypervolume. IEEE Transaction on Evolutionary Computation, Vol.10, No. 1, pp 29–38 (2006)</ref>.

== Алгоритм Hypervolume by Slicing Objectives (HSO) ==

Под ''Objectives'' в названии данного алгоритма имеются в виду координаты пространства <tex>\mathbb{R}^d</tex>.

Если алгоритм LebMeasure по очереди рассматривает все точки, то алгоритм HSO<ref>Zhou X., Sun C., Mao N., Li W.: Generalization of HSO algortihms for computing hypervolume for mutiobjective optimization problems, IEEE Congress on Evolutionary Computation (2007)</ref> рассматривает по очереди все координаты, сводя задачу к меньшей на единицу размерности.

Изначально все точки сортируются по первой координате. Значения этой координаты используются для ''расслоения'' (разрезания) всего множества на <tex>n</tex> частей, внутри каждой из которых при движении вдоль первой координаты форма ''разреза'' (гиперплоскости), проходящего перпендикулярно оси первой координаты, не меняется. Таким образом, для подсчета объема каждой части необходимо найти объем разреза, и умножить на длину части вдоль первой координаты. При этом получившийся разрез имеет на единицу меньшую размерность.
Заметим, что после сортировки и расслоения первая часть содержит все исходные <tex>n</tex> точек множества <tex>X</tex>, вторая - все, помимо точки с минимальной первой координатой, вдоль которой происходит расслоение, и т.д., а последняя часть содержит только одну точку с максимальной этой координатой.

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

Описанный процесс можно реализовать как в виде рекурсивной процедуры, расслаивающей множество вдоль очередной координаты и вызывающей себя рекурсивно для каждой части и следующей координаты, так и нерекурсивно, если поддерживать множество всех текущих частей и на очередной итерации разбивать их все вдоль очередной координаты.

Время работы алгоритма напрямую зависит от суммарного количества частей, на которые будет разбито исходное множество. По аналогичным предыдущему алгоритму рассуждениям легко показывается, что частей не более <tex>n^d</tex>. Тем не менее, существует более точная оценка.

{{Теорема
|statement=Суммарное количество частей, полученных алгоритмом HSO, не превышает <center><tex>f(n, d) = {n + d - 2 \choose d - 1}</tex>.</center>

|proof=
База:<center><tex>f(n, 1) = 1</tex>.</center>

Также, из алгоритма расслоения следует, что из части с <tex>n</tex> точками будет получено по одной части для каждого <tex>i = 1...n</tex>, то есть: <center><tex>f(n, d) = \sum \limits_{i = 1}^{n}f(i, d - 1)</tex>.</center>

Докажем полезные леммы:

{{Лемма
|statement=<center><tex>{x \choose y} = {x - r \choose y} + \sum \limits_{i = 1}^{r} {x - i \choose y - 1}</tex>.</center>
|proof=Индукция по <tex>r</tex>. База: <center><tex>{x \choose y} = {x - 1 \choose y} + {x - 1 \choose y - 1}</tex>.</center>
Переход: предполагаем, что утверждение верно для <tex>r = j</tex>. Раскрываем первое слагаемое в правой части предположения аналогично базе и получаем:
<center><tex>{x \choose y} = {x - j - 1 \choose y} + {x - (j + 1) \choose y - 1} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1} = {x - (j + 1) \choose y} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1}</tex>.</center>
}}

{{Лемма
|statement=<center><tex>{x \choose y} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
|proof=Подставляем в предыдущую лемму <tex>r = x - y</tex> и получаем:
<center><tex>{x \choose y} = {y \choose y} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = {x - (x - y + 1) \choose y - 1} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
}}

продолжаем доказательство теоремы по индукции:
<center><tex>f(n, d) = {n + d - 2 \choose d - 1} = \sum \limits_{i = 1}^{n} {n + d - 2 - i \choose d - 2} = \sum \limits_{k = n}^{1} {n + d + 2 - (n - k + 1) \choose d - 2} = \sum \limits_{k = 1}^{n} {k + d - 3 \choose d - 2} = \sum \limits_{k = 1}^{n}{k + (d - 1) - 2 \choose (d - 1) - 1} = \sum \limits_{k = 1}^{n}{f(k, d - 1)}</tex></center>

}}

Таким образом, на этот алгоритм получена более низкая оценка времени работы, чем на предыдущий, и это подтверждается на практике <ref name ="while"/>.

== Сведение к задаче KMP (Klee's Measure Problem) ==

Задача KMP состоит в нахождении объема объединения прямоугольных гиперпараллелепипедов в <tex>d</tex>-мерном пространстве. Как показано в описании алгоритма IEA, исходная задача легко сводится к этой, если каждой точке поставить в соответствие гиперкуб с одной вершиной в центре координат, а противоположной - в этой точке.

Существует несколько алгоритмов решения задачи KMP<ref name = "overmars">Overmars, M.H., Yap, C.K.: New upper bounds in Klee’s measure problem. SIAM Journal on Computing 20(6), pp 1034–1045. (1991)</ref><ref>N. Beume.: S-Metric calculation by considering dominated hypervolume as Klee’s measure problem. Evolutionary Computation, 17(4), pp 477–492. (2009)</ref><ref>N. Beume, G. Rudolph.: Faster S-metric calculation by considering dominated hypervolume as Klee’s measure problem. In Proc. Second International Conference on Computational Intelligence (IASTED ’06), pp 233–238. (2006)
</ref>, самый оптимальный из которых использует идеи сканирующей гиперплоскости, [[Статистики на отрезках. Корневая эвристика|корневой эвристики]] и КД-дерево, и позволяет решать задачу за время <tex>O(n^{d/2}\log n)</tex>. Рассмотрим его.

Для начала изложим идею сканирующей гиперплоскости. Рассмотрим всевозможные значения <tex>d</tex>-ой координаты у точек множества <tex>X</tex>.
Возьмем гиперплоскость, перпендикулярную оси <tex>d</tex>-ой координаты, будем ее двигать по этим значениям от минимального до максимального и рассматривать проекцию всех гиперкубов на эту гиперплоскость. Заметим, что между точками проекция меняться не будет, а в каждой точке будут появляться или исчезать проекции некоторых гиперкубов. Поэтому, если уметь эффективно пересчитывать объем проекции при каждом появлении/исчезновении и прибавлять это значение, умноженное на расстояние между последовательными значениями <tex>d</tex>-ой координаты, то можно легко посчитать суммарный гиперобъем. Как раз для эффективного пересчета используется такая структура, как КД-дерево.

Теперь изложим идею алгоритма на примере трехмерного пространства.
В трехмерном пространстве сканирующая гиперплоскость превращается в обычную плоскость. Будем считать, что она всегда располагается перпендикулярно оси <tex>OZ</tex> и двигается от минимальных значений по этой оси до максимальных. В каждый момент проекция всех параллелепипедов на эту плоскость представляет собой множество прямоугольников и необходимо уметь обрабатывать следующие три запроса:
# добавить прямоугольник
# удалить прямоугольник
# посчитать площадь объединения всех прямоугольников

{{Определение
|definition=КД-деревом называется [[Дерево поиска, наивная реализация|сбалансированное двоичное дерево]], где каждой вершине <tex>\alpha</tex> соответствует некоторая область <tex>C_{\alpha}</tex> <tex>d</tex>-мерного пространства, удовлетворяющая следующим свойствам:
# Корневой вершине <tex>root</tex> соответствует <tex>C_{root}</tex> - все пространство.
# Каждое <tex>C_{\alpha}</tex> является гиперпараллелепипедом (возможно бесконечным в некоторых направлениях).
# Для каждой вершины <tex>\alpha</tex> с детьми <tex>\beta</tex> и <tex>\gamma</tex> верно, что области <tex>C_{\beta}</tex> и <tex>C_{\gamma}</tex> не пересекаются и <tex>C_{\alpha} = C_{\beta} \cup C_{\gamma}</tex>.
}}

Будем хранить прямоугольники в КД-дереве следующим образом:
# Для каждого листа <tex>\alpha</tex> будем хранить все прямоугольники, которые пересекают внутреннюю часть области <tex>C_{\alpha}</tex>.
# Для каждой внутренней вершины <tex>\alpha</tex> будем считать число <tex>Cover_{\alpha}</tex> - количество прямоугольников, которые полностью покрывают область <tex>C_{\alpha}</tex>, но не полностью покрывают область родителя этой вершины <tex>C_{father(\alpha)}</tex>.

Также для каждой вершины будем хранить значение площади пересечения области этой вершины со всеми прямоугольниками <tex>M</tex>, которое определяется следующим образом:
# Для каждого листа <tex>M</tex> равно площади пересечения всех прямоугольников с областью этого листа.
# Для каждой внутренней вершины <tex>\alpha</tex> если <tex>Cover_{\alpha} > 0</tex>, тогда <tex>M</tex> равно площади всей области <tex>C_{\alpha}</tex>, иначе <tex>M</tex> равно сумме этих значений для двух сыновей этой вершины.

Таким образом, значение для корневой вершины <tex>M_{root}</tex> и является искомой площадью объединения всех прямоугольников.

Рассмотрим операцию добавления прямоугольника (box) в КД-дерево.

procedure Insert(box, <tex>\alpha</tex>)
if <tex>\alpha</tex> - лист then
добавить box в эту вершину
пересчитать <tex>M_{\alpha}</tex>
elseif box покрывает область <tex>C_{\alpha}</tex> then
<tex>Cover_{\alpha} = Cover_{\alpha} + 1</tex>
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
elseif box пересекает область <tex>C_{\alpha}</tex> then
Insert(box, leftson(<tex>\alpha</tex>))
Insert(box, rightson(<tex>\alpha</tex>))
if <tex>Cover_{\alpha} > 0 </tex> then
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
else
<tex>M_{\alpha} = M_{leftson(\alpha)} + M_{rightson(\alpha)}</tex>

Аналогично устроена операция удаления прямоугольника - она выполняет обратные действия.

Рассмотрим теперь как именно будет устроено разделение пространства на области КД-деревом для эффективного его использования.
Возьмем проекции всех параллелепипедов на плоскость, перпендикулярную оси <tex>OZ</tex> - множество прямоугольников <tex>V</tex>. Рассмотрим все <tex>x</tex>-координаты границ прямоугольников и разобьем ось <tex>OX</tex> на <tex>2\sqrt{n}</tex> промежутков так, чтобы в каждом из них было не более <tex>\sqrt{n}</tex> границ. Получившиеся <tex>2\sqrt{n}</tex> полос разобьем горизонтальными линиями на итоговые области по отдельности следующим образом. Для каждой полосы у всех прямоугольников, имеющих вертикальную границу в этой полосе, возьмем все горизонтальные границы (их будет не более <tex>2\sqrt{n}</tex>), а среди всех горизонтальных границ, пересекающих область, возьмем каждую <tex>\sqrt{n}</tex>-ую. Таким образом, в итоге получим не более <tex>4\sqrt{n}</tex> областей в каждой вертикальной полосе.

Полученные области обладают следующими свойствами:
# Всего <tex>O(n)</tex> областей.
# Каждый прямоугольник частично покрывает не более <tex>O(\sqrt{n})</tex> областей.
# Никакая вершина прямоугольника не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти утверждения несложно следуют из алгоритма разбиения на области.

Каждая из этих областей будет одним из листов КД-дерева. Для построения всего дерева сначала объединим все области для одной вертикальной полосы, сливая соседние области в одну вершину дерева. После этого аналогично будем сливать соседние вертикальные области и в итоге получим корневую вершину, представляющую всю плоскость.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n)</tex> вершин.
# Каждый прямоугольник хранится в <tex>O(\sqrt{n})</tex> листьях.
# Каждый прямоугольник влияет на <tex>O(\sqrt{n} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина прямоугольника не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти свойства следуют из свойств разбиения на области и логарифмической высоты дерева.

Наконец, получаем итоговый алгоритм. Идем сканирующей плоскостью, перпендикулярной оси <tex>OZ</tex>, добавляя или удаляя прямоугольники из КД-дерева. При этом каждый раз пересчитываем площадь объединения всех текущих прямоугольников и, складывая эти величины, умноженные на соответствующую разницу <tex>z</tex>-координат, получаем общий объем.

Каждая операция добавления или удаления прямоугольников требует изменения <tex>O(\sqrt{n} \log n)</tex> полей <tex>Cover</tex>, а также операций добавления или удаления прямоугольника в <tex>O(\sqrt{n})</tex> листьев, каждая из которых требует <tex>O(\log n)</tex> времени. Итого, на одну операцию добавления или удаления требуется <tex>O(\sqrt{n} \log n)</tex> времени, суммарно таких операций <tex>O(n)</tex>, и получается, что весь алгоритм работает за <tex>O(n \sqrt{n} \log n)</tex>. При этом на хранение всей информации требуется <tex>O(n\sqrt{n}) </tex> памяти - каждый лист может хранить <tex>O(\sqrt{n})</tex> прямоугольников.

Все идеи этого алгоритма обобщаются с трехмерного пространства на <tex>d</tex>-мерное следующим образом.

Построение разбиения на области в <tex>d</tex>-мерном пространстве вместо двухмерного ведется сначала аналогичным образом, только после того, как в двухмерном пространстве получаются итоговые прямоугольные области, здесь продолжаются аналогичные разбиения вдоль третьей координаты - берутся все границы третьей координаты гиперпараллелепипедов, имеющие границу по первой или второй координате строго внутри области (их получается не более <tex>O(\sqrt{n})</tex>), а также каждую <tex>O(\sqrt{n})</tex>-ую границу по третьей координате гиперепараллелепипедов, пересекающих эту область, и так далее.

Это разбиение на области имеет следующие аналогичные свойства:
# Всего <tex>O(n^{d/2})</tex> областей.
# Каждый гиперпараллелепипед частично покрывает не более <tex>O(n^{(d-1)/2})</tex> областей.
# Никакая вершина гиперпараллелепипеда не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Общее дерево строится аналогичным образом и оно в итоге состоит из <tex>d</tex> уровней, каждый из которых имеет высоту <tex>O(\log n)</tex> в дереве - самый верхний уровень содержит разбиение по первой координате, следующий - по второй и так далее.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n^{d/2})</tex> вершин.
# Каждый гиперпараллелепипед хранится в <tex>O(n^{(d-1)/2})</tex> листьях.
# Каждый гиперпараллелепипед влияет на <tex>O(n^{(d-1)/2} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина гиперпараллелепипеда не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Операции в этом дереве выполняются аналогично трехмерному случаю и каждая операция добавления или удаления в КД-дерево размерности <tex>d</tex> занимает <tex>O(n^{(d-1)/2} \log n)</tex> времени. Выполняя сканирование гиперплоскостью вдоль одной из координат получаем КД-дерево размерности <tex>d-1</tex>, в котором необходимо провести <tex>O(n)</tex> операций, то есть суммарное время работы алгоритма равно: <tex>O(n \times n^{(d-1-1)/2} \log n) = O(n^{d/2} \log n)</tex>. При этом хранение КД-дерева требует <tex>O(n^{d/2})</tex> памяти.

== Источники ==
<references/>

Алгоритмы точного вычисления гиперобъема

2012-06-20T11:48:56Z

Joshik:

{{В разработке}}

== Постановка задачи ==
<tex>x = (x_1, x_2, ..., x_d; x_i \ge 0) \in \mathbb{R}^d</tex> - точка в <tex>d</tex>-мерном пространстве.

Точка <tex>x</tex> доминирует точку <tex>y</tex> (<tex>x \succ y</tex>), если <tex>\forall i : x_i \ge y_i, \exists j : x_j > y_j</tex>.

<tex>X = \{x^1, x^2, ..., x^n\} \subset \mathbb{R}^d</tex> - множество из <tex>n</tex> точек в <tex>d</tex>-мерном пространстве таких, что <tex>\nexists i \neq j : x_i \succ x_j</tex> - никакая точка не доминируется другой точкой из этого множества.

<tex>S(X) = \mu (\bigcup \limits_{x \in X} \{y | x \succ y\}) </tex> - гиперобъем множества <tex>X</tex>, где <tex>\mu</tex> - мера Лебега.

В частности, если <tex>X = \{x\}</tex>, то <tex>S(X) = \prod \limits_{i=1}^{d} x_i</tex>.

Задача: найти точное значение гиперобъема <tex>S(X)</tex> множества из <tex>n</tex> точек <tex>d</tex>-мерного пространоства.

== Алгоритм включения-исключения (Inclusion-Exclusion Algorithm, IEA) ==

Самый простой алгоритм нахождения гиперобъема базируется на идее комбинаторной [[формула включения-исключения|формулы включения-искючения]].
Все множество <tex>X</tex> представляется в виде объединения <tex>n</tex> гиперкубов (<tex>X^i</tex>), соответствующих отдельным точкам <tex>x^i</tex>.

После этого объем всего множества вычисляется по формуле: <center> <tex> S(X) = \sum \limits_{I \in 2^n} (-1)^{|I|+1} S(\bigcap \limits_{ j \in I} X^j) </tex> </center>

Объем пересечения гиперкубов легко считается как произведение по каждой координате минимального значения этой координаты среди всех точек, которым соответствуют гиперкубы.

Таким образом, в этом алгоритме перебираются все подмножества точек множества <tex>X</tex>, для каждого множества находится гиперобъем пересечения соответствующих гиперкубов, и он прибавляется с соответствующим знаком к ответу.
Время работы этого алгоритма составляет <tex>O(n 2^n)</tex>.

== Алгоритм LebMeasure ==

Алгоритм LebMeasure<ref>K. Bringmann, T. Friedrich.: An Efficient Algorithm for Computing Hypervolume Contributions. Evolutionary Computation, Vol. 18, No. 3, pp 383-402, MIT Press. (2010)</ref> обрабатывает точки множества <tex>X</tex> по очереди. Для каждой очередной точки <tex>x</tex> находится объем некоторого максимального по включению гиперкуба, доминируемого эксклюзивно только этой точкой <tex>x</tex> и она заменяется на некоторое множество ''порожденных'' точек, которые доминируют оставшуюся область, доминировшуюся до этого точкой <tex>x</tex>.

Например, если изначально было четыре точки в трехмерном пространстве <tex>(6, 7, 4), (9, 5, 5), (1, 9, 3), (4, 1, 9)</tex>, то точка <tex>(6, 7, 4)</tex> эксклюзивно доминирует куб с одним концов в <tex>(6, 7, 4)</tex>, а другим - в <tex>(4, 5, 3)</tex>. После добавления объема этого куба к ответу, точка <tex>(6, 7, 4)</tex> порождает три точки: <tex>(4, 7, 4), (6, 5, 3), (6, 7, 3)</tex>. При этом точка <tex>(6, 5, 4)</tex> доминируется точкой <tex>(9, 5, 5)</tex> и сразу удаляется из множества <tex>X</tex>.

Обработка точек продолжается, пока все точки не будут обработаны. Таким образом, время работы алгоритма напрямую зависит от количества порожденных точек. Легко заметить, что таких точек не более, чем <tex>n^d</tex>, поскольку каждая координата каждой порожденной точки равна соответствующей координате некоторой точки исходного множества <tex>X</tex>.

К сожалению, эта верхняя оценка является достижимой. Например, если исходное множество <tex>X</tex> имеет вид: <tex>(1, n, n, ..., n), (2, n - 1, n - 1, ..., n - 1), ..., (n, 1, 1, ..., 1)</tex>, и точки обрабатываются именно в этом порядке, то всего будет обработано <tex>n^{d-1}</tex> точка, что показано в [1]. Правда, если в этом примере точки обрабатываются в обратном порядке, то суммарное количество обработанных точек линейно зависит от <tex>n</tex> и <tex>d</tex>. Тем не менее, существуют примеры, для которых любой порядок обработки приводит к экспоненциальной зависимости числа порожденных точек от размерности пространства <tex>d</tex> и эта зависимость близка к <tex>n^d</tex> <ref name="while">
While, L., Hingston, P., Barone, L., Huband, S.: A Faster Algorithm for Calculating Hypervolume. IEEE Transaction on Evolutionary Computation, Vol.10, No. 1, pp 29–38 (2006)</ref>.

== Алгоритм Hypervolume by Slicing Objectives (HSO) ==

Под ''Objectives'' в названии данного алгоритма имеются в виду координаты пространства <tex>\mathbb{R}^d</tex>.

Если алгоритм LebMeasure по очереди рассматривает все точки, то алгоритм HSO<ref>Zhou X., Sun C., Mao N., Li W.: Generalization of HSO algortihms for computing hypervolume for mutiobjective optimization problems, IEEE Congress on Evolutionary Computation (2007)</ref> рассматривает по очереди все координаты, сводя задачу к меньшей на единицу размерности.

Изначально все точки сортируются по первой координате. Значения этой координаты используются для ''расслоения'' (разрезания) всего множества на <tex>n</tex> частей, внутри каждой из которых при движении вдоль первой координаты форма ''разреза'' (гиперплоскости), проходящего перпендикулярно оси первой координаты, не меняется. Таким образом, для подсчета объема каждой части необходимо найти объем разреза, и умножить на длину части вдоль первой координаты. При этом получившийся разрез имеет на единицу меньшую размерность.
Заметим, что после сортировки и расслоения первая часть содержит все исходные <tex>n</tex> точек множества <tex>X</tex>, вторая - все, помимо точки с минимальной первой координатой, вдоль которой происходит расслоение, и т.д., а последняя часть содержит только одну точку с максимальной этой координатой.

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

Описанный процесс можно реализовать как в виде рекурсивной процедуры, расслаивающей множество вдоль очередной координаты и вызывающей себя рекурсивно для каждой части и следующей координаты, так и нерекурсивно, если поддерживать множество всех текущих частей и на очередной итерации разбивать их все вдоль очередной координаты.

Время работы алгоритма напрямую зависит от суммарного количества частей, на которые будет разбито исходное множество. По аналогичным предыдущему алгоритму рассуждениям легко показывается, что частей не более <tex>n^d</tex>. Тем не менее, существует более точная оценка.

{{Теорема
|statement=Суммарное количество частей, полученных алгоритмом HSO, не превышает <center><tex>f(n, d) = {n + d - 2 \choose d - 1}</tex>.</center>

|proof=
База:<center><tex>f(n, 1) = 1</tex>.</center>

Также, из алгоритма расслоения следует, что из части с <tex>n</tex> точками будет получено по одной части для каждого <tex>i = 1...n</tex>, то есть: <center><tex>f(n, d) = \sum \limits_{i = 1}^{n}f(i, d - 1)</tex>.</center>

Докажем полезные леммы:

{{Лемма
|statement=<center><tex>{x \choose y} = {x - r \choose y} + \sum \limits_{i = 1}^{r} {x - i \choose y - 1}</tex>.</center>
|proof=Индукция по <tex>r</tex>. База: <center><tex>{x \choose y} = {x - 1 \choose y} + {x - 1 \choose y - 1}</tex>.</center>
Переход: предполагаем, что утверждение верно для <tex>r = j</tex>. Раскрываем первое слагаемое в правой части предположения аналогично базе и получаем:
<center><tex>{x \choose y} = {x - j - 1 \choose y} + {x - (j + 1) \choose y - 1} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1} = {x - (j + 1) \choose y} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1}</tex>.</center>
}}

{{Лемма
|statement=<center><tex>{x \choose y} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
|proof=Подставляем в предыдущую лемму <tex>r = x - y</tex> и получаем:
<center><tex>{x \choose y} = {y \choose y} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = {x - (x - y + 1) \choose y - 1} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
}}

продолжаем доказательство теоремы по индукции:
<center><tex>f(n, d) = {n + d - 2 \choose d - 1} = \sum \limits_{i = 1}^{n} {n + d - 2 - i \choose d - 2} = \sum \limits_{k = n}^{1} {n + d + 2 - (n - k + 1) \choose d - 2} = \sum \limits_{k = 1}^{n} {k + d - 3 \choose d - 2} = \sum \limits_{k = 1}^{n}{k + (d - 1) - 2 \choose (d - 1) - 1} = \sum \limits_{k = 1}^{n}{f(k, d - 1)}</tex></center>

}}

Таким образом, на этот алгоритм получена более низкая оценка времени работы, чем на предыдущий, и это подтверждается на практике <ref name ="while"/>.

== Сведение к задаче KMP (Klee's Measure Problem) ==

Задача KMP состоит в нахождении объема объединения прямоугольных гиперпараллелепипедов в <tex>d</tex>-мерном пространстве. Как показано в описании алгоритма IEA, исходная задача легко сводится к этой, если каждой точке поставить в соответствие гиперкуб с одной вершиной в центре координат, а противоположной - в этой точке.

Существует несколько алгоритмов решения задачи KMP<ref name = "overmars">Overmars, M.H., Yap, C.K.: New upper bounds in Klee’s measure problem. SIAM Journal on Computing 20(6), pp 1034–1045. (1991)</ref><ref>N. Beume.: S-Metric calculation by considering dominated hypervolume as Klee’s measure problem. Evolutionary Computation, 17(4), pp 477–492. (2009)</ref><ref>N. Beume, G. Rudolph.: Faster S-metric calculation by considering dominated hypervolume as Klee’s measure problem. In Proc. Second International Conference on Computational Intelligence (IASTED ’06), pp 233–238. (2006)
</ref>, самый оптимальный из которых использует идеи сканирующей гиперплоскости, [[Статистики на отрезках. Корневая эвристика|корневой эвристики]] и КД-дерево, и позволяет решать задачу за время <tex>O(n^{d/2}\log n)</tex>. Рассмотрим его.

Для начала изложим идею сканирующей гиперплоскости. Рассмотрим всевозможные значения <tex>d</tex>-ой координаты у точек множества <tex>X</tex>.
Возьмем гиперплоскость, перпендикулярную оси <tex>d</tex>-ой координаты, будем ее двигать по этим значениям от минимального до максимального и рассматривать проекцию всех гиперкубов на эту гиперплоскость. Заметим, что между точками проекция меняться не будет, а в каждой точке будут появляться или исчезать проекции некоторых гиперкубов. Поэтому, если уметь эффективно пересчитывать объем проекции при каждом появлении/исчезновении и прибавлять это значение, умноженное на расстояние между последовательными значениями <tex>d</tex>-ой координаты, то можно легко посчитать суммарный гиперобъем. Как раз для эффективного пересчета используется такая структура, как КД-дерево.

Теперь изложим идею алгоритма на примере трехмерного пространства.
В трехмерном пространстве сканирующая гиперплоскость превращается в обычную плоскость. Будем считать, что она всегда располагается перепендикулярно оси <tex>OZ</tex> и двигается от минимальных значений по этой оси до максимальных. В каждый момент проекция всех параллелепипедов на эту плоскость представляет собой множество прямоугольников и необходимо уметь обрабатывать следующие три запроса:
# добавить прямоугольник
# удалить прямоугольник
# посчитать площадь объединения всех прямоугольников

{{Определение
|definition=КД-деревом называется [[Дерево поиска, наивная реализация|сбалансированное двоичное дерево]], где каждой вершине <tex>\alpha</tex> соответствует некоторая область <tex>C_{\alpha}</tex> <tex>d</tex>-мерного пространства, удовлетворяющая следующим свойствам:
# Корневой вершине <tex>root</tex> соответствует <tex>C_{root}</tex> - все пространство.
# Каждое <tex>C_{\alpha}</tex> является гиперпараллелепипедом (возможно бесконечным в некоторых направлениях).
# Для каждой вершины <tex>\alpha</tex> с детьми <tex>\beta</tex> и <tex>\gamma</tex> верно, что области <tex>C_{\beta}</tex> и <tex>C_{\gamma}</tex> не пересекаются и <tex>C_{\alpha} = C_{\beta} \cup C_{\gamma}</tex>.
}}

Будем хранить прямоугольники в КД-дереве следующим образом:
# Для каждого листа <tex>\alpha</tex> будем хранить все прямоугольники, которые пересекают внутреннюю часть области <tex>C_{\alpha}</tex>.
# Для каждой внутренней вершины <tex>\alpha</tex> будем считать число <tex>Cover_{\alpha}</tex> - количество прямоугольников, которые полностью покрывают область <tex>C_{\alpha}</tex>, но не полностью покрывают область родителя этой вершины <tex>C_{father(\alpha)}</tex>.

Также для каждой вершины будем хранить значение площади пересечения области этой вершины со всеми прямоугольниками <tex>M</tex>, которое определяется следующим образом:
# Для каждого листа <tex>M</tex> равно площади пересечения всех прямоугольников с областью этого листа.
# Для каждой внутренней вершины <tex>\alpha</tex> если <tex>Cover_{\alpha} > 0</tex>, тогда <tex>M</tex> равно площади всей области <tex>C_{\alpha}</tex>, иначе <tex>M</tex> равно сумме этих значений для двух сыновей этой вершины.

Таким образом, значение для корневой вершины <tex>M_{root}</tex> и является искомой площадью объединения всех прямоугольников.

Рассмотрим операцию добавления прямоугольника (box) в КД-дерево.

procedure Insert(box, <tex>\alpha</tex>)
if <tex>\alpha</tex> - лист then
добавить box в эту вершину
пересчитать <tex>M_{\alpha}</tex>
elseif box покрывает область <tex>C_{\alpha}</tex> then
<tex>Cover_{\alpha} = Cover_{\alpha} + 1</tex>
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
elseif box пересекает область <tex>C_{\alpha}</tex> then
Insert(box, leftson(<tex>\alpha</tex>))
Insert(box, rightson(<tex>\alpha</tex>))
if <tex>Cover_{\alpha} > 0 </tex> then
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
else
<tex>M_{\alpha} = M_{leftson(\alpha)} + M_{rightson(\alpha)}</tex>

Аналогично устроена операция удаления прямоугольника - она выполняет обратные действия.

Рассмотрим теперь как именно будет устроено разделение пространства на области КД-деревом для эффективного его использования.
Возьмем проекции всех параллелепипедов на плоскость, перепендикулярную оси <tex>OZ</tex> - множество прямоугольников <tex>V</tex>. Рассмотрим все <tex>x</tex>-координаты границ прямоугольников и разобьем ось <tex>OX</tex> на <tex>2\sqrt{n}</tex> промежутков так, чтобы в каждом из них было не более <tex>\sqrt{n}</tex> границ. Получившиеся <tex>2\sqrt{n}</tex> полос разобьем горизонтальными линиями на итоговые области по отдельности следующим образом. Для каждой полосы у всех прямоугольников, имеющих вертикальную границу в этой полосе, возьмем все горизонтальные границы (их будет не более <tex>2\sqrt{n}</tex>), а среди всех горизонтальных границ, пересекающих область, возьмем каждую <tex>\sqrt{n}</tex>-ую. Таким образом, в итоге получим не более <tex>4\sqrt{n}</tex> областей в каждой верикальной полосе.

Полученные области обладают следующими свойствами:
# Всего <tex>O(n)</tex> областей.
# Каждый прямоугольник частично покрывает не более <tex>O(\sqrt{n})</tex> областей.
# Никакая вершина прямоугольника не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти утверждения несложно следуют из алгоритма разбиения на области.

Каждая из этих областей будет одним из листов КД-дерева. Для построения всего дерева сначала объединим все области для одной вертикальной полосы, сливая соседние области в одну вершину дерева. После этого аналогично будем сливать соседние вертикальные области и в итоге получим корневую вершину, представляющую всю плоскость.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n)</tex> вершин.
# Каждый прямоугольник хранится в <tex>O(\sqrt{n})</tex> листьях.
# Каждый прямоугольник влияет на <tex>O(\sqrt{n} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина прямоугольника не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти свойства следуют из свойств разбиения на области и логарифмической высоты дерева.

Наконец, получаем итоговый алгоритм. Идем сканирующей плоскостью, перепендикулярной оси <tex>OZ</tex>, добавляя или удаляя прямоугольники из КД-дерева. При этом каждый раз пересчитываем площадь объединения всех текущих прямоугольников и, складывая эти величины, умноженные на соответствующую разницу <tex>z</tex>-координат, получаем общий объем.

Каждая операция добавления или удаления прямоугольников требует изменения <tex>O(\sqrt{n} \log n)</tex> полей <tex>Cover</tex>, а также операций добавления или удаления прямоугольника в <tex>O(\sqrt{n})</tex> листьев, каждая из которых требует <tex>O(\log n)</tex> времени. Итого, на одну операцию добавления или удаления требуется <tex>O(\sqrt{n} \log n)</tex> времени, суммарно таких операций <tex>O(n)</tex>, и получается, что весь алгоритм работает за <tex>O(n \sqrt{n} \log n)</tex>. При этом на хранение всей информации требуется <tex>O(n\sqrt{n}) </tex> памяти - каждый лист может хранить <tex>O(\sqrt{n})</tex> прямоугольников.

Все идеи этого алгоритма обобщаются с трехмерного пространства на <tex>d</tex>-мерное следующим образом.

Построение разбиения на области в <tex>d</tex>-мерном пространстве вместо двухмерного ведется сначала аналогичным образом, только после того, как в двухмерном пространстве получаются итоговые прямоугольные области, здесь продолжаются аналогичные разбиения вдоль третьей координаты - берутся все границы третьей координаты гиперпараллелепипедов, имеющие границу по первой или второй координате строго внутри области (их получается не более <tex>O(\sqrt{n})</tex>), а также каждую <tex>O(\sqrt{n})</tex>-ую границу по третьей координате гиперепараллелепипедов, пересекающих эту область, и так далее.

Это разбиение на области имеет следующие аналогичные свойства:
# Всего <tex>O(n^{d/2})</tex> областей.
# Каждый гиперпараллелепипед частично покрывает не более <tex>O(n^{(d-1)/2})</tex> областей.
# Никакая вершина гиперпараллелепипеда не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Общее дерево строится аналогичным образом и оно в итоге состоит из <tex>d</tex> уровней, каждый из которых имеет высоту <tex>O(\log n)</tex> в дереве - самый верхний уровень содержит разбиение по первой координате, следующий - по второй и так далее.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n^{d/2})</tex> вершин.
# Каждый гиперпараллелепипед хранится в <tex>O(n^{(d-1)/2})</tex> листьях.
# Каждый гиперпараллелепипед влияет на <tex>O(n^{(d-1)/2} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина гиперпараллелепипеда не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Операции в этом дереве выполняются аналогично трехмерному случаю и каждая операция добавления или удаления в КД-дерево размерности <tex>d</tex> занимает <tex>O(n^{(d-1)/2} \log n)</tex> времени. Выполняя сканирование гиперплоскостью вдоль одной из координат получаем КД-дерево размерности <tex>d-1</tex>, в котором необходимо провести <tex>O(n)</tex> операций, то есть суммарное время работы алгоритма равно: <tex>O(n \times n^{(d-1-1)/2} \log n) = O(n^{d/2} \log n)</tex>. При этом хранение КД-дерева требует <tex>O(n^{d/2})</tex> памяти.

== Источники ==
<references/>

Алгоритмы точного вычисления гиперобъема

2012-06-18T06:43:17Z

Joshik: орфографические ошибки

{{В разработке}}

== Постановка задачи ==
<tex>x = (x_1, x_2, ..., x_d; x_i \ge 0) \in R^d</tex> - точка в <tex>d</tex>-мерном пространстве.

Точка <tex>x</tex> доминирует точку <tex>y</tex> (<tex>x \succ y</tex>), если <tex>\forall i : x_i \ge y_i, \exists j : x_j > y_j</tex>.

<tex>X = (x^1, x^2, ..., x^n) \subset R^d</tex> - множество из <tex>n</tex> точек в <tex>d</tex>-мерном пространстве таких, что <tex>\nexists i \neq j : x_i \succ x_j</tex> - никакая точка не доминируется другой точкой из этого множества.

<tex>S(X) = \mu (\bigcup \limits_{x \in X} \{y | x \succ y\}) </tex> - гиперобъем множества <tex>X</tex>.

В частности, если <tex>X = \{x\}</tex>, то <tex>S(X) = \prod \limits_{i=1}^{d} x_i</tex>.

Задача: найти точное значение гиперобъема <tex>S(X)</tex> множества из <tex>n</tex> точек <tex>d</tex>-мерного пространоства.

== Алгоритм включения-исключения (Inclusion-Exclusion Algorithm, IEA) ==

Самый простой алгоритм нахождения гиперобъема базируется на идее комбинаторной [[формула включения-исключения|формулы включения-искючения]].
Все множество <tex>X</tex> представляется в виде объединения <tex>n</tex> гиперкубов (<tex>X^i</tex>), соответствующих отдельным точкам <tex>x^i</tex>.

После этого объем всего множества вычисляется по формуле: <center> <tex> S(X) = \sum \limits_{I \in 2^n} (-1)^{|I|+1} S(\bigcap \limits_{ j \in I} X^j) </tex> </center>

Объем пересечения гиперкубов легко считается как произведение по каждой координате минимального значения этой координаты среди всех точек, которым соответствуют гиперкубы.

Таким образом, в этом алгоритме перебираются все подмножества точек множества <tex>X</tex>, для каждого множества находится гиперобъем пересечения соответствующих гиперкубов, и он прибавляется с соответствующим знаком к ответу.
Время работы этого алгоритма составляет <tex>O(n 2^n)</tex>.

== Алгоритм LebMeasure ==

Алгоритм LebMeasure обрабатывает точки множества <tex>X</tex> по очереди. Для каждой очередной точки <tex>x</tex> находится объем некоторого максимального по включению гиперкуба, доминируемого эксклюзивно только этой точкой <tex>x</tex> и она заменяется на некоторое множество ''порожденных'' точек, которые доминируют оставшуюся область, доминировшуюся до этого точкой <tex>x</tex>.

Например, если изначально было четыре точки в трехмерном пространстве <tex>(6, 7, 4), (9, 5, 5), (1, 9, 3), (4, 1, 9)</tex>, то точка <tex>(6, 7, 4)</tex> эксклюзивно доминирует куб с одним концов в <tex>(6, 7, 4)</tex>, а другим - в <tex>(4, 5, 3)</tex>. После добавления объема этого куба к ответу, точка <tex>(6, 7, 4)</tex> порождает три точки: <tex>(4, 7, 4), (6, 5, 3), (6, 7, 3)</tex>. При этом точка <tex>(6, 5, 4)</tex> доминируется точкой <tex>(9, 5, 5)</tex> и сразу удаляется из множества <tex>X</tex>.

Обработка точек продолжается, пока все точки не будут обработаны. Таким образом, время работы алгоритма напрямую зависит от количества порожденных точек. Легко заметить, что таких точек не более, чем <tex>n^d</tex>, поскольку каждая координата каждой порожденной точки равна соответствующей координате некоторой точки исходного множества <tex>X</tex>.

К сожалению, эта верхняя оценка является достижимой. Например, если исходное множество <tex>X</tex> имеет вид: <tex>(1, n, n, ..., n), (2, n - 1, n - 1, ..., n - 1), ..., (n, 1, 1, ..., 1)</tex>, и точки обрабатываются именно в этом порядке, то всего будет обработано <tex>n^{d-1}</tex> точка, что показано в [1]. Правда, если в этом примере точки обрабатываются в обратном порядке, то суммарное количество обработанных точек линейно зависит от <tex>n</tex> и <tex>d</tex>. Тем не менее, существуют примеры, для которых любой порядок обработки приводит к экспоненциальной зависимости числа порожденных точек от размерности пространства <tex>d</tex> и эта зависимость близка к <tex>n^d</tex> [1].

== Алгоритм Hypervolume by Slicing Objectives (HSO) ==

Под ''Objectives'' в названии данного алгоритма имеются в виду координаты пространства <tex>R^d</tex>.

Если алгоритм LebMeasure по очереди рассматривает все точки, то алгоритм HSO рассматривает по очереди все координаты, сводя задачу к меньшей на единицу размерности.

Изначально все точки сортируются по первой координате. Значения этой координаты используются для ''расслоения'' (разрезания) всего множества на <tex>n</tex> частей, внутри каждой из которых при движении вдоль первой координаты форма ''разреза'' (гиперплоскости), проходящего перпендикулярно оси первой координаты, не меняется. Таким образом, для подсчета объема каждой части необходимо найти объем разреза, и умножить на длину части вдоль первой координаты. При этом получившийся разрез имеет на единицу меньшую размерность.
Заметим, что после сортировки и расслоения первая часть содержит все исходные <tex>n</tex> точек множества <tex>X</tex>, вторая - все, помимо точки с минимальной первой координатой, вдоль которой происходит расслоение, и т.д., а последняя часть содержит только одну точку с максимальной этой координатой.

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

Описанный процесс можно реализовать как в виде рекурсивной процедуры, расслаивающей множество вдоль очередной координаты и вызывающей себя рекурсивно для каждой части и следующей координаты, так и нерекурсивно, если поддерживать множество всех текущих частей и на очередной итерации разбивать их все вдоль очередной координаты.

Время работы алгоритма напрямую зависит от суммарного количества частей, на которые будет разбито исходное множество. По аналогичным предыдущему алгоритму рассуждениям легко показывается, что частей не более <tex>n^d</tex>. Тем не менее, существует более точная оценка.

{{Теорема
|statement=Суммарное количество частей, полученных алгоритмом HSO, не превышает <center><tex>f(n, d) = {n + d - 2 \choose d - 1}</tex>.</center>

|proof=
База:<center><tex>f(n, 1) = 1</tex>.</center>

Также, из алгоритма расслоения следует, что из части с <tex>n</tex> точками будет получено по одной части для каждого <tex>i = 1...n</tex>, то есть: <center><tex>f(n, d) = \sum \limits_{i = 1}^{n}f(i, d - 1)</tex>.</center>

Докажем полезные леммы:

{{Лемма
|statement=<center><tex>{x \choose y} = {x - r \choose y} + \sum \limits_{i = 1}^{r} {x - i \choose y - 1}</tex>.</center>
|proof=Индукция по <tex>r</tex>. База: <center><tex>{x \choose y} = {x - 1 \choose y} + {x - 1 \choose y - 1}</tex>.</center>
Переход: предполагаем, что утверждение верно для <tex>r = j</tex>. Раскрываем первое слагаемое в правой части предположения аналогично базе и получаем:
<center><tex>{x \choose y} = {x - j - 1 \choose y} + {x - (j + 1) \choose y - 1} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1} = {x - (j + 1) \choose y} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1}</tex>.</center>
}}

{{Лемма
|statement=<center><tex>{x \choose y} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
|proof=Подставляем в предыдущую лемму <tex>r = x - y</tex> и получаем:
<center><tex>{x \choose y} = {y \choose y} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = {x - (x - y + 1) \choose y - 1} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
}}

продолжаем доказательство теоремы по индукции:
<center><tex>f(n, d) = {n + d - 2 \choose d - 1} = \sum \limits_{i = 1}^{n} {n + d - 2 + i \choose d - 2} = \sum \limits_{k = n}^{1} {n + d + 2 - (n - k + 1) \choose d - 2} = \sum \limits_{k = 1}^{n} {k + d - 3 \choose d - 2} = \sum \limits_{k = 1}^{n}{k + (d - 1) - 2 \choose (d - 1) - 1} = \sum \limits_{k = 1}^{n}{f(k, d - 1)}</tex></center>

}}

Таким образом, на этот алгоритм получена более низкая оценка времени работы, чем на предыдущий, и это подтверждается на практике [1].

== Сведение к задаче KMP (Klee's Measure Problem) ==

Задача KMP состоит в нахождении объема объединения прямоугольных гиперпараллелепипедов в <tex>d</tex>-мерном пространстве. Как показано в описании алгоритма IEA, исходная задача легко сводится к этой, если каждой точке поставить в соответствие гиперкуб с одной вершиной в центре координат, а противоположной - в этой точке.

Существует несколько алгоритмов решения задачи KMP, самый оптимальный из которых использует идеи сканирующей гиперплоскости, [[Статистики на отрезках. Корневая эвристика|корневой эвристики]] и КД-дерево, и позволяет решать задачу за время <tex>O(n^{d/2}\log n)</tex>. Рассмотрим его.

Для начала изложим идею сканирующей гиперплоскости. Рассмотрим всевозможные значения <tex>d</tex>-ой координаты у точек множества <tex>X</tex>.
Возьмем гиперплоскость, перпендикулярную оси <tex>d</tex>-ой координаты, будем ее двигать по этим значениям от минимального до максимального и рассматривать проекцию всех гиперкубов на эту гиперплоскость. Заметим, что между точками проекция меняться не будет, а в каждой точке будут появляться или исчезать проекции некоторых гиперкубов. Поэтому, если уметь эффективно пересчитывать объем проекции при каждом появлении/исчезновении и прибавлять это значение, умноженное на расстояние между последовательными значениями <tex>d</tex>-ой координаты, то можно легко посчитать суммарный гиперобъем. Как раз для эффективного пересчета используется такая структура, как КД-дерево.

Теперь изложим идею алгоритма на примере трехмерного пространства.
Будем считать, что сканирующая плоскость двигается перепендикулярно оси <tex>OZ</tex>. В каждый момент проекция всех параллелепипедов на эту плоскость предствляет собой множество прямоугольников и необходимо уметь обрабатывать следующие три запроса:
# добавить прямоугольник
# удалить прямоугольник
# посчитать площадь объединения всех прямоугольников

{{Определение
|definition=КД-деревом называется [[Дерево поиска, наивная реализация|сбалансированное двоичное дерево]], где каждой вершине <tex>\alpha</tex> соответствует некоторая область <tex>C_{\alpha}</tex> <tex>d</tex>-мерного пространства, удовлетворяющая следующим свойствам:
# Корневой вершине <tex>root</tex> соответствует <tex>C_{root}</tex> - все пространство.
# Каждое <tex>C_{\alpha}</tex> является гиперпараллелепипедом (возможно бесконечным в некоторых направлениях).
# Для каждой вершины <tex>\alpha</tex> с детьми <tex>\beta</tex> и <tex>\gamma</tex> верно, что области <tex>C_{\beta}</tex> и <tex>C_{\gamma}</tex> не пересекаются и <tex>C_{\alpha} = C_{\beta} \cup C_{\gamma}</tex>.
}}

Будем хранить прямоугольники в КД-дереве следующим образом:
# Для каждого листа <tex>\alpha</tex> будем хранить все прямоугольники, которые пересекают внутреннюю часть области <tex>C_{\alpha}</tex>.
# Для каждой внутренней вершины <tex>\alpha</tex> будем считать число <tex>Cover_{\alpha}</tex> - количество прямоугольников, которые полностью покрывают область <tex>C_{\alpha}</tex>, но не полностью покрывают область родителя этой вершины <tex>C_{father(\alpha)}</tex>.

Также для каждой вершины будем хранить значение площади пересечения области этой вершины со всеми прямоугольниками <tex>M</tex>, которое определяется следующим образом:
# Для каждого листа <tex>M</tex> равно площади пересечения всех прямоугольников с областью этого листа.
# Для каждой внутренней вершины <tex>\alpha</tex> если <tex>Cover_{\alpha} > 0</tex>, тогда <tex>M</tex> равно площади всей области <tex>C_{\alpha}</tex>, иначе <tex>M</tex> равно сумме этих значений для двух сыновей этой вершины.

Таким образом, значение для корневой вершины <tex>M_{root}</tex> и является искомой площадью объединения всех прямоугольников.

Рассмотрим операцию добавления прямоугольника (box) в КД-дерево.

procedure Insert(box, <tex>\alpha</tex>)
if <tex>\alpha</tex> - лист then
добавить box в эту вершину
пересчитать <tex>M_{\alpha}</tex>
elseif box покрывает область <tex>C_{\alpha}</tex> then
<tex>Cover_{\alpha} = Cover_{\alpha} + 1</tex>
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
elseif box пересекает область <tex>C_{\alpha}</tex> then
Insert(box, leftson(<tex>\alpha</tex>))
Insert(box, rightson(<tex>\alpha</tex>))
if <tex>Cover_{\alpha} > 0 </tex> then
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
else
<tex>M_{\alpha} = M_{leftson(\alpha)} + M_{rightson(\alpha)}</tex>

Аналогично устроена операция удаления прямоугольника - она выполняет обратные действия.

Рассмотрим теперь как именно будет устроено разделение пространства на области КД-деревом для эффективного его использования.
Возьмем проекции всех параллелепипедов на плоскость, перепендикулярную оси <tex>OZ</tex> - множество прямоугольников <tex>V</tex>. Рассмотрим все <tex>x</tex>-координаты границ прямоугольников и разобьем ось <tex>OX</tex> на <tex>2\sqrt{n}</tex> промежутков так, чтобы в каждом из них было не более <tex>\sqrt{n}</tex> границ. Получившиеся <tex>2\sqrt{n}</tex> полос разобьем горизонтальными линиями на итоговые области по отдельности следующим образом. Для каждой полосы у всех прямоугольников, имеющих вертикальную границу в этой полосе, возьмем все горизонтальные границы (их будет не более <tex>2\sqrt{n}</tex>), а среди всех горизонтальных границ, пересекающих область, возьмем каждую <tex>\sqrt{n}</tex>-ую. Таким образом, в итоге получим не более <tex>4\sqrt{n}</tex> областей в каждой верикальной полосе.

Полученные области обладают следующими свойствами:
# Всего <tex>O(n)</tex> областей.
# Каждый прямоугольник частично покрывает не более <tex>O(\sqrt{n})</tex> областей.
# Никакая вершина прямоугольника не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти утверждения несложно следуют из алгоритма разбиения на области.

Каждая из этих областей будет одним из листов КД-дерева. Для построения всего дерева сначала объединим все области для одной вертикальной полосы, сливая соседние области в одну вершину дерева. После этого аналогично будем сливать соседние вертикальные области и в итоге получим корневую вершину, представляющую всю плоскость.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n)</tex> вершин.
# Каждый прямоугольник хранится в <tex>O(\sqrt{n})</tex> листьях.
# Каждый прямоугольник влияет на <tex>O(\sqrt{n} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина прямоугольника не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти свойства следуют из свойств разбиения на области и логарифмической высоты дерева.

Наконец, получаем итоговый алгоритм. Идем сканирующей плоскостью, перепендикулярной оси <tex>OZ</tex>, добавляя или удаляя прямоугольники из КД-дерева. При этом каждый раз пересчитываем площадь объединения всех текущих прямоугольников и, складывая эти величины, умноженные на соответствующую разницу <tex>z</tex>-координат, получаем общий объем.

Каждая операция добавления или удаления прямоугольников требует изменения <tex>O(\sqrt{n} \log n)</tex> полей <tex>Cover</tex>, а также операций добавления или удаления прямоугольника в <tex>O(\sqrt{n})</tex> листьев, каждая из которых требует <tex>O(\log n)</tex> времени. Итого, на одну операцию добавления или удаления требуется <tex>O(\sqrt{n} \log n)</tex> времени, суммарно таких операций <tex>O(n)</tex>, и получается, что весь алгоритм работает за <tex>O(n \sqrt{n} \log n)</tex>. При этом на хранение всей информации требуется <tex>O(n\sqrt{n}) </tex> памяти - каждый лист может хранить <tex>O(\sqrt{n})</tex> прямоугольников.

Все идеи этого алгоритма обобщаются с трехмерного пространства на <tex>d</tex>-мерное следующим образом.

Построение разбиения на области в <tex>d</tex>-мерном пространстве вместо двухмерного ведется сначала аналогичным образом, только после того, как в двухмерном пространстве получаются итоговые прямоугольные области, здесь продолжаются аналогичные разбиения вдоль третьей координаты - берутся все границы третьей координаты гиперпараллелепипедов, имеющие границу по первой или второй координате строго внутри области (их получается не более <tex>O(\sqrt{n})</tex>), а также каждую <tex>O(\sqrt{n})</tex>-ую границу по третьей координате гиперепараллелепипедов, пересекающих эту область, и так далее.

Это разбиение на области имеет следующие аналогичные свойства:
# Всего <tex>O(n^{d/2})</tex> областей.
# Каждый гиперпараллелепипед частично покрывает не более <tex>O(n^{(d-1)/2})</tex> областей.
# Никакая вершина гиперпараллелепипеда не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Общее дерево строится аналогичным образом и оно в итоге состоит из <tex>d</tex> уровней, каждый из которых имеет высоту <tex>O(\log n)</tex> в дереве - самый верхний уровень содержит разбиение по первой координате, следующий - по второй и так далее.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n^{d/2})</tex> вершин.
# Каждый гиперпараллелепипед хранится в <tex>O(n^{(d-1)/2})</tex> листьях.
# Каждый гиперпараллелепипед влияет на <tex>O(n^{(d-1)/2} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина гиперпараллелепипеда не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Операции в этом дереве выполняются аналогично трехмерному случаю и каждая операция добавления или удаления в КД-дерево размерности <tex>d</tex> занимает <tex>O(n^{(d-1)/2} \log n)</tex> времени. Выполняя сканирование гиперплоскостью вдоль одной из координат получаем КД-дерево размерности <tex>d-1</tex>, в котором необходимо провести <tex>O(n)</tex> операций, то есть суммарное время работы алгоритма равно: <tex>O(n \times n^{(d-1-1)/2} \log n) = O(n^{d/2} \log n)</tex>. При этом хранение КД-дерева требует <tex>O(n^{d/2})</tex> памяти.

== Источники ==
# While, L., Hingston, P., Barone, L., Huband, S.: A Faster Algorithm for Calculating Hypervolume. IEEE Transaction on Evolutionary Computation, Vol.10, No. 1, pp 29–38 (2006)
# [http://rain.ifmo.ru/~tsarev/teaching/ea-2012/seminars/1320.pdf Zhou X., Sun C., Mao N., Li W.: Generalization of HSO algortihms for computing hypervolume for mutiobjective optimization problems, IEEE Congress on Evolutionary Computation (2007)]
# K. Bringmann, T. Friedrich.: An Efficient Algorithm for Computing Hypervolume Contributions. Evolutionary Computation, Vol. 18, No. 3, pp 383-402, MIT Press. (2010)
# N. Beume.: S-Metric calculation by considering dominated hypervolume as Klee’s measure problem. Evolutionary Computation, 17(4), pp 477–492. (2009)
# N. Beume, G. Rudolph.: Faster S-metric calculation by considering dominated hypervolume as Klee’s measure problem. In Proc. Second International Conference on Computational Intelligence (IASTED ’06), pp 233–238. (2006)
# Overmars, M.H., Yap, C.K.: New upper bounds in Klee’s measure problem. SIAM Journal on Computing 20(6), pp 1034–1045. (1991)
# Klee, V.: Can the measure of S[ai, bi] be computed in less than O(n log n) steps? In: American Mathematical Monthly. Volume 84, pp 284–285. (1977)

Алгоритмы точного вычисления гиперобъема

2012-06-17T23:20:59Z

Joshik:

{{В разработке}}

== Постановка задачи ==
<tex>x = (x_1, x_2, ..., x_d; x_i \ge 0) \in R^d</tex> - точка в <tex>d</tex>-мерном пространстве.

Точка <tex>x</tex> доминирует точку <tex>y</tex> (<tex>x \succ y</tex>), если <tex>\forall i : x_i \ge y_i, \exists j : x_j > y_j</tex>.

<tex>X = (x^1, x^2, ..., x^n) \subset R^d</tex> - множество из <tex>n</tex> точек в <tex>d</tex>-мерном пространстве таких, что <tex>\nexists i \neq j : x_i \succ x_j</tex> - никакая точка не доминируется другой точкой из этого множества.

<tex>S(X) = \mu (\bigcup \limits_{x \in X} \{y | y \succ x\}) </tex> - гиперобъем множества <tex>X</tex>.

В частности, если <tex>X = \{x\}</tex>, то <tex>S(X) = \prod \limits_{i=1}^{d} x_i</tex>.

Задача: найти точное значение гиперобъема <tex>S(X)</tex> множества из <tex>n</tex> точек <tex>d</tex>-мерного пространоства.

== Алгоритм включения-исключения (Inclusion-Exclusion Algorithm, IEA) ==

Самый простой алгоритм нахождения гиперобъема базируется на идее комбинаторной [[формула включения-исключения|формулы включения-искючения]].
Все множество <tex>X</tex> представляется в виде объединения <tex>n</tex> гиперкубов (<tex>X^i</tex>), соответствующих отдельным точкам <tex>x^i</tex>.

После этого объем всего множества вычисляется по формуле: <center> <tex> S(X) = \sum \limits_{I \in 2^n} (-1)^{|I|+1} S(\bigcap \limits_{ j \in I} X^j) </tex> </center>

Объем пересечения гиперкубов легко считается как произведение по каждой координате минимального значения этой координаты среди всех точек, которым соответствуют гиперкубы.

Таким образом, в этом алгоритме перебираются все подмножества точек множества <tex>X</tex>, для каждого множества находится гиперобъем пересечения соответствующих гиперкубов и он прибавляется с соответствующим знаком к ответу.
Время работы этого алгоритма составляет <tex>O(n 2^n)</tex>.

== Алгоритм LebMeasure ==

Алгоритм LebMeasure обрабатывает точки множества <tex>X</tex> по очереди. Для каждой очередной точки <tex>x</tex> находится объем некоторого максимального по включению гиперкуба, доминируемого эксклюзивно только этой точкой <tex>x</tex> и она заменяется на некоторое множество ''порожденных'' точек, которые доминируют оставшуюся область, доминировшуюся до этого точкой <tex>x</tex>.

Например, если изначально было четыре точки в трехмерном пространстве <tex>(6, 7, 4), (9, 5, 5), (1, 9, 3), (4, 1, 9)</tex>, то точка <tex>(6, 7, 4)</tex> эксклюзивно доминирует куб с одним концов в <tex>(6, 7, 4)</tex>, а другим - в <tex>(4, 5, 3)</tex>. После добавления объема этого куба в ответу, точка <tex>(6, 7, 4)</tex> порождает три точки: <tex>(4, 7, 4), (6, 5, 3), (6, 7, 3)</tex>. При этом точка <tex>(6, 5, 4)</tex> доминируется точкой <tex>(9, 5, 5)</tex> и сразу удаляется из множества <tex>X</tex>.

Обработка точек продолжается, пока все точки не будут обработаны. Таким образом, время работы алгоритмы напрямую зависит от количества порожденных точек. Легко заметить, что таких точек не более, чем <tex>n^d</tex>, поскольку каждая координата каждой порожденной точки равна соответствующей координате некоторой точки исходного множества <tex>X</tex>.

К сожалению, эта верхняя оценка является достижимой. Например, если исходное множество <tex>X</tex> имеет вид: <tex>(1, n, n, ..., n), (2, n - 1, n - 1, ..., n - 1), ..., (n, 1, 1, ..., 1)</tex>, и точки обрабатываются в этом порядке, то всего будет обработано <tex>n^{d-1}</tex> точка, что показано в [1]. Правда, если в этом примере точки обрабатываются в обратном порядке, то суммарное количество обработанных точек линейно зависит от <tex>n</tex> и <tex>d</tex>. Тем не менее, существуют примеры, для которых любой порядок обработки приводит к экспоненциальной зависимости числа порожденных точек от размерности пространства <tex>d</tex> и близок к <tex>n^d</tex> [1].

== Алгоритм Hypervolume by Slicing Objectives (HSO) ==

Под ''Objectives'' в названии данного алгоритма имеются в виду координаты пространства <tex>R^d</tex>.

Если алгоритм LebMeasure по очереди рассматривает все точки, то алгоритм HSO рассматривает по очереди все координаты, сводя задачу к меньшей на единицу размерности.

Изначально все точки сортируются по первой координате. Значения этой координаты используются для ''расслоения''(разрезания) всего множества на <tex>n</tex> частей, внутри каждой из которых при движении вдоль первой координаты форма ''разреза'' перпендикулярно оси первой координаты не меняется. Таким образом, для подсчета объема каждой части необходимо найти объем разреза и умножить длину части вдоль первой координаты. При этом получившийся разрез имеет на единицу меньшую размерность.
Заметим, что после сортировки и расслоения первая часть содержит все <tex>n</tex> точек, вторая - все, помимо точки с минимальной координатой, вдоль которой происходит расслоение и т.д., а последняя часть содержит только одну точку с максимальной этой координатой.

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

Описанный процесс можно реализовать как в виде рекурсивной процедуры, расслаивающей множество вдоль очередной координаты и вызывающей себя рекурсивно для каждой части и следующей координаты, так и нерекурсивно, если поддерживать множество всех текущих частей и на очередной итерации разбивать их все вдоль очередной координаты.

Время работы алгоритма напрямую зависит от суммарного количества частей, на которые будет разбито исходное множество. По аналогичным предыдущему алгоритму рассуждениям легко показывается, что частей не более <tex>n^d</tex>. Тем не менее, существует более точная оценка.

{{Теорема
|statement=Суммарное количество частей, полученных алгоритмом HSO, не превышает <center><tex>f(n, d) = {n + d - 2 \choose d - 1}</tex>.</center>

|proof=
База:<center><tex>f(n, 1) = 1</tex>.</center>

Также, из алгоритма расслоения следует, что из части с <tex>n</tex> точками будет получено по одной части для каждого <tex>i = 1...n</tex>, то есть: <center><tex>f(n, d) = \sum \limits_{i = 1}^{n}f(i, d - 1)</tex>.</center>

Докажем полезные леммы:

{{Лемма
|statement=<center><tex>{x \choose y} = {x - r \choose y} + \sum \limits_{i = 1}^{r} {x - i \choose y - 1}</tex>.</center>
|proof=Индукция по <tex>r</tex>. База: <center><tex>{x \choose y} = {x - 1 \choose y} + {x - 1 \choose y - 1}</tex>.</center>
Переход: предполагаем, что утверждение верно для <tex>r = j</tex>. Раскрываем первое слагаемое в правой части предположения, аналогично базе и получаем:
<center><tex>{x \choose y} = {x - j - 1 \choose y} + {x - (j + 1) \choose y - 1} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1} = {x - (j + 1) \choose y} + \sum \limits_{i = 1}^{j}{x - i \choose y - 1}</tex>.</center>
}}

{{Лемма
|statement=<center><tex>{x \choose y} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
|proof=Подставляем в предыдущую лемму <tex>r = x - y</tex> и получаем:
<center><tex>{x \choose y} = {y \choose y} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = {x - (x - y + 1) \choose y - 1} + \sum \limits_{i = 1}^{x - y} {x - i \choose y - 1} = \sum \limits_{i = 1}^{x - y + 1} {x - i \choose y - 1}</tex>.</center>
}}

продолжаем доказательство теоремы по индукции:
<tex>f(n, d) = {n + d - 2 \choose d - 1} = \sum \limits_{i = 1}^{n} {n + d - 2 + i \choose d - 2} = \sum \limits_{k = n}^{1} {n + d + 2 - (n - k + 1) \choose d - 2} = \sum \limits_{k = 1}^{n} {k + d - 3 \choose d - 2} = \sum \limits_{k = 1}^{n}{k + (d - 1) - 2 \choose (d - 1) - 1} = \sum \limits_{k = 1}^{n}{f(k, d - 1)}</tex>

}}

Таким образом, на этот алгоритм получена более низкая оценка времени работы, чем на предыдущий, и это подтверждается на практике [1].

== Сведение к задаче KMP (Klee's Measure Problem) ==

Задача KMP состоит в нахождении объема объединения прямоугольных гиперпараллелепипедов в <tex>d</tex>-мерном пространстве. Как показано в описании алгоритма IEA, исходная задача легко сводится к этой, если каждой точке поставить в соответствие гиперкуб с одной вершиной в центре координат, а противоположной - в этой точке.

Существует несколько алгоритмов решения задачи KMP, самый оптимальный из которых использует идеи сканирующей гиперплоскости, [[Статистики на отрезках. Корневая эвристика|корневой эвристики]] и КД-дерево, и позволяет решать задачу за время <tex>O(n^{d/2}\log n)</tex>. Рассмотрим его.

Для начала изложим идею сканирующей гиперплоскости. Рассмотрим всевозможные значения <tex>d</tex>-ой координаты у точек множества <tex>X</tex>.
Возьмем гиперплоскость, перпендикулярную оси <tex>d</tex>-ой координаты, будем ее двигать по этим значениям от минимального до максимального и рассматривать проекцию всех гиперкубов на эту гиперплоскость. Заметим, что между точками проекция меняться не будет, а в каждой точке будут появляться или исчезать проекции некоторых гиперкубов. Поэтому, если уметь эффективно пересчитывать объем проекции при каждом появлении/исчезновении и прибавлять это значение, умноженное на расстояние между последовательными значениями <tex>d</tex>-ой координаты, то можно легко посчитать суммарный гиперобъем. Как раз для эффективного пересчета используется такая структура, как КД-дерево.

Теперь изложим идею алгоритма на примере трехмерного пространства.
Будем считать, что сканирующая плоскость двигается перепендикулярно оси <tex>OZ</tex>. В каждый момент проекция всех параллелепипедов на эту плоскость предствляет собой множество прямоугольников и необходимо уметь обрабатывать следующие три запроса:
# добавить прямоугольник
# удалить прямоугольник
# посчитать площадь объединения всех прямоугольников

{{Определение
|definition=КД-деревом называется [[Дерево поиска, наивная реализация|сбалансированное двоичное дерево]], где каждой вершине <tex>\alpha</tex> соответствует некоторая область <tex>C_{\alpha}</tex> <tex>d</tex>-мерного пространстве, удовлетворяющая следующим свойствам:
# Корневой вершине <tex>root</tex> соответствует <tex>C_{root}</tex> - все пространство.
# Каждое <tex>C_{\alpha}</tex> является гиперпараллелепипедом (возможно бесконечным в некоторых направлениях).
# Для каждой вершины <tex>\alpha</tex> с детьми <tex>\beta</tex> и <tex>\gamma</tex> верно, что области <tex>C_{\beta}</tex> и <tex>C_{\gamma}</tex> не пересекаются и <tex>C_{\alpha} = C_{\beta} \cup C_{\gamma}</tex>.
}}

Будем хранить прямоугольники в КД-дереве следующим образом:
# Для каждого листа <tex>\alpha</tex> будем хранить все прямоугольники, которые пересекают внутреннюю часть области <tex>C_{\alpha}</tex>.
# Для каждой внутренней вершины <tex>\alpha</tex> будем считать число <tex>Cover_{\alpha}</tex> - количество прямоугольников, которые полностью покрывают область <tex>C_{\alpha}</tex>, но не полностью покрывают область родителя этой вершины <tex>C_{father(\alpha)}</tex>.

Также для каждой вершины будем хранить значение площади пересечения области этой вершины со всеми прямоугольниками <tex>M</tex>, которое определяется следующим образом:
# Для каждого листа <tex>M</tex> равно площади пересечения всех прямоугольников с областью этого листа.
# Для каждой внутренней вершины <tex>\alpha</tex> если <tex>Cover_{\alpha} > 0</tex>, тогда <tex>M</tex> равно площади всех области <tex>C_{\alpha}</tex>, иначе <tex>M</tex> равно сумме этих значений для двух сыновей этой вершины.

Таким образом, значение для корневой вершины <tex>M_{root}</tex> и является искомой площадью объединения всех прямоугольников.

Рассмотрим операцию добавления прямоугольника (box) в КД-дерево.

procedure Insert(box, <tex>\alpha</tex>)
if <tex>\alpha</tex> - лист then
добавить box в эту вершину
пересчитать <tex>M_{\alpha}</tex>
elseif box покрывает область <tex>C_{\alpha}</tex> then
<tex>Cover_{\alpha} = Cover_{\alpha} + 1</tex>
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
elseif box пересекает область <tex>C_{\alpha}</tex> then
Insert(box, leftson(<tex>\alpha</tex>))
Insert(box, rightson(<tex>\alpha</tex>))
if <tex>Cover_{\alpha} > 0 </tex> then
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
else
<tex>M_{\alpha} = M_{leftson(\alpha)} + M_{rightson(\alpha)}</tex>

Аналогично устроена операция удаления прямоугольника - она выполняет обратные действия.

Рассмотрим теперь как именно будет устроено разделение пространства на области КД-деревом для эффективного его использования.
Возьмем проекции всех параллелепипедов на плоскость, перепендикулярную оси <tex>OZ</tex> - множество прямоугольников <tex>V</tex>. Рассмотрим все <tex>x</tex>-координаты границ прямоугольников и разобьем ось <tex>OX</tex> на <tex>2\sqrt{n}</tex> промежутков так, чтобы в каждом из них было не более <tex>\sqrt{n}</tex> границ. Получившиемя <tex>2\sqrt{n}</tex> полос разобьем горизонтальными линиями на итоговые области по отдельности следующим образом. Для каждой полосы у всех прямоугольников, имеющих вертикальную границу в этой полосе возьмем все горизонтальные границы (их будет не более <tex>2\sqrt{n}</tex>), а среди всех горизонтальных границ, пересекающих область, возьмем каждую <tex>\sqrt{n}</tex>-ую. Таким образом, в итоге получим не более <tex>4\sqrt{n}</tex> областей в каждой верикальной полосе.

Полученные области обладают следующими свойствами:
# Всего <tex>O(n)</tex> областей.
# Каждый прямоугольник частично покрывает не более <tex>O(\sqrt{n})</tex> областей.
# Никакая вершина прямоугольника не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти утверждения несложно следуют из алгоритма разбиения на области.

Каждая из этих областей будет одним из листов КД-дерева. Для построения всего дерева сначала объединим все области для одной вертикальной полосы, сливая соседние области в одну вершину. После этого аналогично будем сливать соседние вертикальные области и в итоге получим корневую вершину, представляющую всю плоскость.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n)</tex> вершин.
# Каждый прямоугольник хранится в <tex>O(\sqrt{n})</tex> листьях.
# Каждый прямоугольник влияет на <tex>O(\sqrt{n} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина прямоугольника не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти свойства следуют из свойств разбиения на области и логарифмической высоты дерева.

Наконец, получаем итоговый алгоритм. Идем сканирующей плоскостью, перепендикулярной оси <tex>OZ</tex>, добавляя или удаляя прямоугольники из КД-дерева. При этом каждый раз пересчитывается площадь объединения всех текущих прямоугольников и, складывая эти величины, умноженные на соответствующую разницу <tex>z</tex>-координат получаем общий объем.

Каждая операция добавления или удаления прямоугольников требует изменения <tex>O(\sqrt{n} \log n)</tex> полей <tex>Cover</tex>, а также операций добавления или удаления прямоугольника в <tex>O(\sqrt{n})</tex> листьев, каждая из которых требует <tex>O(\log n)</tex> времени. Итого, на одну операцию добавления или удаления требуется <tex>O(\sqrt{n} \log n)</tex> времени, суммарно таких операций <tex>O(n)</tex>, и получается, что весь алгоритм работает за <tex>O(n \sqrt{n} \log n)</tex>. При этом на хранение всей информации требуется <tex>O(n\sqrt{n}) </tex> памяти - каждый лист может хранить <tex>O(\sqrt{n})</tex> прямоугольников.

Все идеи этого алгоритма обобщаются с трехмерного пространства на <tex>d</tex>-мерное следующим образом.

Построение разбиения на области в <tex>d</tex>-мерном пространстве вместо двухмерного ведется сначала аналогичным образом, только после того, как в двухмерном пространстве получаются итоговые прямоугольные области, здесь продолжаются аналогичные разбиения вдоль третьей координаты - берутся все границы третьей координаты гиперпараллелепипедов, имеющие границу по первой или второй координате строго внутри области (их получается не более <tex>O(\sqrt{n})</tex>), а также каждую <tex>O(\sqrt{n})</tex>-ую границу по третьей координате гиперепараллелепипедов, пересекающих эту область, и так далее.

Это разбиение на области имеет следующие аналогичные свойства:
# Всего <tex>O(n^{d/2})</tex> областей.
# Каждый гиперпараллелепипед частично покрывает не более <tex>O(n^{(d-1)/2})</tex> областей.
# Никакая вершина гиперпараллелепипеда не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Общее дерево строится аналогичным образом и оно в итоге состоит из <tex>d</tex> уровней, каждый из которых имеет высоту <tex>O(\log n)</tex> в дереве - самый верхний уровень содержит разбиение по первой координате, следующий - по второй и так далее.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n^{d/2})</tex> вершин.
# Каждый гиперпараллелепипед хранится в <tex>O(n^{(d-1)/2})</tex> листьях.
# Каждый гиперпараллелепипед влияет на <tex>O(n^{(d-1)/2} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина гиперпараллелепипеда не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Операции в этом дереве выполняются аналогично трехмерному случаю и каждая операция добавления или удаления в КД-дерево размерности <tex>d</tex> занимает <tex>O(n^{(d-1)/2} \log n)</tex> времени. Выполняя сканирование гиперплоскостью вдоль одной из координат получаем КД-дерево размерности <tex>d-1</tex>, в котором необходимо провести <tex>O(n)</tex> операций, то есть суммарное время работы алгоритма равно: <tex>O(n \times n^{(d-1-1)/2} \log n) = O(n^{d/2} \log n)</tex>. При этом хранение КД-дерева требует <tex>O(n^{d/2})</tex> памяти.

== Источники ==
# While, L., Hingston, P., Barone, L., Huband, S.: A Faster Algorithm for Calculating Hypervolume. IEEE Transaction on Evolutionary Computation, Vol.10, No. 1, pp 29–38 (2006)
# [http://rain.ifmo.ru/~tsarev/teaching/ea-2012/seminars/1320.pdf Zhou X., Sun C., Mao N., Li W.: Generalization of HSO algortihms for computing hypervolume for mutiobjective optimization problems, IEEE Congress on Evolutionary Computation (2007)]
# K. Bringmann, T. Friedrich.: An Efficient Algorithm for Computing Hypervolume Contributions. Evolutionary Computation, Vol. 18, No. 3, pp 383-402, MIT Press. (2010)
# N. Beume.: S-Metric calculation by considering dominated hypervolume as Klee’s measure problem. Evolutionary Computation, 17(4), pp 477–492. (2009)
# N. Beume, G. Rudolph.: Faster S-metric calculation by considering dominated hypervolume as Klee’s measure problem. In Proc. Second International Conference on Computational Intelligence (IASTED ’06), pp 233–238. (2006)
# Overmars, M.H., Yap, C.K.: New upper bounds in Klee’s measure problem. SIAM Journal on Computing 20(6), pp 1034–1045. (1991)
# Klee, V.: Can the measure of S[ai, bi] be computed in less than O(n log n) steps? In: American Mathematical Monthly. Volume 84, pp 284–285. (1977)

Алгоритмы точного вычисления гиперобъема

2012-06-17T22:35:37Z

Joshik:

{{В разработке}}

== Постановка задачи ==
<tex>x = (x_1, x_2, ..., x_d; x_i \ge 0) \in R^d</tex> - точка в <tex>d</tex>-мерном пространстве.

Точка <tex>x</tex> доминирует точку <tex>y</tex> (<tex>x \succ y</tex>), если <tex>\forall i : x_i \ge y_i, \exists j : x_j > y_j</tex>.

<tex>X = (x^1, x^2, ..., x^n) \subset R^d</tex> - множество из <tex>n</tex> точек в <tex>d</tex>-мерном пространстве таких, что <tex>\nexists i \neq j : x_i \succ x_j</tex> - никакая точка не доминируется другой точкой из этого множества.

<tex>S(X) = \mu (\bigcup \limits_{x \in X} \{y | y \succ x\}) </tex> - гиперобъем множества <tex>X</tex>.

В частности, если <tex>X = \{x\}</tex>, то <tex>S(X) = \prod \limits_{i=1}^{d} x_i</tex>.

Задача: найти точное значение гиперобъема <tex>S(X)</tex> множества из <tex>n</tex> точек <tex>d</tex>-мерного пространоства.

== Алгоритм включения-исключения (Inclusion-Exclusion Algorithm, IEA) ==

Самый простой алгоритм нахождения гиперобъема базируется на идее комбинаторной [[формула включения-исключения|формулы включения-искючения]].
Все множество <tex>X</tex> представляется в виде объединения <tex>n</tex> гиперкубов (<tex>X^i</tex>), соответствующих отдельным точкам <tex>x^i</tex>.

После этого объем всего множества вычисляется по формуле: <center> <tex> S(X) = \sum \limits_{I \in 2^n} (-1)^{|I|+1} S(\bigcap \limits_{ j \in I} X^j) </tex> </center>

Объем пересечения гиперкубов легко считается как произведение по каждой координате минимального значения этой координаты среди всех точек, которым соответствуют гиперкубы.

Таким образом, в этом алгоритме перебираются все подмножества точек множества <tex>X</tex>, для каждого множества находится гиперобъем пересечения соответствующих гиперкубов и он прибавляется с соответствующим знаком к ответу.
Время работы этого алгоритма составляет <tex>O(n 2^n)</tex>.

== Алгоритм LebMeasure ==

Алгоритм LebMeasure обрабатывает точки множества <tex>X</tex> по очереди. Для каждой очередной точки <tex>x</tex> находится объем некоторого максимального по включению гиперкуба, доминируемого эксклюзивно только этой точкой <tex>x</tex> и она заменяется на некоторое множество ''порожденных'' точек, которые доминируют оставшуюся область, доминировшуюся до этого точкой <tex>x</tex>.

Например, если изначально было четыре точки в трехмерном пространстве <center><tex>(6, 7, 4), (9, 5, 5), (1, 9, 3), (4, 1, 9)</tex>, </center> то точка <tex>(6, 7, 4)</tex> эксклюзивно доминирует куб с одним концов в <tex>(6, 7, 4)</tex>, а другим - в <tex>(4, 5, 3)</tex>. После добавления объема этого куба в ответу, точка <tex>(6, 7, 4)</tex> порождает три точки: <center><tex>(4, 7, 4), (6, 5, 3), (6, 7, 3)</tex>.</center>При этом, точка <tex>(6, 5, 4)</tex> доминируется точкой <tex>(9, 5, 5)</tex> и сразу удаляется из множества <tex>X</tex>.

Обработка точек продолжается, пока все точки не будут обработаны. Таким образом, время работы алгоритмы напрямую зависит от количества порожденных точек. Легко заметить, что таких точек не более, чем <tex>n^d</tex>, поскольку каждая координата каждой порожденной точки равна соответствующей координате некоторой точки исходного множества <tex>X</tex>.

К сожалению, эта верхняя оценка является достижимой. Например, если исходное множество <tex>X</tex> имеет вид: <tex>(1, n, n, ..., n), (2, n - 1, n - 1, ..., n - 1), ..., (n, 1, 1, ..., 1)</tex>, и точки обрабатываются в этом порядке, то всего будет обработано <tex>n^{d-1}</tex> точка, что показано в [1]. Правда, если в этом примере точки обрабатываются в обратном порядке, то суммарное количество обработанных точек линейно зависит от <tex>n</tex> и <tex>d</tex>. Тем не менее, существуют примеры, для которых любой порядок обработки приводит к экспоненциальной зависимости числа порожденных точек от размерности пространства <tex>d</tex> и близок к <tex>n^d</tex> [1].

== Алгоритм Hypervolume by Slicing Objectives (HSO) ==

Под ''Objectives'' в названии данного алгоритма имеются в виду координаты пространства <tex>R^d</tex>.

Если алгоритм LebMeasure по очереди рассматривает все точки, то алгоритм HSO рассматривает по очереди все координаты, сводя задачу к меньшей на единицу размерности.

Изначально все точки сортируются по первой координате. Значения этой координаты используются для ''расслоения''(разрезания) всего множества на <tex>n</tex> частей, внутри каждой из которых при движении вдоль первой координаты форма ''разреза'' перпендикулярно оси первой координаты не меняется. Таким образом, для подсчета объема каждой части необходимо найти объем разреза и умножить длину части вдоль первой координаты. При этом получившийся разрез имеет на единицу меньшую размерность.
Заметим, что после сортировки и расслоения первая часть содержит все <tex>n</tex> точек, вторая - все, помимо точки с минимальной координатой, вдоль которой происходит расслоение и т.д., а последняя часть содержит только одну точку с максимальной этой координатой.

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

Описанный процесс можно реализовать как в виде рекурсивной процедуры, расслаивающей множество вдоль очередной координаты и вызывающей себя рекурсивно для каждой части и следующей координаты, так и нерекурсивно, если поддерживать множество всех текущих частей и на очередной итерации разбивать их все вдоль очередной координаты.

Время работы алгоритма напрямую зависит от суммарного количества частей, на которые будет разбито исходное множество. По аналогичным предыдущему алгоритму рассуждениям легко показывается, что частей не более <tex>n^d</tex>. Тем не менее, существует более точная оценка.

{{Теорема
|statement=Суммарное количество частей, полученных алгоритмом HSO, не превышает <center><tex>f(n, d) = {n + d - 2 \choose d - 1}</tex>.</center>

|proof=
База: <center><tex>f(n, 1) = 1</tex>.</center>

Также, из алгоритма расслоения следует, что из части с <tex>n</tex> точками будет получено по одной части для каждого <tex>i = 1...n</tex>, то есть: <center><tex>f(n, d) = \sum \limits_{i = 1}^{n}f(i, d - 1)</tex>.</center>

Докажем полезные леммы:

{{Лемма
|statement=<center><tex>{x \choose y} = {x - r \choose y} + \sum \limits_{i = 1}{r} {x - i \choose y - 1}</tex>.</center>
|proof=Индукция по <tex>r</tex>. База: <center><tex>{x \choose y} = {x - 1 \choose y} + {x - 1 \choose y - 1}</tex>.</center>

}}

2)

}}

Таким образом, на этот алгоритм получена более низкая оценка времени работы, чем на предыдущий, и это подтверждается на практике [1].

== Сведение к задаче KMP (Klee's Measure Problem) ==

Задача KMP состоит в нахождении объема объединения прямоугольных гиперпараллелепипедов в <tex>d</tex>-мерном пространстве. Как показано в описании алгоритма IEA, исходная задача легко сводится к этой, если каждой точке поставить в соответствие гиперкуб с одной вершиной в центре координат, а противоположной - в этой точке.

Существует несколько алгоритмов решения задачи KMP, самый оптимальный из которых использует идеи сканирующей гиперплоскости, [[Статистики на отрезках. Корневая эвристика|корневой эвристики]] и КД-дерево, и позволяет решать задачу за время <tex>O(n^{d/2}\log n)</tex>. Рассмотрим его.

Для начала изложим идею сканирующей гиперплоскости. Рассмотрим всевозможные значения <tex>d</tex>-ой координаты у точек множества <tex>X</tex>.
Возьмем гиперплоскость, перпендикулярную оси <tex>d</tex>-ой координаты, будем ее двигать по этим значениям от минимального до максимального и рассматривать проекцию всех гиперкубов на эту гиперплоскость. Заметим, что между точками проекция меняться не будет, а в каждой точке будут появляться или исчезать проекции некоторых гиперкубов. Поэтому, если уметь эффективно пересчитывать объем проекции при каждом появлении/исчезновении и прибавлять это значение, умноженное на расстояние между последовательными значениями <tex>d</tex>-ой координаты, то можно легко посчитать суммарный гиперобъем. Как раз для эффективного пересчета используется такая структура, как КД-дерево.

Теперь изложим идею алгоритма на примере трехмерного пространства.
Будем считать, что сканирующая плоскость двигается перепендикулярно оси <tex>OZ</tex>. В каждый момент проекция всех параллелепипедов на эту плоскость предствляет собой множество прямоугольников и необходимо уметь обрабатывать следующие три запроса:
# добавить прямоугольник
# удалить прямоугольник
# посчитать площадь объединения всех прямоугольников

{{Определение
|definition=КД-деревом называется [[Дерево поиска, наивная реализация|сбалансированное двоичное дерево]], где каждой вершине <tex>\alpha</tex> соответствует некоторая область <tex>C_{\alpha}</tex> <tex>d</tex>-мерного пространстве, удовлетворяющая следующим свойствам:
# Корневой вершине <tex>root</tex> соответствует <tex>C_{root}</tex> - все пространство.
# Каждое <tex>C_{\alpha}</tex> является гиперпараллелепипедом (возможно бесконечным в некоторых направлениях).
# Для каждой вершины <tex>\alpha</tex> с детьми <tex>\beta</tex> и <tex>\gamma</tex> верно, что области <tex>C_{\beta}</tex> и <tex>C_{\gamma}</tex> не пересекаются и <tex>C_{\alpha} = C_{\beta} \cup C_{\gamma}</tex>.
}}

Будем хранить прямоугольники в КД-дереве следующим образом:
# Для каждого листа <tex>\alpha</tex> будем хранить все прямоугольники, которые пересекают внутреннюю часть области <tex>C_{\alpha}</tex>.
# Для каждой внутренней вершины <tex>\alpha</tex> будем считать число <tex>Cover_{\alpha}</tex> - количество прямоугольников, которые полностью покрывают область <tex>C_{\alpha}</tex>, но не полностью покрывают область родителя этой вершины <tex>C_{father(\alpha)}</tex>.

Также для каждой вершины будем хранить значение площади пересечения области этой вершины со всеми прямоугольниками <tex>M</tex>, которое определяется следующим образом:
# Для каждого листа <tex>M</tex> равно площади пересечения всех прямоугольников с областью этого листа.
# Для каждой внутренней вершины <tex>\alpha</tex> если <tex>Cover_{\alpha} > 0</tex>, тогда <tex>M</tex> равно площади всех области <tex>C_{\alpha}</tex>, иначе <tex>M</tex> равно сумме этих значений для двух сыновей этой вершины.

Таким образом, значение для корневой вершины <tex>M_{root}</tex> и является искомой площадью объединения всех прямоугольников.

Рассмотрим операцию добавления прямоугольника (box) в КД-дерево.

procedure Insert(box, <tex>\alpha</tex>)
if <tex>\alpha</tex> - лист then
добавить box в эту вершину
пересчитать <tex>M_{\alpha}</tex>
elseif box покрывает область <tex>C_{\alpha}</tex> then
<tex>Cover_{\alpha} = Cover_{\alpha} + 1</tex>
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
elseif box пересекает область <tex>C_{\alpha}</tex> then
Insert(box, leftson(<tex>\alpha</tex>))
Insert(box, rightson(<tex>\alpha</tex>))
if <tex>Cover_{\alpha} > 0 </tex> then
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
else
<tex>M_{\alpha} = M_{leftson(\alpha)} + M_{rightson(\alpha)}</tex>

Аналогично устроена операция удаления прямоугольника - она выполняет обратные действия.

Рассмотрим теперь как именно будет устроено разделение пространства на области КД-деревом для эффективного его использования.
Возьмем проекции всех параллелепипедов на плоскость, перепендикулярную оси <tex>OZ</tex> - множество прямоугольников <tex>V</tex>. Рассмотрим все <tex>x</tex>-координаты границ прямоугольников и разобьем ось <tex>OX</tex> на <tex>2\sqrt{n}</tex> промежутков так, чтобы в каждом из них было не более <tex>\sqrt{n}</tex> границ. Получившиемя <tex>2\sqrt{n}</tex> полос разобьем горизонтальными линиями на итоговые области по отдельности следующим образом. Для каждой полосы у всех прямоугольников, имеющих вертикальную границу в этой полосе возьмем все горизонтальные границы (их будет не более <tex>2\sqrt{n}</tex>), а среди всех горизонтальных границ, пересекающих область, возьмем каждую <tex>\sqrt{n}</tex>-ую. Таким образом, в итоге получим не более <tex>4\sqrt{n}</tex> областей в каждой верикальной полосе.

Полученные области обладают следующими свойствами:
# Всего <tex>O(n)</tex> областей.
# Каждый прямоугольник частично покрывает не более <tex>O(\sqrt{n})</tex> областей.
# Никакая вершина прямоугольника не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти утверждения несложно следуют из алгоритма разбиения на области.

Каждая из этих областей будет одним из листов КД-дерева. Для построения всего дерева сначала объединим все области для одной вертикальной полосы, сливая соседние области в одну вершину. После этого аналогично будем сливать соседние вертикальные области и в итоге получим корневую вершину, представляющую всю плоскость.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n)</tex> вершин.
# Каждый прямоугольник хранится в <tex>O(\sqrt{n})</tex> листьях.
# Каждый прямоугольник влияет на <tex>O(\sqrt{n} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина прямоугольника не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти свойства следуют из свойств разбиения на области и логарифмической высоты дерева.

Наконец, получаем итоговый алгоритм. Идем сканирующей плоскостью, перепендикулярной оси <tex>OZ</tex>, добавляя или удаляя прямоугольники из КД-дерева. При этом каждый раз пересчитывается площадь объединения всех текущих прямоугольников и, складывая эти величины, умноженные на соответствующую разницу <tex>z</tex>-координат получаем общий объем.

Каждая операция добавления или удаления прямоугольников требует изменения <tex>O(\sqrt{n} \log n)</tex> полей <tex>Cover</tex>, а также операций добавления или удаления прямоугольника в <tex>O(\sqrt{n})</tex> листьев, каждая из которых требует <tex>O(\log n)</tex> времени. Итого, на одну операцию добавления или удаления требуется <tex>O(\sqrt{n} \log n)</tex> времени, суммарно таких операций <tex>O(n)</tex>, и получается, что весь алгоритм работает за <tex>O(n \sqrt{n} \log n)</tex>. При этом на хранение всей информации требуется <tex>O(n\sqrt{n}) </tex> памяти - каждый лист может хранить <tex>O(\sqrt{n})</tex> прямоугольников.

Все идеи этого алгоритма обобщаются с трехмерного пространства на <tex>d</tex>-мерное следующим образом.

Построение разбиения на области в <tex>d</tex>-мерном пространстве вместо двухмерного ведется сначала аналогичным образом, только после того, как в двухмерном пространстве получаются итоговые прямоугольные области, здесь продолжаются аналогичные разбиения вдоль третьей координаты - берутся все границы третьей координаты гиперпараллелепипедов, имеющие границу по первой или второй координате строго внутри области (их получается не более <tex>O(\sqrt{n})</tex>), а также каждую <tex>O(\sqrt{n})</tex>-ую границу по третьей координате гиперепараллелепипедов, пересекающих эту область, и так далее.

Это разбиение на области имеет следующие аналогичные свойства:
# Всего <tex>O(n^{d/2})</tex> областей.
# Каждый гиперпараллелепипед частично покрывает не более <tex>O(n^{(d-1)/2})</tex> областей.
# Никакая вершина гиперпараллелепипеда не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Общее дерево строится аналогичным образом и оно в итоге состоит из <tex>d</tex> уровней, каждый из которых имеет высоту <tex>O(\log n)</tex> в дереве - самый верхний уровень содержит разбиение по первой координате, следующий - по второй и так далее.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n^{d/2})</tex> вершин.
# Каждый гиперпараллелепипед хранится в <tex>O(n^{(d-1)/2})</tex> листьях.
# Каждый гиперпараллелепипед влияет на <tex>O(n^{(d-1)/2} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина гиперпараллелепипеда не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Операции в этом дереве выполняются аналогично трехмерному случаю и каждая операция добавления или удаления в КД-дерево размерности <tex>d</tex> занимает <tex>O(n^{(d-1)/2} \log n)</tex> времени. Выполняя сканирование гиперплоскостью вдоль одной из координат получаем КД-дерево размерности <tex>d-1</tex>, в котором необходимо провести <tex>O(n)</tex> операций, то есть суммарное время работы алгоритма равно: <tex>O(n \times n^{(d-1-1)/2} \log n) = O(n^{d/2} \log n)</tex>. При этом хранение КД-дерева требует <tex>O(n^{d/2})</tex> памяти.

== Источники ==
# While, L., Hingston, P., Barone, L., Huband, S.: A Faster Algorithm for Calculating Hypervolume. IEEE Transaction on Evolutionary Computation, Vol.10, No. 1, pp 29–38 (2006)
# Zhou X., Sun C., Mao N., Li W.: Generalization of HSO algortihms for computing hypervolume for mutiobjective optimization problems, IEEE Congress on Evolutionary Computation (2007)
# K. Bringmann, T. Friedrich.: An Efficient Algorithm for Computing Hypervolume Contributions. Evolutionary Computation, Vol. 18, No. 3, pp 383-402, MIT Press. (2010)
# N. Beume.: S-Metric calculation by considering dominated hypervolume as Klee’s measure problem. Evolutionary Computation, 17(4), pp 477–492. (2009)
# N. Beume, G. Rudolph.: Faster S-metric calculation by considering dominated hypervolume as Klee’s measure problem. In Proc. Second International Conference on Computational Intelligence (IASTED ’06), pp 233–238. (2006)
# Overmars, M.H., Yap, C.K.: New upper bounds in Klee’s measure problem. SIAM Journal on Computing 20(6), pp 1034–1045. (1991)
# Klee, V.: Can the measure of S[ai, bi] be computed in less than O(n log n) steps? In: American Mathematical Monthly. Volume 84, pp 284–285. (1977)

Алгоритмы точного вычисления гиперобъема

2012-06-17T22:19:38Z

Joshik:

{{В разработке}}

== Постановка задачи ==
<tex>x = (x_1, x_2, ..., x_d; x_i \ge 0) \in R^d</tex> - точка в <tex>d</tex>-мерном пространстве.

Точка <tex>x</tex> доминирует точку <tex>y</tex> (<tex>x \succ y</tex>), если <tex>\forall i : x_i \ge y_i, \exists j : x_j > y_j</tex>.

<tex>X = (x^1, x^2, ..., x^n) \subset R^d</tex> - множество из <tex>n</tex> точек в <tex>d</tex>-мерном пространстве таких, что <tex>\nexists i \neq j : x_i \succ x_j</tex> - никакая точка не доминируется другой точкой из этого множества.

<tex>S(X) = \mu (\bigcup \limits_{x \in X} \{y | y \succ x\}) </tex> - гиперобъем множества <tex>X</tex>.

В частности, если <tex>X = \{x\}</tex>, то <tex>S(X) = \prod \limits_{i=1}^{d} x_i</tex>.

Задача: найти точное значение гиперобъема <tex>S(X)</tex> множества из <tex>n</tex> точек <tex>d</tex>-мерного пространоства.

== Алгоритм включения-исключения (Inclusion-Exclusion Algorithm, IEA) ==

Самый простой алгоритм нахождения гиперобъема базируется на идее комбинаторной [[формула включения-исключения|формулы включения-искючения]].
Все множество <tex>X</tex> представляется в виде объединения <tex>n</tex> гиперкубов (<tex>X^i</tex>), соответствующих отдельным точкам <tex>x^i</tex>.

После этого объем всего множества вычисляется по формуле: <center> <tex> S(X) = \sum \limits_{I \in 2^n} (-1)^{|I|+1} S(\bigcap \limits_{ j \in I} X^j) </tex> </center>

Объем пересечения гиперкубов легко считается как произведение по каждой координате минимального значения этой координаты среди всех точек, которым соответствуют гиперкубы.

Таким образом, в этом алгоритме перебираются все подмножества точек множества <tex>X</tex>, для каждого множества находится гиперобъем пересечения соответствующих гиперкубов и он прибавляется с соответствующим знаком к ответу.
Время работы этого алгоритма составляет <tex>O(n 2^n)</tex>.

== Алгоритм LebMeasure ==

Алгоритм LebMeasure обрабатывает точки множества <tex>X</tex> по очереди. Для каждой очередной точки <tex>x</tex> находится объем некоторого максимального по включению гиперкуба, доминируемого эксклюзивно только этой точкой <tex>x</tex> и она заменяется на некоторое множество ''порожденных'' точек, которые доминируют оставшуюся область, доминировшуюся до этого точкой <tex>x</tex>.

Например, если изначально было четыре точки в трехмерном пространстве <center><tex>(6, 7, 4), (9, 5, 5), (1, 9, 3), (4, 1, 9)</tex>, </center> то точка <tex>(6, 7, 4)</tex> эксклюзивно доминирует куб с одним концов в <tex>(6, 7, 4)</tex>, а другим - в <tex>(4, 5, 3)</tex>. После добавления объема этого куба в ответу, точка <tex>(6, 7, 4)</tex> порождает три точки: <center><tex>(4, 7, 4), (6, 5, 3), (6, 7, 3)</tex>.</center>При этом, точка <tex>(6, 5, 4)</tex> доминируется точкой <tex>(9, 5, 5)</tex> и сразу удаляется из множества <tex>X</tex>.

Обработка точек продолжается, пока все точки не будут обработаны. Таким образом, время работы алгоритмы напрямую зависит от количества порожденных точек. Легко заметить, что таких точек не более, чем <tex>n^d</tex>, поскольку каждая координата каждой порожденной точки равна соответствующей координате некоторой точки исходного множества <tex>X</tex>.

К сожалению, эта верхняя оценка является достижимой. Например, если исходное множество <tex>X</tex> имеет вид: <tex>(1, n, n, ..., n), (2, n - 1, n - 1, ..., n - 1), ..., (n, 1, 1, ..., 1)</tex>, и точки обрабатываются в этом порядке, то всего будет обработано <tex>n^{d-1}</tex> точка, что показано в [1]. Правда, если в этом примере точки обрабатываются в обратном порядке, то суммарное количество обработанных точек линейно зависит от <tex>n</tex> и <tex>d</tex>. Тем не менее, существуют примеры, для которых любой порядок обработки приводит к экспоненциальной зависимости числа порожденных точек от размерности пространства <tex>d</tex> и близок к <tex>n^d</tex> [1].

== Алгоритм Hypervolume by Slicing Objectives (HSO) ==

Под ''Objectives'' в названии данного алгоритма имеются в виду координаты пространства <tex>R^d</tex>.

Если алгоритм LebMeasure по очереди рассматривает все точки, то алгоритм HSO рассматривает по очереди все координаты, сводя задачу к меньшей на единицу размерности.

Изначально все точки сортируются по первой координате. Значения этой координаты используются для ''расслоения''(разрезания) всего множества на <tex>n</tex> частей, внутри каждой из которых при движении вдоль первой координаты форма ''разреза'' перпендикулярно оси первой координаты не меняется. Таким образом, для подсчета объема каждой части необходимо найти объем разреза и умножить длину части вдоль первой координаты. При этом получившийся разрез имеет на единицу меньшую размерность.
Заметим, что после сортировки и расслоения первая часть содержит все <tex>n</tex> точек, вторая - все, помимо точки с минимальной координатой, вдоль которой происходит расслоение и т.д., а последняя часть содержит только одну точку с максимальной этой координатой.

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

Описанный процесс можно реализовать как в виде рекурсивной процедуры, расслаивающей множество вдоль очередной координаты и вызывающей себя рекурсивно для каждой части и следующей координаты, так и нерекурсивно, если поддерживать множество всех текущих частей и на очередной итерации разбивать их все вдоль очередной координаты.

Время работы алгоритма напрямую зависит от суммарного количества частей, на которые будет разбито исходное множество. По аналогичным предыдущему алгоритму рассуждениям легко показывается, что частей не более <tex>n^d</tex>. Тем не менее, существует более точная оценка.

{{Теорема
|statement=Суммарное количество частей, полученных алгоритмом HSO, не превышает <center><tex>f(n, d) = {n + d - 2 \choose d - 1}</tex>.</center>

|proof=
База: <center><tex>f(n, 1) = 1</tex>.</center>

Также, из алгоритма расслоения следует, что из части с <tex>n</tex> точками будет получено по одной части для каждого <tex>i = 1...n</tex>, то есть: <center><tex>f(n, d) = \sum \limits_{i = 1}^{n}f(i, d - 1)</tex>.</center>

1) Покажем, что: <center><tex>{x \choose y} = {x - r \choose y} + \sum \limits_{i = 1}{r} {x - i \choose y - 1}</tex>.</center>
Индукция по <tex>r</tex>. База: <center><tex>{x \choose y} = {x - 1 \choose y} + {x - 1 \choose y - 1}</tex>.</center>

2)

}}

Таким образом, на этот алгоритм получена более низкая оценка времени работы, чем на предыдущий, и это подтверждается на практике [1].

== Сведение к задаче KMP (Klee's Measure Problem) ==

Задача KMP состоит в нахождении объема объединения прямоугольных гиперпараллелепипедов в <tex>d</tex>-мерном пространстве. Как показано в описании алгоритма IEA, исходная задача легко сводится к этой, если каждой точке поставить в соответствие гиперкуб с одной вершиной в центре координат, а противоположной - в этой точке.

Существует несколько алгоритмов решения задачи KMP, самый оптимальный из которых использует идеи сканирующей гиперплоскости, [[Статистики на отрезках. Корневая эвристика|корневой эвристики]] и КД-дерево, и позволяет решать задачу за время <tex>O(n^{d/2}\log n)</tex>. Рассмотрим его.

Для начала изложим идею сканирующей гиперплоскости. Рассмотрим всевозможные значения <tex>d</tex>-ой координаты у точек множества <tex>X</tex>.
Возьмем гиперплоскость, перпендикулярную оси <tex>d</tex>-ой координаты, будем ее двигать по этим значениям от минимального до максимального и рассматривать проекцию всех гиперкубов на эту гиперплоскость. Заметим, что между точками проекция меняться не будет, а в каждой точке будут появляться или исчезать проекции некоторых гиперкубов. Поэтому, если уметь эффективно пересчитывать объем проекции при каждом появлении/исчезновении и прибавлять это значение, умноженное на расстояние между последовательными значениями <tex>d</tex>-ой координаты, то можно легко посчитать суммарный гиперобъем. Как раз для эффективного пересчета используется такая структура, как КД-дерево.

Теперь изложим идею алгоритма на примере трехмерного пространства.
Будем считать, что сканирующая плоскость двигается перепендикулярно оси <tex>OZ</tex>. В каждый момент проекция всех параллелепипедов на эту плоскость предствляет собой множество прямоугольников и необходимо уметь обрабатывать следующие три запроса:
# добавить прямоугольник
# удалить прямоугольник
# посчитать площадь объединения всех прямоугольников

{{Определение
|definition=КД-деревом называется [[Дерево поиска, наивная реализация|сбалансированное двоичное дерево]], где каждой вершине <tex>\alpha</tex> соответствует некоторая область <tex>C_{\alpha}</tex> <tex>d</tex>-мерного пространстве, удовлетворяющая следующим свойствам:
# Корневой вершине <tex>root</tex> соответствует <tex>C_{root}</tex> - все пространство.
# Каждое <tex>C_{\alpha}</tex> является гиперпараллелепипедом (возможно бесконечным в некоторых направлениях).
# Для каждой вершины <tex>\alpha</tex> с детьми <tex>\beta</tex> и <tex>\gamma</tex> верно, что области <tex>C_{\beta}</tex> и <tex>C_{\gamma}</tex> не пересекаются и <tex>C_{\alpha} = C_{\beta} \cup C_{\gamma}</tex>.
}}

Будем хранить прямоугольники в КД-дереве следующим образом:
# Для каждого листа <tex>\alpha</tex> будем хранить все прямоугольники, которые пересекают внутреннюю часть области <tex>C_{\alpha}</tex>.
# Для каждой внутренней вершины <tex>\alpha</tex> будем считать число <tex>Cover_{\alpha}</tex> - количество прямоугольников, которые полностью покрывают область <tex>C_{\alpha}</tex>, но не полностью покрывают область родителя этой вершины <tex>C_{father(\alpha)}</tex>.

Также для каждой вершины будем хранить значение площади пересечения области этой вершины со всеми прямоугольниками <tex>M</tex>, которое определяется следующим образом:
# Для каждого листа <tex>M</tex> равно площади пересечения всех прямоугольников с областью этого листа.
# Для каждой внутренней вершины <tex>\alpha</tex> если <tex>Cover_{\alpha} > 0</tex>, тогда <tex>M</tex> равно площади всех области <tex>C_{\alpha}</tex>, иначе <tex>M</tex> равно сумме этих значений для двух сыновей этой вершины.

Таким образом, значение для корневой вершины <tex>M_{root}</tex> и является искомой площадью объединения всех прямоугольников.

Рассмотрим операцию добавления прямоугольника (box) в КД-дерево.

procedure Insert(box, <tex>\alpha</tex>)
if <tex>\alpha</tex> - лист then
добавить box в эту вершину
пересчитать <tex>M_{\alpha}</tex>
elseif box покрывает область <tex>C_{\alpha}</tex> then
<tex>Cover_{\alpha} = Cover_{\alpha} + 1</tex>
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
elseif box пересекает область <tex>C_{\alpha}</tex> then
Insert(box, leftson(<tex>\alpha</tex>))
Insert(box, rightson(<tex>\alpha</tex>))
if <tex>Cover_{\alpha} > 0 </tex> then
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
else
<tex>M_{\alpha} = M_{leftson(\alpha)} + M_{rightson(\alpha)}</tex>

Аналогично устроена операция удаления прямоугольника - она выполняет обратные действия.

Рассмотрим теперь как именно будет устроено разделение пространства на области КД-деревом для эффективного его использования.
Возьмем проекции всех параллелепипедов на плоскость, перепендикулярную оси <tex>OZ</tex> - множество прямоугольников <tex>V</tex>. Рассмотрим все <tex>x</tex>-координаты границ прямоугольников и разобьем ось <tex>OX</tex> на <tex>2\sqrt{n}</tex> промежутков так, чтобы в каждом из них было не более <tex>\sqrt{n}</tex> границ. Получившиемя <tex>2\sqrt{n}</tex> полос разобьем горизонтальными линиями на итоговые области по отдельности следующим образом. Для каждой полосы у всех прямоугольников, имеющих вертикальную границу в этой полосе возьмем все горизонтальные границы (их будет не более <tex>2\sqrt{n}</tex>), а среди всех горизонтальных границ, пересекающих область, возьмем каждую <tex>\sqrt{n}</tex>-ую. Таким образом, в итоге получим не более <tex>4\sqrt{n}</tex> областей в каждой верикальной полосе.

Полученные области обладают следующими свойствами:
# Всего <tex>O(n)</tex> областей.
# Каждый прямоугольник частично покрывает не более <tex>O(\sqrt{n})</tex> областей.
# Никакая вершина прямоугольника не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти утверждения несложно следуют из алгоритма разбиения на области.

Каждая из этих областей будет одним из листов КД-дерева. Для построения всего дерева сначала объединим все области для одной вертикальной полосы, сливая соседние области в одну вершину. После этого аналогично будем сливать соседние вертикальные области и в итоге получим корневую вершину, представляющую всю плоскость.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n)</tex> вершин.
# Каждый прямоугольник хранится в <tex>O(\sqrt{n})</tex> листьях.
# Каждый прямоугольник влияет на <tex>O(\sqrt{n} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина прямоугольника не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти свойства следуют из свойств разбиения на области и логарифмической высоты дерева.

Наконец, получаем итоговый алгоритм. Идем сканирующей плоскостью, перепендикулярной оси <tex>OZ</tex>, добавляя или удаляя прямоугольники из КД-дерева. При этом каждый раз пересчитывается площадь объединения всех текущих прямоугольников и, складывая эти величины, умноженные на соответствующую разницу <tex>z</tex>-координат получаем общий объем.

Каждая операция добавления или удаления прямоугольников требует изменения <tex>O(\sqrt{n} \log n)</tex> полей <tex>Cover</tex>, а также операций добавления или удаления прямоугольника в <tex>O(\sqrt{n})</tex> листьев, каждая из которых требует <tex>O(\log n)</tex> времени. Итого, на одну операцию добавления или удаления требуется <tex>O(\sqrt{n} \log n)</tex> времени, суммарно таких операций <tex>O(n)</tex>, и получается, что весь алгоритм работает за <tex>O(n \sqrt{n} \log n)</tex>. При этом на хранение всей информации требуется <tex>O(n\sqrt{n}) </tex> памяти - каждый лист может хранить <tex>O(\sqrt{n})</tex> прямоугольников.

Все идеи этого алгоритма обобщаются с трехмерного пространства на <tex>d</tex>-мерное следующим образом.

Построение разбиения на области в <tex>d</tex>-мерном пространстве вместо двухмерного ведется сначала аналогичным образом, только после того, как в двухмерном пространстве получаются итоговые прямоугольные области, здесь продолжаются аналогичные разбиения вдоль третьей координаты - берутся все границы третьей координаты гиперпараллелепипедов, имеющие границу по первой или второй координате строго внутри области (их получается не более <tex>O(\sqrt{n})</tex>), а также каждую <tex>O(\sqrt{n})</tex>-ую границу по третьей координате гиперепараллелепипедов, пересекающих эту область, и так далее.

Это разбиение на области имеет следующие аналогичные свойства:
# Всего <tex>O(n^{d/2})</tex> областей.
# Каждый гиперпараллелепипед частично покрывает не более <tex>O(n^{(d-1)/2})</tex> областей.
# Никакая вершина гиперпараллелепипеда не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Общее дерево строится аналогичным образом и оно в итоге состоит из <tex>d</tex> уровней, каждый из которых имеет высоту <tex>O(\log n)</tex> в дереве - самый верхний уровень содержит разбиение по первой координате, следующий - по второй и так далее.

Построенное КД-дерево обладает следующими свойствами:
# В нем <tex>O(n^{d/2})</tex> вершин.
# Каждый гиперпараллелепипед хранится в <tex>O(n^{(d-1)/2})</tex> листьях.
# Каждый гиперпараллелепипед влияет на <tex>O(n^{(d-1)/2} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина гиперпараллелепипеда не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее гиперпараллелепипедов.

Операции в этом дереве выполняются аналогично трехмерному случаю и каждая операция добавления или удаления в КД-дерево размерности <tex>d</tex> занимает <tex>O(n^{(d-1)/2} \log n)</tex> времени. Выполняя сканирование гиперплоскостью вдоль одной из координат получаем КД-дерево размерности <tex>d-1</tex>, в котором необходимо провести <tex>O(n)</tex> операций, то есть суммарное время работы алгоритма равно: <tex>O(n \times n^{(d-1-1)/2} \log n) = O(n^{d/2} \log n)</tex>. При этом хранение КД-дерева требует <tex>O(n^{d/2})</tex> памяти.

Алгоритмы точного вычисления гиперобъема

2012-06-17T21:49:06Z

Joshik:

{{В разработке}}

== Постановка задачи ==
<tex>x = (x_1, x_2, ..., x_d; x_i \ge 0) \in R^d</tex> - точка в <tex>d</tex>-мерном пространстве.

Точка <tex>x</tex> доминирует точку <tex>y</tex> (<tex>x \succ y</tex>), если <tex>\forall i : x_i \ge y_i, \exists j : x_j > y_j</tex>.

<tex>X = (x^1, x^2, ..., x^n) \subset R^d</tex> - множество из <tex>n</tex> точек в <tex>d</tex>-мерном пространстве таких, что <tex>\nexists i \neq j : x_i \succ x_j</tex> - никакая точка не доминируется другой точкой из этого множества.

<tex>S(X) = \mu (\bigcup \limits_{x \in X} \{y | y \succ x\}) </tex> - гиперобъем множества <tex>X</tex>.

В частности, если <tex>X = \{x\}</tex>, то <tex>S(X) = \prod \limits_{i=1}^{d} x_i</tex>.

Задача: найти точное значение гиперобъема <tex>S(X)</tex> множества из <tex>n</tex> точек <tex>d</tex>-мерного пространоства.

== Алгоритм включения-исключения (Inclusion-Exclusion Algorithm, IEA) ==

Самый простой алгоритм нахождения гиперобъема базируется на идее комбинаторной [[формула включения-исключения|формулы включения-искючения]].
Все множество <tex>X</tex> представляется в виде объединения <tex>n</tex> гиперкубов (<tex>X^i</tex>), соответствующих отдельным точкам <tex>x^i</tex>.

После этого объем всего множества вычисляется по формуле: <center> <tex> S(X) = \sum \limits_{I \in 2^n} (-1)^{|I|+1} S(\bigcap \limits_{ j \in I} X^j) </tex> </center>

Объем пересечения гиперкубов легко считается как произведение по каждой координате минимального значения этой координаты среди всех точек, которым соответствуют гиперкубы.

Таким образом, в этом алгоритме перебираются все подмножества точек множества <tex>X</tex>, для каждого множества находится гиперобъем пересечения соответствующих гиперкубов и он прибавляется с соответствующим знаком к ответу.
Время работы этого алгоритма составляет <tex>O(n 2^n)</tex>.

== Алгоритм LebMeasure ==

Алгоритм LebMeasure обрабатывает точки множества <tex>X</tex> по очереди. Для каждой очередной точки <tex>x</tex> находится объем некоторого максимального по включению гиперкуба, доминируемого эксклюзивно только этой точкой <tex>x</tex> и она заменяется на некоторое множество ''порожденных'' точек, которые доминируют оставшуюся область, доминировшуюся до этого точкой <tex>x</tex>.

Например, если изначально было четыре точки в трехмерном пространстве <center><tex>(6, 7, 4), (9, 5, 5), (1, 9, 3), (4, 1, 9)</tex>, </center> то точка <tex>(6, 7, 4)</tex> эксклюзивно доминирует куб с одним концов в <tex>(6, 7, 4)</tex>, а другим - в <tex>(4, 5, 3)</tex>. После добавления объема этого куба в ответу, точка <tex>(6, 7, 4)</tex> порождает три точки: <center><tex>(4, 7, 4), (6, 5, 3), (6, 7, 3)</tex>.</center>При этом, точка <tex>(6, 5, 4)</tex> доминируется точкой <tex>(9, 5, 5)</tex> и сразу удаляется из множества <tex>X</tex>.

Обработка точек продолжается, пока все точки не будут обработаны. Таким образом, время работы алгоритмы напрямую зависит от количества порожденных точек. Легко заметить, что таких точек не более, чем <tex>n^d</tex>, поскольку каждая координата каждой порожденной точки равна соответствующей координате некоторой точки исходного множества <tex>X</tex>.

К сожалению, эта верхняя оценка является достижимой. Например, если исходное множество <tex>X</tex> имеет вид: <tex>(1, n, n, ..., n), (2, n - 1, n - 1, ..., n - 1), ..., (n, 1, 1, ..., 1)</tex>, и точки обрабатываются в этом порядке, то всего будет обработано <tex>n^{d-1}</tex> точка, что показано в [1]. Правда, если в этом примере точки обрабатываются в обратном порядке, то суммарное количество обработанных точек линейно зависит от <tex>n</tex> и <tex>d</tex>. Тем не менее, существуют примеры, для которых любой порядок обработки приводит к экспоненциальной зависимости числа порожденных точек от размерности пространства <tex>d</tex> и близок к <tex>n^d</tex> [1].

== Алгоритм Hypervolume by Slicing Objectives (HSO) ==

Под ''Objectives'' в названии данного алгоритма имеются в виду координаты пространства <tex>R^d</tex>.

Если алгоритм LebMeasure по очереди рассматривает все точки, то алгоритм HSO рассматривает по очереди все координаты, сводя задачу к меньшей на единицу размерности.

Изначально все точки сортируются по первой координате. Значения этой координаты используются для ''расслоения''(разрезания) всего множества на <tex>n</tex> частей, внутри каждой из которых при движении вдоль первой координаты форма ''разреза'' перпендикулярно оси первой координаты не меняется. Таким образом, для подсчета объема каждой части необходимо найти объем разреза и умножить длину части вдоль первой координаты. При этом получившийся разрез имеет на единицу меньшую размерность.
Заметим, что после сортировки и расслоения первая часть содержит все <tex>n</tex> точек, вторая - все, помимо точки с минимальной координатой, вдоль которой происходит расслоение и т.д., а последняя часть содержит только одну точку с максимальной этой координатой.

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

Описанный процесс можно реализовать как в виде рекурсивной процедуры, расслаивающей множество вдоль очередной координаты и вызывающей себя рекурсивно для каждой части и следующей координаты, так и нерекурсивно, если поддерживать множество всех текущих частей и на очередной итерации разбивать их все вдоль очередной координаты.

Время работы алгоритма напрямую зависит от суммарного количества частей, на которые будет разбито исходное множество. По аналогичным предыдущему алгоритму рассуждениям легко показывается, что частей не более <tex>n^d</tex>. Тем не менее, существует более точная оценка.

{{Теорема
|statement=Суммарное количество частей, полученных алгоритмом HSO, не превышает <center><tex>f(n, d) = {n + d - 2 \choose d - 1}</tex>.</center>

|proof=
База: <center><tex>f(n, 1) = 1</tex>.</center>

Также, из алгоритма расслоения следует, что из части с <tex>n</tex> точками будет получено по одной части для каждого <tex>i = 1...n</tex>, то есть: <center><tex>f(n, d) = \sum \limits_{i = 1}^{n}f(i, d - 1)</tex>.</center>

1) Покажем, что: <center><tex>{x \choose y} = {x - r \choose y} + \sum \limits_{i = 1}{r} {x - i \choose y - 1}</tex>.</center>
Индукция по <tex>r</tex>. База: <center><tex>{x \choose y} = {x - 1 \choose y} + {x - 1 \choose y - 1}</tex>.</center>

2)

}}

Таким образом, на этот алгоритм получена более низкая оценка времени работы, чем на предыдущий, и это подтверждается на практике [1].

== Сведение к задаче KMP (Klee's Measure Problem) ==

Задача KMP состоит в нахождении объема объединения прямоугольных гиперпараллелепипедов в <tex>d</tex>-мерном пространстве. Как показано в описании алгоритма IEA, исходная задача легко сводится к этой, если каждой точке поставить в соответствие гиперкуб с одной вершиной в центре координат, а противоположной - в этой точке.

Существует несколько алгоритмов решения задачи KMP, самый оптимальный из которых использует идеи сканирующей гиперплоскости, [[Статистики на отрезках. Корневая эвристика|корневой эвристики]] и КД-дерево, и позволяет решать задачу за время <tex>O(n^{d/2}\log n)</tex>. Рассмотрим его.

Для начала изложим идею сканирующей гиперплоскости. Рассмотрим всевозможные значения <tex>d</tex>-ой координаты у точек множества <tex>X</tex>.
Возьмем гиперплоскость, перпендикулярную оси <tex>d</tex>-ой координаты, будем ее двигать по этим значениям от минимального до максимального и рассматривать проекцию всех гиперкубов на эту гиперплоскость. Заметим, что между точками проекция меняться не будет, а в каждой точке будут появляться или исчезать проекции некоторых гиперкубов. Поэтому, если уметь эффективно пересчитывать объем проекции при каждом появлении/исчезновении и прибавлять это значение, умноженное на расстояние между последовательными значениями <tex>d</tex>-ой координаты, то можно легко посчитать суммарный гиперобъем. Как раз для эффективного пересчета используется такая структура, как КД-дерево.

Теперь изложим идею алгоритма на примере трехмерного пространства.
Будем считать, что сканирующая плоскость двигается перепендикулярно оси <tex>OZ</tex>. В каждый момент проекция всех параллелепипедов на эту плоскость предствляет собой множество прямоугольников и необходимо уметь обрабатывать следующие три запроса:
# добавить прямоугольник
# удалить прямоугольник
# посчитать площадь объединения всех прямоугольников

{{Определение
|definition=КД-деревом называется [[Дерево поиска, наивная реализация|сбалансированное двоичное дерево]], где каждой вершине <tex>\alpha</tex> соответствует некоторая область <tex>C_{\alpha}</tex> <tex>d</tex>-мерного пространстве, удовлетворяющая следующим свойствам:
# Корневой вершине <tex>root</tex> соответствует <tex>C_{root}</tex> - все пространство.
# Каждое <tex>C_{\alpha}</tex> является гиперпараллелепипедом (возможно бесконечным в некоторых направлениях).
# Для каждой вершины <tex>\alpha</tex> с детьми <tex>\beta</tex> и <tex>\gamma</tex> верно, что области <tex>C_{\beta}</tex> и <tex>C_{\gamma}</tex> не пересекаются и <tex>C_{\alpha} = C_{\beta} \cup C_{\gamma}</tex>.
}}

Будем хранить прямоугольники в КД-дереве следующим образом:
# Для каждого листа <tex>\alpha</tex> будем хранить все прямоугольники, которые пересекают внутреннюю часть области <tex>C_{\alpha}</tex>.
# Для каждой внутренней вершины <tex>\alpha</tex> будем считать число <tex>Cover_{\alpha}</tex> - количество прямоугольников, которые полностью покрывают область <tex>C_{\alpha}</tex>, но не полностью покрывают область родителя этой вершины <tex>C_{father(\alpha)}</tex>.

Также для каждой вершины будем хранить значение площади пересечения области этой вершины со всеми прямоугольниками <tex>M</tex>, которое определяется следующим образом:
# Для каждого листа <tex>M</tex> равно площади пересечения всех прямоугольников с областью этого листа.
# Для каждой внутренней вершины <tex>\alpha</tex> если <tex>Cover_{\alpha} > 0</tex>, тогда <tex>M</tex> равно площади всех области <tex>C_{\alpha}</tex>, иначе <tex>M</tex> равно сумме этих значений для двух сыновей этой вершины.

Таким образом, значение для корневой вершины <tex>M_{root}</tex> и является искомой площадью объединения всех прямоугольников.

Рассмотрим операцию добавления прямоугольника (box) в КД-дерево.

procedure Insert(box, <tex>\alpha</tex>)
if <tex>\alpha</tex> - лист then
добавить box в эту вершину
пересчитать <tex>M_{\alpha}</tex>
elseif box покрывает область <tex>C_{\alpha}</tex> then
<tex>Cover_{\alpha} = Cover_{\alpha} + 1</tex>
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
elseif box пересекает область <tex>C_{\alpha}</tex> then
Insert(box, leftson(<tex>\alpha</tex>))
Insert(box, rightson(<tex>\alpha</tex>))
if <tex>Cover_{\alpha} > 0 </tex> then
<tex>M_{\alpha} = </tex> объем области <tex>C_{\alpha}</tex>
else
<tex>M_{\alpha} = M_{leftson(\alpha)} + M_{rightson(\alpha)}</tex>

Аналогично устроена операция удаления прямоугольника - она выполняет обратные действия.

Рассмотрим теперь как именно будет устроено разделение пространства на области КД-деревом для эффективного его использования.
Возьмем проекции всех параллелепипедов на плоскость, перепендикулярную оси <tex>OZ</tex> - множество прямоугольников <tex>V</tex>. Рассмотрим все <tex>x</tex>-координаты границ прямоугольников и разобьем ось <tex>OX</tex> на <tex>2\sqrt{n}</tex> промежутков так, чтобы в каждом из них было не более <tex>\sqrt{n}</tex> границ. Получившиемя <tex>2\sqrt{n}</tex> полос разобьем горизонтальными линиями на итоговые области по отдельности следующим образом. Для каждой полосы у всех прямоугольников, имеющих вертикальную границу в этой полосе возьмем все горизонтальные границы (их будет не более <tex>2\sqrt{n}</tex>), а среди всех горизонтальных границ, пересекающих область, возьмем каждую <tex>\sqrt{n}</tex>-ую. Таким образом, в итоге получим не более <tex>4\sqrt{n}</tex> областей в каждой верикальной полосе.

Полученные области обладают следующими свойствами:
# Всего <tex>O(n)</tex> областей.
# Каждый прямоугольник частично покрывает не более <tex>O(\sqrt{n})</tex> областей.
# Никакая вершина прямоугольника не лежит внутри какой-нибудь области.
# Для каждой области есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти утверждения несложно следуют из алгоритма разбиения на области.

Каждая из этих областей будет одним из листов КД-дерева. Для построения всего дерева сначала объединим все области для одной вертикальной полосы, сливая соседние области в одну вершину. После этого аналогично будем сливать соседние вертикальные области и в итоге получим корневую вершину, представляющую всю плоскость.

Построенное дерево обладает следующими свойствами:
# В нем <tex>O(n)</tex> областей.
# Каждый прямоугольник хранится в <tex>O(\sqrt{n})</tex> листьях.
# Каждый прямоугольник влияет на <tex>O(\sqrt{n} \log n)</tex> значений <tex>Cover</tex>.
# Никакая вершина прямоугольника не лежит внутри области какого-нибудь листа.
# Для области каждого листа есть не более <tex>O(\sqrt{n})</tex> частично покрывающих ее прямоугольников.

Все эти свойства следуют из свойств разбиения на области и логарифмической высоты дерева.

Наконец, получаем итоговый алгоритм. Идем сканирующей плоскостью, перепендикулярной оси <tex>OZ</tex>, добавляя или удаляя прямоугольники из КД-дерева. При этом каждый раз пересчитывается площадь объединения всех текущих прямоугольников и, складывая эти величины, умноженные на соответствующую разницу <tex>z</tex>-координат получаем общий объем.

Каждая операция добавления или удаления прямоугольников требует изменения <tex>O(\sqrt{n} \log n)</tex> полей <tex>Cover</tex>, а также операций добавления или удаления прямоугольника в <tex>O(\sqrt{n})</tex> листьев, каждая из которых требует <tex>O(\log n)</tex> времени. Итого, на одну операцию добавления или удаления требуется <tex>O(\sqrt{n} \log n)</tex> времени, суммарно таких операций <tex>O(n)</tex>, и получается, что весь алгоритм работает за <tex>O(n \sqrt{n} \log n)</tex>. При этом на хранение всей информации требуется <tex>O(n\sqrt{n}) </tex> памяти - каждый лист может хранить <tex>O(\sqrt{n})</tex> прямоугольников.

Алгоритмы точного вычисления гиперобъема

2012-06-17T20:35:16Z

Joshik:

Алгоритмы точного вычисления гиперобъема

2012-06-17T18:54:57Z

Joshik:

Алгоритмы точного вычисления гиперобъема

2012-06-17T18:15:08Z

Joshik:

Алгоритмы точного вычисления гиперобъема

2012-06-17T17:49:16Z

Joshik:

Алгоритмы точного вычисления гиперобъема

2012-06-17T17:40:26Z

Joshik:

Алгоритмы точного вычисления гиперобъема

2012-06-17T17:20:49Z

Joshik:

Алгоритмы точного вычисления гиперобъема

2012-06-17T16:59:19Z

Joshik: Новая страница: «{{В разработке}} == Постановка задачи == <tex>x = (x_1, x_2, ..., x_d; x_i \ge 0) \in R^d</tex> - точка в <tex>d</tex>-мерн...»

Класс PCP

2010-06-01T12:42:54Z

Joshik:

==Определение==
Классом '''PCP[r(n), q(n)]''' ('''PCP''' - Probabilistically Checkable Proof), где <tex>n</tex> - длина входного слова, называется множество языков, распознаваемых машиной <tex>V^{\pi}(x)</tex>, обладающей следующими свойствами:

# Время работы <tex>V^{\pi}(x)</tex> ограничено сверху некоторым полиномом от длины <tex>x</tex>.
# <tex>\pi</tex> - некоторая строка, выступающая в качестве средства доказательства (аналогично P в [[Класс IP|интерактивном протоколе доказательства]]). Очевидно, ее длина не превосходит 2<sup>poly(x)</sup>, так как только к такому множеству позиций сможет обратиться <tex>V</tex>.
# <tex>V</tex> - [[Вероятностная машина Тьюринга|вероятностная машина Тьюринга]], обращающаяся к случайной ленте не более <tex>r(n)</tex> раз.
# <tex>V</tex> обращается к строке <tex>\pi</tex> не более <tex>q(n)</tex> раз.
# <tex>x \in L \Rightarrow \exists \pi : Pr(V^{\pi}(x)=1)=1 \ </tex> , где <tex>Pr(V^{\pi}(x)=1)</tex> - вероятность того, что <tex>V</tex> допустит <tex>x</tex>.
# <tex>\ x \notin L \Rightarrow \forall \pi : Pr(V^{\pi}(x)=1)\le \frac{1}{2} </tex>.

==Свойства==

# '''PCP[0, 0] = [[P|P]]''' (по определению '''P''' - нет случайности и обращений к <tex>\pi</tex>)
# '''PCP[log(n), 0] = [[P|P]]''' (логарифмическое число обращений к случайной ленте не помогают, так как можно за полиномиальное время перебрать всевозможные результаты обращений)
# '''PCP[0, log(n)] = [[P|P]]''' (логарифмическое число обращений к строке <tex>\pi</tex> также не помогают, так как можно аналогичным образом перебрать всевозможные результаты обращений за полиномиальное время)
# '''PCP[poly(n), 0] = [[Сложностные_классы_RP_и_coRP|coRP]]''' (по определению '''coRP''')
# '''PCP[0, poly(n)] = [[NP|NP]]''' (по определению '''NP''' на языке сертификатов)
# '''PCP[log(n), O(1)] = [[NP|NP]]''' ([[PCP-теорема]])

==Теорема==

'''[[GNI]]''' ∈ '''PCP[poly(n), O(1)]'''

==Доказательство==

Алгоритм работы машины '''V''' аналогичен алгоритму работы при доказательстве того, что '''[[GNI]]''' ∈ '''IP'''.

В строке <tex>\pi</tex> для каждого графа (введена некоторая нумерация графов) записано, кому из данных графов <tex>G_1, G_2</tex> он изоморфен.

Класс PCP

2010-06-01T12:25:54Z

Joshik: /* Свойства */

==Определение==
Классом '''PCP[r(n), q(n)]''' ('''PCP''' - Probabilistically Checkable Proof), где <tex>n</tex> - длина входного слова, называется множество языков, распознаваемых машиной <tex>V^{\pi}(x)</tex>, обладающей следующими свойствами:

1) Время работы <tex>V^{\pi}(x)</tex> ограничено сверху некоторым полиномом от длины <tex>x</tex>

2) <tex>\pi</tex> - некоторая строка, выступающая в качестве средства доказательства (аналогично P в [[Класс IP|интерактивном протоколе доказательства]]). Очевидно, ее длина не превосходит 2<sup>poly(x)</sup>, так как только к такому множеству позиций сможет обратиться <tex>V</tex>

3) <tex>V</tex> - [[Вероятностная машина Тьюринга|вероятностная машина Тьюринга]], обращающаяся к случайной ленте не более <tex>r(n)</tex> раз

4) <tex>V</tex> обращается к строке <tex>\pi</tex> не более <tex>q(n)</tex> раз

5) <tex>x \in L \Rightarrow \exists \pi : Pr(V^{\pi}(x)=1)=1 \ </tex> , где <tex>Pr(V^{\pi}(x)=1)</tex> - вероятность того, что <tex>V</tex> допустит <tex>x</tex>

6) <tex>\ x \notin L \Rightarrow \forall \pi : Pr(V^{\pi}(x)=1)\le \frac{1}{2} </tex>

==Свойства==

1) '''PCP[0, 0] = [[P|P]]''' (по определению '''P''' - нет случайности и обращений к <tex>\pi</tex>)

2) '''PCP[log(n), 0] = [[P|P]]''' (логарифмическое число обращений к случайной ленте не помогают, так как можно за полиномиальное время перебрать всевозможные результаты обращений)

3) '''PCP[0, log(n)] = [[P|P]]''' (логарифмическое число обращений к строке <tex>\pi</tex> также не помогают, так как можно аналогичным образом перебрать всевозможные результаты обращений за полиномиальное время)

4) '''PCP[poly(n), 0] = [[Сложностные_классы_RP_и_coRP|coRP]]''' (по определению '''coRP''')

5) '''PCP[0, poly(n)] = [[NP|NP]]''' (по определению '''NP''' на языке сертификатов)

6) '''PCP[log(n), O(1)] = [[NP|NP]]''' ([[PCP-теорема]])

Теория сложности (старая трешовая версия)

2010-05-30T09:57:45Z

Joshik:

== Лекция 1 ==
*[[Класс DSPACE]]
*[[Класс DTIME]]
*[[Теорема о емкостной иерархии]]
*[[Теорема о временной иерархии]]
*[[Сведение по Карпу]]
*[[Сведение по Куку]]
*[[Класс P]]
*[[Класс NP]]
*[[Класс coNP]]

== Практика 1 ==
*[[Сведение по Куку задачи факторизации к языку из NP]]

== Лекция 2 ==
*[[Теорема Кука]]

== Практика 2 ==
*[[Понятие NP-трудной и NP-полной задачи]]
*[[NP-полнота задачи BH1N]]
*[[NP-полнота задачи о выполнимости булевой формулы в форме КНФ]]
*[[NP-полнота задачи о выполнимости булевой формулы в форме 3-КНФ]]
*[[NP-полнота задачи о клике]]
*[[NP-полнота задачи о независимом множестве]]
*[[NP-полнота задачи о вершинном покрытии]]

== Лекция 3 ==
*[[Теорема Ладнера]]
*[[Теорема Левина]]
*[[Теорема Бейкера-Гилла-Соловэя]]

== Практика 3 ==
*[[NP-полнота задач о гамильтоновом цикле и пути в графах]]
*[[NP-полнота задачи о сумме подмножества]]
*[[NP-полнота задачи о рюкзаке]]

== Практика, которой на самом деле не было ==
*[[NP-полнота задачи о раскраске графа]]

== Лекция 4 ==
*[[Класс PS]]
*[[Теорема Сэвича]]
*[[PS-полнота задачи Generalized geography]]

== Лекция 6 ==
*Классы [[L]], [[NL]], [[NL-полнота|NLC]]
*[[NL-полнота задачи о достижимости в графе]]
*[[Классы EXP, NEXP. Полнота языков EXP и NEXP]]
*[[Теорема о связи вопросов EXP=NEXP и P=NP]]
*[[Теорема Иммермана]]

== Практика 6 ==
*[[Классы Sigma_i и Pi_i]]
*[[Класс PH]]
*[[Полиномиальная иерархия]]
*[[Теоремы о коллапсе полиномиальной иерархии]]

== Лекция 7 ==
*[[Теорема Карпа-Липтона]]

== Практика 7 ==
*[[Вероятностная машина Тьюринга]]
*[[Класс ZPP]]
*[[Сложностные классы RP и coRP]]
*[[Сложностный класс PP]]
*[[Сложностный класс BPP]]
*[[Уменьшение ошибки в классе RP, сильное и слабое определение]]

== Лекция 8 ==
*[[Теорема о включении BPP в P/poly]]
*[[Теорема Лаутемана]]
*[[Теорема Валианта-Вазирани]]

== Практика 8 ==
*[[Лемма Шварца-Зиппеля]]

== Лекция 9 ==
*[[Класс IP|Класс IP]]
*[[GNI|Принадлежность проблемы GNI классу IP]]
*[[Sharp SAT|#SAT]]

== Лекция 10 ==
*[[Теорема Шамира]]
*[[Семейство универсальных попарно независимых хеш-функций|Семейство универсальных попарно независимых хеш-функций]]
*[[Теорема Голдвассера, Сипсера]]

== Лекция 11 ==
*[[Абсолютная секретность]]
*[[Лемма о невозможности существования вычислительно безопасных шифров в случае P = NP]]
*[[Односторонние функции и псевдослучайные генераторы]]
*[[Доказательства с нулевым разглашением]]

== Лекция 12 ==
*[[Кубит]]
*[[Унитарные операторы]]
*[[Квантовый логический элемент NOT]]
*[[Преобразование Адамара]]
*[[Квантовый логический элемент CNOT]]
*[[Квантовый логический элемент Тоффоли]]
*[[Квантовая схема]]

== Лекция 13 ==
*[[Класс PCP]]

NL-полнота задачи о достижимости в графе

2010-04-06T11:48:20Z

Joshik:

== Формулировка задачи==
Даны ориентированный граф <tex> G = \langle V, E \rangle </tex> и две вершины <tex> s, t</tex> в нем. Необходимо проверить, правда ли, что в графе <tex> G </tex> существует путь из вершины <tex> s </tex> в вершину <tex> t </tex>. Эту задачу принято называть <tex> st-connectivity </tex> или <tex> STCON </tex>.

== Утверждение ==
Задача <tex> STCON </tex> [[NL-полнота|NL-полна]].

== Доказательство ==
Для доказательства [[NL-полнота|NL-полноты]] необходимо показать, что эта задача NL-трудная и принадлежит [[NL|классу NL]].

=== Доказательство принадлежности задачи STCON классу NL ===

Для доказательства необходимо предъявить алгоритм для недетерминированной машины Тьюринга, который использует конечное число переменных, каждая из которых занимает <tex> O(log n) </tex> памяти, где <tex> n </tex> - размер входа для задачи и за время порядка <tex> O(poly(n)) </tex> решает эту задачу.

Алгоритм:

1. Начиная с вершины <tex> s </tex> недетерминированно переходит в одну из вершин, смежных с ней. (Очевидно, для этого необходимо конечное число переменных)

2. Проверяет, правда ли, что текущая вершина совпадает с <tex> t </tex>. Если это так, возвращает TRUE.

3. Отдельно считает количество пройденных вершин. Как только это число превышает количество вершин в графе, то алгоритм возвращает FALSE, так как посетил некоторую вершину дважды.

Таким образом в каждый момент алгоритму достаточно хранить текущую вершину, количество посещенных вершин, финальную вершину <tex> t </tex> и некоторое число вспомогательных переменных, для совершения переходов. Все эти переменные принимают значения не более, чем максимальный номер вершины, то есть как раз занимают <tex> O(log n) </tex> памяти.

Так как, если из <tex> s </tex> существует путь в <tex> t </tex>, то он имеет длину не более, чем количество вершин в графе, то алгоритм корректно возвращает FALSE.

=== Доказательство NL-трудности задачи STCON ===

NL-полнота задачи о достижимости в графе

2010-04-06T10:58:06Z

Joshik:

NL-полнота задачи о достижимости в графе

2010-04-06T10:29:04Z

Joshik: Новая страница: «== Формулировка задачи== Даны ориентированный граф <tex> G = \langle V, E \rangle </tex> и две вершины <tex> s, t…»

Теория сложности (старая трешовая версия)

2010-04-06T10:24:21Z

Joshik:

Теорема о ёмкостной иерархии

2010-03-23T14:35:00Z

Joshik:

== Формулировка ==
'''Теорема о емкостной иерархии''' утверждает, что для любых двух [[Конструируемая по памяти функция|конструируемых по памяти функций]] <tex>f</tex> и <tex>g</tex> таких, что <tex> \lim \limits_{n \rightarrow \infty} f(n)/g(n) = 0</tex>, выполняется <tex>DSPACE(g(n)) \ne DSPACE(f(n))</tex>.

== Доказательство ==
Зафиксируем функции <tex>f</tex> и <tex>g</tex>.

Рассмотрим язык <tex>L = \{ \langle m,x \rangle \mid m(\langle m,x \rangle )</tex> не допускает, используя не более <tex> f(|\langle m,x\rangle|)</tex> памяти <tex>\}</tex> и докажем, что <tex>L \notin DSPACE(f)</tex> и <tex>L \in DSPACE(g)</tex>.

Допустим, что <tex>L \in DSPACE(f)</tex>, тогда существует детерминированная машина Тьюринга <tex>m_0</tex> такая, что <tex>L(m_0)=L</tex>.

Рассмотрим выход машины <tex>m_0</tex> на входе <tex>\langle m_0,x\rangle</tex>.

Пусть <tex>m_0</tex> допускает <tex>\langle m_0,x\rangle</tex>. Тогда <tex>\langle m_0,x\rangle \in L</tex>, но в <tex>L</tex> по определению не может быть пары <tex>\langle m,x\rangle</tex>, которую допускает <tex>m</tex>. Таким образом, <tex>m_0</tex> не может допускать <tex>\langle m_0,x\rangle</tex>.

Если <tex>m_0</tex> не допускает <tex>\langle m_0,x\rangle</tex>, то <tex>\langle m_0,x\rangle</tex> не принадлежит языку <tex>L</tex>. Из определения это значит, что либо <tex>m_0</tex> допускает <tex>\langle m_0,x\rangle</tex>, либо не допускает, используя памяти больше <tex>f(|\langle m_0,x\rangle|)</tex>. Но <tex>m_0</tex> выбрана таким образом, что на любом входе <tex>x</tex> она использует не более <tex>f(|x|)</tex> памяти. Получаем противоречие.

Следовательно, такой машины не существует. Таким образом, <tex>L \notin DSPACE(f)</tex>.

<tex>L \in DSPACE(g)</tex>, так как языку <tex>L</tex> можно сопоставить машину Тьюринга <tex>m_0</tex>, распознающую <tex>L</tex> и такую, что на любом входе <tex>\langle m_1,x\rangle \in L</tex> <tex>m_0</tex> будет работать аналогично <tex>m_1</tex>. Если <tex>m_1</tex> завершила работу, используя не более <tex>f(|\langle m_1,x\rangle|)</tex> памяти, и не допустила, то <tex>m_0</tex> допускает <tex>\langle m_1,x\rangle</tex>. В другом случае не допускает. Любая такая машина использует памяти не более <tex>f(|\langle m_1,x\rangle|)</tex>. По условию теоремы <tex>\lim \limits_{n \rightarrow \infty} f(n)/g(n) = 0</tex>, поэтому начиная с некоторого <tex>n</tex>, <tex>m_1</tex> будет использовать памяти не более <tex>g(|\langle m_1,x\rangle|)</tex>.

Таким образом получили, что <tex>L \in DSPACE(g(n)) \setminus DSPACE(f(n))</tex>. Следовательно, <tex>DSPACE(g(n)) \neq DSPACE(f(n))</tex>, что и требовалось доказать.