Об ускорении поиска повторов ОСАМ-методом — различия между версиями

Материал из YourcmcWiki
Перейти к: навигация, поиск
м
м (Дальнейшие изыскания)
 
(не показано 17 промежуточных версий этого же участника)
Строка 5: Строка 5:
 
Имеем следующий алгоритм поиска приближённых повторов в ДНК-последовательностях:
 
Имеем следующий алгоритм поиска приближённых повторов в ДНК-последовательностях:
  
* Скользим по последовательности окном размера W и каждую букву заменяем на количество пуринов (A, G, U — аденин, гуанин, урацил) / пиримидинов (C, T — цитозин, тимин) в этом окне (W-окрестности буквы). Полученную числовую последовательность называем «профилем последовательности».
+
* Скользим по последовательности окном размера W и каждую букву заменяем на количество пуринов (A, G, U — аденин, гуанин, урацил) / пиримидинов (C, T — цитозин, тимин) в этом окне (W-окрестности буквы). Полученную числовую последовательность называем «профилем последовательности». Здесь могут быть вариации — скажем, можно считать число других пар букв — например, G и C.
 
* Дальше, начиная с каждой из (N-A)/S позиций, где N — длина последовательности, A — второе окно, «окно аппроксимации», S — шаг аппроксимации (расстояние между соседними начальными позициями), выбираем A чисел. Приняв полученный вектор длины A за функцию дискретного аргумента, осуществляем дискретное спектральное преобразование этой функции в соответствии с одним из классических ортогональных базисов (например, с базисом Фурье), и вычисляем K спектральных коэффициентов. Преобразование может осуществляться как на сетке размера А, так и на сетке меньшего размера. Учитывая то, что каждое значение «включает в себя» информацию о соседних, так как является скользящей суммой — уменьшение сетки влияет на результат очень незначительно. Получаем в каждой из (N-A)/S позиций K коэффициентов, характеризующих часть последовательности длиной A.
 
* Дальше, начиная с каждой из (N-A)/S позиций, где N — длина последовательности, A — второе окно, «окно аппроксимации», S — шаг аппроксимации (расстояние между соседними начальными позициями), выбираем A чисел. Приняв полученный вектор длины A за функцию дискретного аргумента, осуществляем дискретное спектральное преобразование этой функции в соответствии с одним из классических ортогональных базисов (например, с базисом Фурье), и вычисляем K спектральных коэффициентов. Преобразование может осуществляться как на сетке размера А, так и на сетке меньшего размера. Учитывая то, что каждое значение «включает в себя» информацию о соседних, так как является скользящей суммой — уменьшение сетки влияет на результат очень незначительно. Получаем в каждой из (N-A)/S позиций K коэффициентов, характеризующих часть последовательности длиной A.
 
* Получаем вместо длинной строки кучу векторов по K коэффициентов. Проделываем это с обеими сравниваемыми последовательностями (или с одной, если сравниваем её саму с собой).
 
* Получаем вместо длинной строки кучу векторов по K коэффициентов. Проделываем это с обеими сравниваемыми последовательностями (или с одной, если сравниваем её саму с собой).
* Сравниваем все вектора попарно с помощью какой-нибудь метрики (например, евклидовой, максимума модуля разности или суммы модулей разностей координат) и если два вектора оказываются ближе друг к другу, чем на заданное число ε, считаем, что области последовательностей, соответствующие этим векторам, похожи друг на друга.
+
* Сравниваем все вектора попарно с помощью какой-нибудь метрики (например, евклидовой, максимума модуля разности или суммы модулей разностей координат) и если два вектора оказываются ближе друг к другу, чем на заданное число ε, считаем, что куски последовательностей длины A, соответствующие этим векторам, похожи друг на друга.
 +
* Визуально матрица сравнений отображается обычным изображением, в котором точки одного цвета соответствуют «близким» векторам, а точки другого цвета — отсутствию близости.
 +
* Далее можно заметить, что более длинные повторы в последовательностях соответствуют отрезкам, параллельным главной диагонали. Соответственно, найдя все отрезки заданной длины на матрице сравнений, можно найти приблизительные координаты длинных прямых повторов.
 +
* Любопытным свойством алгоритма является то, что на основе уже сгенерированных спектральных векторов мы можем найти также и обратные повторы. Ибо смена знака всех коэффициентов разложения с чётными номерами (если считать с 1) соответствует разложению f(-x), так как чётные коэффициенты — это либо полиномы нечётной степени, либо синусы, а и те, и те являются нечётными функциями. А f(-x) в случае использования сетки от −1 до 1 соответствует обращению строки, из которой сгенерирован раскладываемый вектор. Соответственно, для поиска обратных повторов нужно лишь немного поменять способ вычисления метрики и заметить, что более длинные обратные повторы соответствуют отрезкам, параллельным побочной диагонали матрицы.
  
Теперь представляем, что хотим это ускорить. Для простоты, за расстояние между векторами берём максимум модуля разности, а базис используем тригонометрический (Фурье). Хотим, по сути, уйти от сложности O(N²). Для этого можно применить или создать какую-нибудь индексную структуру.
+
Итак, в итоге основная сложность алгоритма складывается из сложности двух вещей — получения координат разложения и сравнения их друг с другом.
  
Но здесь нас поджидает две проблемы:
+
Теперь представляем, что алгоритм нужно ускорить. Для простоты, за расстояние между векторами берём максимум модуля разности, а базис используем тригонометрический (Фурье). Например, хотим уйти от сложности O(N²) (или N*M) при сравнении. Для этого можно применить или создать какую-нибудь индексную структуру.
 +
 
 +
Но здесь нас поджидают две проблемы:
  
 
==== Природа спектральных векторов ====
 
==== Природа спектральных векторов ====
  
Первая проблема — природа коэффициентов разложения. Их распределение похоже на нормальное, и что ещё хуже, первые коэффициенты сильно больше последних по модулю. Сильно — это примерно на 1-2 порядка. Причины две — во-первых, в реальных данных низкочастотные гармоники всегда больше, во-вторых, в нашем случае высокочастотные гармоники приглушены нами же — исходный вектор длины K => разница между двумя соседними значениями в нём по модулю не больше A/K, так как разница между двумя соседними значениями в профиле последовательности либо 0, либо 1, либо −1. Следовательно, максимум i-го коэффициента примерно (K+1-i)*A/K. И всё это — распределено околонормально. Таким образом, разбиение пространства на какие-то области (построение индекса) уже становится проблемой — «полезных» измерений с большим относительно размера окрестности поиска диапазоном мало (3-5), бОльшая часть спектральных векторов распределена в центре всей области, центр маленький. То есть, при поиске приходится просматривать довольно много лишних векторов.
+
Первая проблема — природа коэффициентов разложения. Они распределены весьма плотно по сравнению с типичным размером окрестности поиска, а распределение похоже на нормальное, и что ещё хуже, первые коэффициенты сильно больше последних по модулю. Сильно — это примерно на 1-2 порядка. Причины две — во-первых, в реальных данных низкочастотные гармоники всегда оказываются больше, во-вторых, в нашем случае высокочастотные гармоники приглушены нами же — разница между двумя соседними значениями в раскладываемом векторе значений не больше 1, если сетка полная (размер сетки = A), и не больше отношения длины окна A к размеру сетки A/L, если она меньше A, так как разница между двумя соседними значениями в профиле последовательности равна либо 0, либо 1, либо −1. Следовательно, максимум i-го коэффициента разложения равен примерно (K+1-i)*A/L. И всё это — распределено около-нормально (чему, видимо, тоже есть обоснование — там же везде суммы «случайных» значений, а закон больших чисел никто не отменял).
 +
 
 +
Таким образом, разбиение пространства на какие-то области (построение индекса) уже становится проблемой — «полезных» измерений с большим относительно размера окрестности поиска диапазоном мало (3-5), бОльшая часть спектральных векторов распределена в центре всей области, центр маленький. То есть, при поиске приходится просматривать довольно много лишних векторов.
  
 
==== Простота линейного поиска ====
 
==== Простота линейного поиска ====
Строка 22: Строка 29:
 
Линейный поиск очень прост по сравнению с поиском по индексу — сложность одного шага поиска по индексу сама по себе гораздо больше сложности одного шага «тупого» поиска, так как, скорее всего, включает ветвления, переходы по дереву, запоминание списка найденных элементов и тому подобные операции. Поэтому, даже выигрывая в числе шагов, мы получим меньшее ускорение, чем хотели бы.
 
Линейный поиск очень прост по сравнению с поиском по индексу — сложность одного шага поиска по индексу сама по себе гораздо больше сложности одного шага «тупого» поиска, так как, скорее всего, включает ветвления, переходы по дереву, запоминание списка найденных элементов и тому подобные операции. Поэтому, даже выигрывая в числе шагов, мы получим меньшее ускорение, чем хотели бы.
  
Кроме того, с очевидной оптимизацией линейный поиск становится ещё проще. При попарном сравнении всех векторов не нужно делать их сравнивать полностью! Их достаточно сравнивать с искомой окрестностью последовательно от первой координаты к последней K-ой. Учитывая природу коэффициентов разложения — это примерно соответствует порядку убывания абсолютных значений коэффициентов, а значит, в среднем и разности.
+
Кроме того, с очевидной оптимизацией линейный поиск становится ещё проще. При попарном сравнении всех векторов не нужно их сравнивать полностью! Их достаточно сравнивать с искомой окрестностью последовательно от первой координаты к последней K-ой. Учитывая природу коэффициентов разложения — это примерно соответствует порядку убывания абсолютных значений коэффициентов, а значит, в среднем и разности.
  
 
То есть, большая часть неподходящих нам векторов отсекается вообще по 1-ой координате. Относительно мало — по второй, ещё меньше — по третьей, и т. п. В реальности получается, что по коэффициентам дальше 6-го не отсекается почти ничего, а дальше 9-го — совсем ничего. В то же время, для полной проверки ''подходящего'' вектора нужно сделать в точности K сравнений и не меньше.
 
То есть, большая часть неподходящих нам векторов отсекается вообще по 1-ой координате. Относительно мало — по второй, ещё меньше — по третьей, и т. п. В реальности получается, что по коэффициентам дальше 6-го не отсекается почти ничего, а дальше 9-го — совсем ничего. В то же время, для полной проверки ''подходящего'' вектора нужно сделать в точности K сравнений и не меньше.
Строка 32: Строка 39:
 
Из-за того, что по высшим гармоникам фильтрации почти не происходит, целесообразно оказывается использовать 3, максимум 5 коэффициентов разложения. В простейшем случае эти 3 коэффициента — среднее, проекция на синус, проекция на косинус.
 
Из-за того, что по высшим гармоникам фильтрации почти не происходит, целесообразно оказывается использовать 3, максимум 5 коэффициентов разложения. В простейшем случае эти 3 коэффициента — среднее, проекция на синус, проекция на косинус.
  
== K-D-N дерево ==
+
== Индексные структуры ==
  
 
В простейшей реализации алгоритма поиска повторов все сгенерированные вектора спектральных координат сравниваются друг с другом попарно, что сильно ухудшает производительность. Порядок сложности такого алгоритма — O(N²). Чтобы его ускорить, нужно применять индексные структуры, позволяющие уменьшить сложность поиска.
 
В простейшей реализации алгоритма поиска повторов все сгенерированные вектора спектральных координат сравниваются друг с другом попарно, что сильно ухудшает производительность. Порядок сложности такого алгоритма — O(N²). Чтобы его ускорить, нужно применять индексные структуры, позволяющие уменьшить сложность поиска.
  
Обычное K-D (от K-Dimensional) дерево представляет собой вид BSP-дерева (BSP = Binary Space Partitioning, разбиение пространства на две части), в котором все разбиения осуществляются гиперплоскостями, перпендикулярными одной из осей координат. Каждый узел K-D дерева характеризуется номером координаты и пороговым значением. Если значение выбранной координаты в точке меньше или равно этого порогового значения, точка относится к левому подпространству, а если больше — к правому.
+
Для индексации многомерных данных придумано очень много разных структур и алгоритмов. Большинство из них, правда, осталось на теоретическом уровне статей и диссертаций :) в реальности обычно используются самые простые вариации. В целом, все многомерные индексы делятся на два класса:
 +
* Производные от R-дерева. R-дерево — наиболее известная и используемая в реальности индексная структура (R — от Region), прямое обобщение B-деревьев на многомерный случай. Узлы индекса — «гиперпрямоугольники» многомерного пространства, которые могут содержать либо точки, либо дочерние узлы. R-дерево эффективно работает для минимальных размерностей (2-3) и поэтому широко используется для хранения пространственных данных, например, положений объектов на картах местности. При увеличении размерности, однако, возникает проблема перекрытия узлов и индекс перестаёт давать ускорение, так как при поиске оказывается нужно заглянуть в почти все узлы. Соответственно, большинство производных от R-дерева структур ставят своей целью именно борьбу с перекрытиями узлов.
 +
* BSP-деревья и их обобщения (BSP = Binary Space Partitioning). Каждый узел представляет собой разбиение пространства на две части, таким образом, части пространства никогда не пересекаются.
 +
 
 +
Существуют индексы, разрабатываемые специально для использования в СУБД — обычно от таких индексов требуется ориентация на дисковое хранение и обновляемость без значительной потери эффективности. Существуют и более простые структуры, предназначенные для использования в памяти и не имеющие возможности обновления. На самом деле, из реально широко используемых обновляемых многомерных индексов можно упомянуть только R-дерево.
 +
 
 +
Существует два основных применения индексов для поиска:
 +
* Поиск ближайшего соседа (точки, наиболее близкой к заданной) или нескольких ближайших соседей. Здесь наиболее эффективны производные R-дерева, так как их узлы содержат границы областей, для каждого узла можно вычислить максимально возможное расстояние от искомого элемента, и соответственно, узлы можно просматривать сразу в нужном порядке.
 +
* Поиск по заданной области, требуемый как раз в нашем случае. В нашем случае не требуется ни внешнее хранение, ни обновляемость индекса (с ней у R-деревьев получше), поэтому более эффективными оказываются разделяющие деревья.
 +
 
 +
=== K-D дерево ===
 +
 
 +
Обычное K-D (от K-Dimensional) дерево представляет собой вид BSP-дерева, в котором все разбиения осуществляются гиперплоскостями, перпендикулярными одной из осей координат. Каждый узел K-D дерева характеризуется номером координаты и пороговым значением. Если значение выбранной координаты в точке меньше или равно этого порогового значения, точка относится к левому подпространству, а если больше — к правому.
 +
 
 +
Выбор разделяющих плоскостей в простейшей и наиболее часто используемой реализации K-D дерева делается просто — выбирается медианное значение нужной координаты по всей выборке индексируемых точек. Таким образом строится сбалансированное K-D дерево. Выбор измерения для разбиения пространства на каждом шаге в простейшем случае производится по кругу, начиная с первого измерения — такой подход даёт наибольшее «равноправие» измерений.
 +
 
 +
== K-D-N дерево ==
 +
 
 +
K-D-N дерево:
 +
* Содержит уровней не больше, чем число измерений.
 +
* Нелистовые узлы содержат полное разбиение набора данных на диапазоны по одному из измерений, плюс вычисленные на этапе построения дерева минимальное и максимальное значение выбранного измерения у всех попавших в узел элементов. Размеры всех диапазонов одинаковы и равны половине размера ожидаемой окрестности, по которой будет производиться поиск. Больше хуже — массив данных очень плотный, и разделение пространства становится неэффективным. Меньше тоже хуже — сильно возрастает количество узлов дерева, которые нужно просматривать при поиске.
 +
* Если полный диапазон по какому-то измерению меньше, чем 1.5 * размер ожидаемой при запросе окрестности поиска, разбиение по этому измерению не производится, потому что оно всё равно будет неэффективно — при поиске почти наверняка придётся заглядывать во все дочерние узлы. Почему 1.5? Потому что это 3*ε, если ищем по окрестности (c-ε ≤ x ≤ c+ε). В этом случае пространство разделится на 3 «полосы» каждая размером по ε, запрашиваемый центр окрестности почти наверняка попадёт в центральную из этих полос (помним, что распределение у нас почти нормальное), соответственно, искомые элементы почти наверняка смогут оказаться в любой из трёх полос.
 +
* Листовые узлы содержат наборы элементов. Технически — списки индексов элементов в изначальном массиве данных. Сами элементы хранить смысла нет, так как нам всё равно нужны индексы, чтобы сопоставлять элементам положения в ДНК-последовательности.
 +
* Дополнительно создаются и хранятся вариационные ряды (списки индексов элементов), упорядоченные по каждой координате.
 +
* При поиске просто спускаемся по дереву, просматривая узлы, в которых могут содержаться подходящие элементы, и всё время запоминаем максимальное и минимальное значение по каждому измерению. Так как индексируемый набор данных известен полностью, запоминаются они в виде номера элемента в вариационном ряде — это позволяет избежать лишних операций с плавающей точкой.
 +
* После поиска нужно произвести окончательную фильтрацию неточного найденного набора элементов. Фильтрация по каждому измерению производится одним из двух способов, оптимальный из которых выбирается сравнением размера фильтруемого набора с количеством элементов, которые придётся отсечь на шаге 2:
 +
** Либо простым сравнением значения каждого элемента с искомой окрестностью — эффективно, если фильтруемый набор мал.
 +
** Либо отсечением кусков вариационного ряда по ассоциативному массиву. Для этого двоичным поиском находится положение минимального и максимального значения измерения искомой окрестности (c-ε, c+ε) в вариационном ряде по нужной координате, а потом участки вариационного ряда от «известного» до «найденного» максимального/минимального значения удаляются из результата. Под «известным» диапазоном понимается найденный на этапе обхода дерева. Технически — во временный массив по номерам элементов из этих участков вариационного ряда сохраняется «1», а фильтруемый набор одним проходом фильтруется с помощью проверки элементов временного массива. На самом деле, при многократном поиске временный массив можно не обнулять, а вместо «1» можно использовать увеличиваемый при каждом следующем поиска счётчик. Эта стратегия эффективна, если фильтруемый набор больше, чем сумма размеров удаляемых участков. Обычно это случается, либо если какое-то измерение почти ничего не фильтрует (то есть, если полный диапазон значений по этому измерению сопоставим с размером окрестности поиска), либо если просто найдено очень много элементов.
 +
* С помощью этих мер достигается эффективность индекса в ситуации, когда пространственный размер массива данных относительно маленький, значимых измерений относительно немного, но при этом нет ограничения на полное число измерений.
 +
 
 +
== Дальнейшие изыскания ==
 +
 
 +
Новая идея, но это будет что-то типа гибрида KD- и M-дерева (K-D-M-дерево? :)). Т.е. юзать то же самое кд-дерево, но очень близкие элементы группировать в гиперкубики и индексировать их центры вместо самих элементов.
  
Выбор разделяющих плоскостей в простейшей и наиболее часто используемой реализации K-D дерева делается просто — выбирается медианное значение нужной координаты по всей выборке индексируемых точек. Таким образом строится сбалансированное К-Д дерево.
+
Можно даже тотально перейти в целые числа.

Текущая версия на 02:08, 22 ноября 2011

О методах ускорения работы алгоритма поиска приближённых повторов в ДНК-последовательностях.

Алгоритм поиска повторов

Имеем следующий алгоритм поиска приближённых повторов в ДНК-последовательностях:

  • Скользим по последовательности окном размера W и каждую букву заменяем на количество пуринов (A, G, U — аденин, гуанин, урацил) / пиримидинов (C, T — цитозин, тимин) в этом окне (W-окрестности буквы). Полученную числовую последовательность называем «профилем последовательности». Здесь могут быть вариации — скажем, можно считать число других пар букв — например, G и C.
  • Дальше, начиная с каждой из (N-A)/S позиций, где N — длина последовательности, A — второе окно, «окно аппроксимации», S — шаг аппроксимации (расстояние между соседними начальными позициями), выбираем A чисел. Приняв полученный вектор длины A за функцию дискретного аргумента, осуществляем дискретное спектральное преобразование этой функции в соответствии с одним из классических ортогональных базисов (например, с базисом Фурье), и вычисляем K спектральных коэффициентов. Преобразование может осуществляться как на сетке размера А, так и на сетке меньшего размера. Учитывая то, что каждое значение «включает в себя» информацию о соседних, так как является скользящей суммой — уменьшение сетки влияет на результат очень незначительно. Получаем в каждой из (N-A)/S позиций K коэффициентов, характеризующих часть последовательности длиной A.
  • Получаем вместо длинной строки кучу векторов по K коэффициентов. Проделываем это с обеими сравниваемыми последовательностями (или с одной, если сравниваем её саму с собой).
  • Сравниваем все вектора попарно с помощью какой-нибудь метрики (например, евклидовой, максимума модуля разности или суммы модулей разностей координат) и если два вектора оказываются ближе друг к другу, чем на заданное число ε, считаем, что куски последовательностей длины A, соответствующие этим векторам, похожи друг на друга.
  • Визуально матрица сравнений отображается обычным изображением, в котором точки одного цвета соответствуют «близким» векторам, а точки другого цвета — отсутствию близости.
  • Далее можно заметить, что более длинные повторы в последовательностях соответствуют отрезкам, параллельным главной диагонали. Соответственно, найдя все отрезки заданной длины на матрице сравнений, можно найти приблизительные координаты длинных прямых повторов.
  • Любопытным свойством алгоритма является то, что на основе уже сгенерированных спектральных векторов мы можем найти также и обратные повторы. Ибо смена знака всех коэффициентов разложения с чётными номерами (если считать с 1) соответствует разложению f(-x), так как чётные коэффициенты — это либо полиномы нечётной степени, либо синусы, а и те, и те являются нечётными функциями. А f(-x) в случае использования сетки от −1 до 1 соответствует обращению строки, из которой сгенерирован раскладываемый вектор. Соответственно, для поиска обратных повторов нужно лишь немного поменять способ вычисления метрики и заметить, что более длинные обратные повторы соответствуют отрезкам, параллельным побочной диагонали матрицы.

Итак, в итоге основная сложность алгоритма складывается из сложности двух вещей — получения координат разложения и сравнения их друг с другом.

Теперь представляем, что алгоритм нужно ускорить. Для простоты, за расстояние между векторами берём максимум модуля разности, а базис используем тригонометрический (Фурье). Например, хотим уйти от сложности O(N²) (или N*M) при сравнении. Для этого можно применить или создать какую-нибудь индексную структуру.

Но здесь нас поджидают две проблемы:

Природа спектральных векторов

Первая проблема — природа коэффициентов разложения. Они распределены весьма плотно по сравнению с типичным размером окрестности поиска, а распределение похоже на нормальное, и что ещё хуже, первые коэффициенты сильно больше последних по модулю. Сильно — это примерно на 1-2 порядка. Причины две — во-первых, в реальных данных низкочастотные гармоники всегда оказываются больше, во-вторых, в нашем случае высокочастотные гармоники приглушены нами же — разница между двумя соседними значениями в раскладываемом векторе значений не больше 1, если сетка полная (размер сетки = A), и не больше отношения длины окна A к размеру сетки A/L, если она меньше A, так как разница между двумя соседними значениями в профиле последовательности равна либо 0, либо 1, либо −1. Следовательно, максимум i-го коэффициента разложения равен примерно (K+1-i)*A/L. И всё это — распределено около-нормально (чему, видимо, тоже есть обоснование — там же везде суммы «случайных» значений, а закон больших чисел никто не отменял).

Таким образом, разбиение пространства на какие-то области (построение индекса) уже становится проблемой — «полезных» измерений с большим относительно размера окрестности поиска диапазоном мало (3-5), бОльшая часть спектральных векторов распределена в центре всей области, центр маленький. То есть, при поиске приходится просматривать довольно много лишних векторов.

Простота линейного поиска

Линейный поиск очень прост по сравнению с поиском по индексу — сложность одного шага поиска по индексу сама по себе гораздо больше сложности одного шага «тупого» поиска, так как, скорее всего, включает ветвления, переходы по дереву, запоминание списка найденных элементов и тому подобные операции. Поэтому, даже выигрывая в числе шагов, мы получим меньшее ускорение, чем хотели бы.

Кроме того, с очевидной оптимизацией линейный поиск становится ещё проще. При попарном сравнении всех векторов не нужно их сравнивать полностью! Их достаточно сравнивать с искомой окрестностью последовательно от первой координаты к последней K-ой. Учитывая природу коэффициентов разложения — это примерно соответствует порядку убывания абсолютных значений коэффициентов, а значит, в среднем и разности.

То есть, большая часть неподходящих нам векторов отсекается вообще по 1-ой координате. Относительно мало — по второй, ещё меньше — по третьей, и т. п. В реальности получается, что по коэффициентам дальше 6-го не отсекается почти ничего, а дальше 9-го — совсем ничего. В то же время, для полной проверки подходящего вектора нужно сделать в точности K сравнений и не меньше.

Например, при K=16 все неподходящие вектора (а их 90 %) дают вклад в сложность, примерно равный (а то и меньший) полной проверке оставшихся 10 %. Соответственно, даже если мы изобретём индекс, который позволит нам сразу отсечь почти все 90 % — мы убьём максимум половину сложности, а в реальности — ещё меньше, потому что многомерный индекс никогда не работает точно.

То есть, индекс никогда не отсечёт сразу всё ненужное, а отсечёт, скажем, 85 % и останется 15 %, которые нужно проверить полностью, сложность их проверки будет в 1.5 раза выше половины изначальной. Соответственно, если даже мы сделаем оптимальный неточный индекс, максимум ускорения, получаемого с его помощью, будет зависеть только от процента подходящих векторов и размерности. Чем больше потенциально подходит — тем хуже, чем больше размерность — тем хуже. Например, в случае 85 %, отсечённых мгновенно (сложность = 0), выигрыш будет всего лишь 25 % = (1 — 0.5*1.5).

Из-за того, что по высшим гармоникам фильтрации почти не происходит, целесообразно оказывается использовать 3, максимум 5 коэффициентов разложения. В простейшем случае эти 3 коэффициента — среднее, проекция на синус, проекция на косинус.

Индексные структуры

В простейшей реализации алгоритма поиска повторов все сгенерированные вектора спектральных координат сравниваются друг с другом попарно, что сильно ухудшает производительность. Порядок сложности такого алгоритма — O(N²). Чтобы его ускорить, нужно применять индексные структуры, позволяющие уменьшить сложность поиска.

Для индексации многомерных данных придумано очень много разных структур и алгоритмов. Большинство из них, правда, осталось на теоретическом уровне статей и диссертаций :) в реальности обычно используются самые простые вариации. В целом, все многомерные индексы делятся на два класса:

  • Производные от R-дерева. R-дерево — наиболее известная и используемая в реальности индексная структура (R — от Region), прямое обобщение B-деревьев на многомерный случай. Узлы индекса — «гиперпрямоугольники» многомерного пространства, которые могут содержать либо точки, либо дочерние узлы. R-дерево эффективно работает для минимальных размерностей (2-3) и поэтому широко используется для хранения пространственных данных, например, положений объектов на картах местности. При увеличении размерности, однако, возникает проблема перекрытия узлов и индекс перестаёт давать ускорение, так как при поиске оказывается нужно заглянуть в почти все узлы. Соответственно, большинство производных от R-дерева структур ставят своей целью именно борьбу с перекрытиями узлов.
  • BSP-деревья и их обобщения (BSP = Binary Space Partitioning). Каждый узел представляет собой разбиение пространства на две части, таким образом, части пространства никогда не пересекаются.

Существуют индексы, разрабатываемые специально для использования в СУБД — обычно от таких индексов требуется ориентация на дисковое хранение и обновляемость без значительной потери эффективности. Существуют и более простые структуры, предназначенные для использования в памяти и не имеющие возможности обновления. На самом деле, из реально широко используемых обновляемых многомерных индексов можно упомянуть только R-дерево.

Существует два основных применения индексов для поиска:

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

K-D дерево

Обычное K-D (от K-Dimensional) дерево представляет собой вид BSP-дерева, в котором все разбиения осуществляются гиперплоскостями, перпендикулярными одной из осей координат. Каждый узел K-D дерева характеризуется номером координаты и пороговым значением. Если значение выбранной координаты в точке меньше или равно этого порогового значения, точка относится к левому подпространству, а если больше — к правому.

Выбор разделяющих плоскостей в простейшей и наиболее часто используемой реализации K-D дерева делается просто — выбирается медианное значение нужной координаты по всей выборке индексируемых точек. Таким образом строится сбалансированное K-D дерево. Выбор измерения для разбиения пространства на каждом шаге в простейшем случае производится по кругу, начиная с первого измерения — такой подход даёт наибольшее «равноправие» измерений.

K-D-N дерево

K-D-N дерево:

  • Содержит уровней не больше, чем число измерений.
  • Нелистовые узлы содержат полное разбиение набора данных на диапазоны по одному из измерений, плюс вычисленные на этапе построения дерева минимальное и максимальное значение выбранного измерения у всех попавших в узел элементов. Размеры всех диапазонов одинаковы и равны половине размера ожидаемой окрестности, по которой будет производиться поиск. Больше хуже — массив данных очень плотный, и разделение пространства становится неэффективным. Меньше тоже хуже — сильно возрастает количество узлов дерева, которые нужно просматривать при поиске.
  • Если полный диапазон по какому-то измерению меньше, чем 1.5 * размер ожидаемой при запросе окрестности поиска, разбиение по этому измерению не производится, потому что оно всё равно будет неэффективно — при поиске почти наверняка придётся заглядывать во все дочерние узлы. Почему 1.5? Потому что это 3*ε, если ищем по окрестности (c-ε ≤ x ≤ c+ε). В этом случае пространство разделится на 3 «полосы» каждая размером по ε, запрашиваемый центр окрестности почти наверняка попадёт в центральную из этих полос (помним, что распределение у нас почти нормальное), соответственно, искомые элементы почти наверняка смогут оказаться в любой из трёх полос.
  • Листовые узлы содержат наборы элементов. Технически — списки индексов элементов в изначальном массиве данных. Сами элементы хранить смысла нет, так как нам всё равно нужны индексы, чтобы сопоставлять элементам положения в ДНК-последовательности.
  • Дополнительно создаются и хранятся вариационные ряды (списки индексов элементов), упорядоченные по каждой координате.
  • При поиске просто спускаемся по дереву, просматривая узлы, в которых могут содержаться подходящие элементы, и всё время запоминаем максимальное и минимальное значение по каждому измерению. Так как индексируемый набор данных известен полностью, запоминаются они в виде номера элемента в вариационном ряде — это позволяет избежать лишних операций с плавающей точкой.
  • После поиска нужно произвести окончательную фильтрацию неточного найденного набора элементов. Фильтрация по каждому измерению производится одним из двух способов, оптимальный из которых выбирается сравнением размера фильтруемого набора с количеством элементов, которые придётся отсечь на шаге 2:
    • Либо простым сравнением значения каждого элемента с искомой окрестностью — эффективно, если фильтруемый набор мал.
    • Либо отсечением кусков вариационного ряда по ассоциативному массиву. Для этого двоичным поиском находится положение минимального и максимального значения измерения искомой окрестности (c-ε, c+ε) в вариационном ряде по нужной координате, а потом участки вариационного ряда от «известного» до «найденного» максимального/минимального значения удаляются из результата. Под «известным» диапазоном понимается найденный на этапе обхода дерева. Технически — во временный массив по номерам элементов из этих участков вариационного ряда сохраняется «1», а фильтруемый набор одним проходом фильтруется с помощью проверки элементов временного массива. На самом деле, при многократном поиске временный массив можно не обнулять, а вместо «1» можно использовать увеличиваемый при каждом следующем поиска счётчик. Эта стратегия эффективна, если фильтруемый набор больше, чем сумма размеров удаляемых участков. Обычно это случается, либо если какое-то измерение почти ничего не фильтрует (то есть, если полный диапазон значений по этому измерению сопоставим с размером окрестности поиска), либо если просто найдено очень много элементов.
  • С помощью этих мер достигается эффективность индекса в ситуации, когда пространственный размер массива данных относительно маленький, значимых измерений относительно немного, но при этом нет ограничения на полное число измерений.

Дальнейшие изыскания

Новая идея, но это будет что-то типа гибрида KD- и M-дерева (K-D-M-дерево? :)). Т.е. юзать то же самое кд-дерево, но очень близкие элементы группировать в гиперкубики и индексировать их центры вместо самих элементов.

Можно даже тотально перейти в целые числа.