CC BY
11
ЧЕБЫШЕВСКИИ СБОРНИК
Том 20. Выпуск 2.
УДК 519.72
DOI 10.22405/2226-8383-2019-20-2-93-107
Новый подход к поиску строковой медианы и визуализация
строковых кластеров1
Д. В. Горбачев, Е. П. Офицеров
Горбачев Дмитрий Викторович — доктор физико-математических наук, профессор, кафедра прикладной математики и информатики, Тульский государственный университет (г. Тула).
e-mail: dvgmail@mail.ru
Офицеров Евгений Петрович — кафедра прикладной математики и информатики, Тульский государственный университет (г. Тула). e-mail: eofitserov@gmail.com
где Б С С* — конечный набор строк над алфавитом С,и С С*, редакционное расстояние Левенштейна. Эта задача имеет важные приложения, например в биоинформатике при анализе белковых последовательностей. Однако известно, что в общем случае для и = С* задача о медиане является NP-cлoжнoй. Поэтому для приближенного решения были предложены эвристические алгоритмы, в частности, жадный алгоритм. Мы предлагаем новый гибкий подход, базирующийся па гладкой аппроксимации расстояния Левенштейна й. В его основе лежит стохастическое кодирование символьных последовательностей и следующая формула для редакционного расстояния:
где минимум берется по всем подпоследовательностям (X' ,Х!{)е равной длины. С одной стороны, стохастическое кодирование расширяет класс, на котором ищется экстремум. Однако наш основной результат показывает, что медиана не меняется. С другой стороны, теперь мы можем воспользоваться гладкими методами оптимизации, если заменить минимум в определении выше его гладким приближением. Мы приводим детали реализации па основе градиентного спуска, включая рекуррентные формулы расчета й. Эффективный расчет приближенной медианы позволяет, например, применить метод ^-средних для кластеризации строк. Мы даем способ визуализации этих кластеров на основе метода стохастического вложения соседей tSNE.
Ключевые слова: редакционное расстояние Левенштейна, строковая медиана, гладкий минимум, градиентный спуск, метод к-средних, визуализация строкового кластера.
Библиография: 19 названий.
1 Результаты исследования опубликованы при финансовой поддержке ТулГУ в рамках научного проекта
Аннотация
Рассматривается следующая задача о поиске медианы набора строк:
d(XbX 2)
mm < -
(Х[ ,хц )eCXi ХХ^ 2
{2ii*' - x!2\\i + № Ы*' i + №Ы*2'|},
,Y":1 ИИ* 2018 28.
Для цитирования:
Д. В. Горбачев, Е. П. Офицеров. Новый подход к поиску строковой медианы и визуализация строковых кластеров // Чебышевский сборник, 2019, т. 20, вып. 2, с. 93-107.
CHEBYSHEVSKII SBORNIK Vol. 20. No. 2.
UDC 519.72 DOI 10.22405/2226-8383-2019-20-2-93-107
New approach to searching for string median and visualization
of string clusters2
D. V. Gorbachev, E. P. Ofitserov
Gorbachev Dmitry Viktorovich — Doctor of physical and mathematical sciences, Professor, Department of Applied Mathematics and Computer Science, Tula State University (Tula). e-mail: dvgmail@mail.ru
Ofitserov Evgenii Petrovich — Department of Applied Mathematics and Computer Science, Tula State University (Tula). e-mail: eofitserov@gmail.com
Abstract
We consider the following problem of searching the median of a string set:
c = argmin d(a, b), aeU beD
where D C G* is the finite set of strings over the alphabet G,U c G*, d is the Levenshtein edit distance. This problem has important applications, e.g., in bioinformatics in the analysis of protein sequences. However, it is known that in the general case for U = G* the median problem is NP-hard. Therefore, for an approximate solution, heuristics algorithms were proposed, in particular, the greedy algorithm. We give a new flexible approach based on smooth Levenshtein distance approximation d. It uses stochastic character encoding of sequences and the following formula for the edit distance:
d(xux2) = minv v {2 -^'lli +1*11 - i + - |X2'||,
where the minimum is taken over all subsequences (X' ,X!')e of equal length. On the one hand, stochastic coding extends the class, over which the extremum is searched. However, our main result shows that the median is the same. On the other hand, now we can use smooth optimization methods, replacing the minimum in the definition above with a smooth approximation. We give implementation details based on the gradient descent, including the recursion formulas for calculating d. Effective calculation of the approximate median allows, e.g., to use the k-means method for string clustering. We propose an approach to visualize such clusters based on tSNE stochastic nesting method.
Keywords: Levenshtein edit distance, string median, smooth minimum, gradient descent, k-means, string cluster visualization.
Bibliography: 19 titles. For citation:
D. V. Gorbachev, E. P. Ofitserov, 2019, "New approach to searching for string median and visualization of string clusters" , Chebyshevskii sbornik, vol. 20, no. 2, pp. 93-107.
2The results of the study are published with the financial support of Tulsu within the framework of the scientific project №:HHP_2018_28.
1. Введение
В этой работе мы приведем теоретическое обоснование метода гладкой аппроксимации редакционного расстояния Левенштейна в задаче поиска строковой медианы. Эти результаты продолжают исследования второго автора, связанные с применением методов дифференцируемого выравнивания в задачах машинного обучения [2, 3, 4, 14, 5]. Мы приведем основные факты численной реализации методами математического программирования. Также дается способ визуализации строкового кластера и его приближенной медианы.
Пусть п е N Сп — {д\,... ,дп] — алфавит, |С|
— и число элементов, С^ множество строк или последовательностей а — а[1] ... а[Ь] над Сп длины Ь — |а| > 0 (при Ь — 0 имеем пустую строку А ), С*п — Уь>0
Через й(а1,а2)'- С*п х С*п ^ Z+ обозначим дискретное редакционное расстояние Левенштейна. Оно определяется как минимальное количество операций замены, вставки и удаления символа, необходимых для преобразования строки а\ в строку а2 [1, гл. 11]. Свойства й и разные варианты подсчета приводятся в пункте 2.1.
Строковой медианой или центроидом конечного набора строк И С С^ на подмножестве и С С*п называется строка с — с(И, и) е и, для которой [12]
с — ащшт У^ й(а,Ъ). (1)
аеи ьев
Основными являются случаи и — О и и — С^- Назовем медиану с(П,И) локальной, а с(В,С*п) — глобальной.
Дискретная экстремальная задача (1) имеет большое прикладное значение (см., например, результаты, упомянутые в [9, 12, 10]). Однако, несмотря на то, что сложность вычисления расстояния Левенштейна по алгоритму Вагнера-Фишера полиномиальная (см. пункт 2.1), доказано, что задача поиска глобальной медианы набора И будет МР-сложной, если количество и длины строк в И могут быть произвольными. Для родственной задачи о множественном выравнивании этот факт был установлен в [19]. МР-сложность задачи (1) в работе [9] доказана при помощи сведения к МР-сложной задаче о нахождении длиннейшей общей подпоследовательности строк из И.
Были предложены различные методы для приближенного нахождения медианы. Самый простой способ состоит в поиске локальной медианы (см., например, [12]). Это задача сложности 0(Ь^пах1В\2) где Ьтах — максимальная длина строки из И. Однако результаты [12] показывают, что с(И, И) может заметно отличаться от глобальной медианы с(В,С^.
В [8] был предложен жадный алгоритм приближенного поиска глобальной медианы (см. пункт 2.2). Он прост в реализации, быстр и дает неплохое приближение, которое имеет смысл использовать для сравнения методов. Сложность жадного алгоритма 0(Lmax|D||G|) за итерацию (которых не более 2Ьтах по лемме 1). Были предложены варианты этого алгоритма (см., например, [11])- Отметим вариант получения приближенной глобальной медианы модификацией локальной медианы, имеющий сложность О^т^ОЦС^ за итерацию [12]. В [10] проблема (1) сводится к задаче целочисленного программирования.
В работе [15] применяется вероятностный подход на основе представления символов стохастическими векторами. Нам потребуется такое представление, поэтому напомним, что вектор и> — ('Ш1,...,'ШП) называется стохастическим (вероятностным) или весовым, если и>\,... е [0,1] и wl + ... + — 1. Множество таких векторов обозначим через Шп.
Мы отсылаем к упомянутым выше работам за более подробным обзором проблемы.
В данной работе предлагается метод приближенного вычисления медианы на основе стохастического кодирования последовательностей и гладкой аппроксимацией расстояния Левенштейна (I. Мы предполагаем, что максимальная длина Ьтах подпоследовательности из И много меньше Подробнее идея состоит в следующем.
Пусть Еп = {ei,..., еП}, где ej = (0,..., 0,1, 0 ..., 0) £ {0,1}П — бинарный вектор с единицей на j-м месте, j = 1,п. Закодируем символ gj £ Gn вектором ej. Тогда произвольная последовательность а = а[1] ... a[L] над алфавитом Gn биективно отображается в бинарную последовательность В = В[1] ... В [L] над алфави том Еп.
Непустую последовательность В также можно интерпретировать как бинарную прямоугольную матрицу со строками В [г] = (В [г, 1],..., В[г, п]), г = 1, L, и элементами
„г , | 1, если а[г] = qj, -
В[г,3 ] = ^' И J = 1,п.
0,
В итоге имеем ЕП = СП и можно не делать различия между бинарными последовательностями (матрицами) и символьными последовательностями.
Метрику Левенштейна можно перенести на бинарные матрицы из ЕП с помощью следствия 1 из пункта 2.1:
d(Bi,B2) = min {2 К - B'i\\i + iBil - |Б11 + |Б2| - \Б21}. (2)
(B' ,B'2)eQBlXB2 2 J
Тогда задача о строковой медиане набора D С ЕП на подмножестве U С ЕП эквивалентно записывается как
S(D,U) = min ^ d(A,B), С(D,U) £ arg5(D,U). (3)
AeU bed
Множество ЕП является множеством крайних точек Следуя [15] вложим дискретное множество ЕП в непрерывное пространство WП последовательностей (матриц) вида X = X[1] ... X[L], где X[г] — стохастический вектор:
П
X[г] £ [0,1]П, ^ X[i,j] = 1, i = 1^, L = \Х| ^ 0 j=i
(L = 0 отвечает отстой матрице Л). Будем называть такие матрицы или последовательности весовыми.
Округлением round (w): ЕП ^ стохастического вектора w £ ШП называется бинарный вектор ej, в котором единица стоит на месте максимального Wj (если таких несколько, то выбирается последний). Соответственно, округлением round (X): W* ^ ЕП весовой матрицы X называется бинарная матрица, которая состоит из округлений строк X.
Формула (2) справедлива для произвольных весовых матриц Х\,Х2 £ WЩ- Поэтому задачу (3) можно рассмотреть на более широком множестве U = W*. Следующая основная теорема показывает, что минимум не изменится и решение (3) можно восстановить, применяя округление.
Теорема 1. Для любого конечного набора D с ЕП имеем
S(D, ЕП) = 5(D, W*), С(D, ЕП) = round (С(D, W*)).
Вычисление S(D, W*) сводится к решению серии задач min ^b^d d(X,B), где X £ WП и \Х\ фиксировано. Прострапство W* непрерывно, однако функция X ^ ^BeD d(X,B) не везде дифференцируема из-за взятия минимума в (2). Это не позволяет применить эффективные методы математического программирования, в частности, градиентного спуска. Чтобы обойти этот момент мы воспользуемся гладкой аппроксимацией d расстояния Левенштейна, основанной на применении гладкого минимума. Для эффективного вычисления функции d и
ее градиента предлагаются рекуррентных формулы, являющиеся аналогами формулы динамического программирования для d. Данные результаты были получены вторым автором в работе [14].
Структура работа следующая. В разделе 2 приведены основные утверждения относительно редакционного расстояния Левенштейна: необходимые для понимания классические результаты (пункт 2.1), сведения о жадном алгоритме (пункт 2.2), доказательство основной теоремы 1 (пункт 2.3). В разделе 3 даются основные факты для дифференцируемого редакционного расстояния (ДРР) (1: определение ДРР (пункт 3.1), рекуррентные формулы для расчета (пункт 3.2), детали градиентного спуска (пункт 3.3). В разделах 4 и 5 приведены краткие сведения о визуализации кластеров строк и возможном развитии темы соответственно.
2. Основные утверждения
2.1. Редакционное расстояние Левенштейна
Мы воспользуемся хорошо известными результатами, связанными с расстоянием Левенштейна d, взятыми из [1, гл. 11].
Для строки а — а[1]а[2] .. .а[Ь] е С*п через
а[1] — а[г1,... ,гц\ — а[г\] ... а[ге] С а, 1 ^ г1 < ... <гц ^ Ь,
обозначается подпоследовательность длины I (пустая при I — 0 или Ь — 0), определяемая индексами 1 — ц,... а[г:]} — а[г]а[г + 1]... а[]] — подстрока от индекса г до ] (пустая при г > а[1: г] — префикс длины г — 0, Ь.
Пусть далее а1 — а1 [1].. .а1[Ь1], а2 — а2[1]...а2[Ь2] и
6(i,j ) = d(ai[l: i],a2[l: j]), i = 0,Li, j = 0,L2. Для вычисления расстояния Левенштейна применяется следующая рекуррентная формула:
d(a\, 0,2) = 5(Li,L2), где 5(г, 0) = г, ¿(0, j) = j и при i,j ^ l
5(i,j)=min[5(i,j - 1) + 1,5(i - 1,j) + 1,S(i - 1,j - 1)+ 9}, в Л1, = ^ (4)
[0, ai[t]= a,2 [j ].
Эта формула лежит в основе алгоритма динамического программирования Вагнера-Фишера и позволяет за 0(LiL2) операций вычислить S(Li,L2), а также оптимальное редакционное предписание — список длины d(ai,a,2) необходимых замен, вставок и удалений.
Расстояние Левенштейна удовлетворяет всем аксиомам метрики и
jjaij - MI ^ d(ai,a,2) ^ max(jaij, |«21}. (5)
Можно дать другое выражение для расстояния Левенштейна, используя термин выравнивание последовательностей [1, гл. 11]. Глобальным выравниванием двух строк называется вставка в строки пробелов (в том числе и па их концах) до равной длины так, чтобы г-й символ или пробел первой строки отвечал г-му символу или пробелу второй строки для всех г = 1, 2,... (пара пробел-пробел не допускается). Приведем пример выравнивания из [1] для строк vintner и writers:
/v_intner \wri_t_ers.
О.
По выравниванию легко получить редакционное предписание и наоборот: совпадающая пара символ-символ пропускается, несовпадающая пара символ-символ отвечает замене, пара пробел-символ — вставке и символ-пробел — удалению. Для примера выше имеем предписание длины 5, которое является оптимальным.
Произвольное выравнивание последовательностей ai, а2 можно построить, задав список пар символ-символ. Символы этих пар образуют подпоследовательности одинаковой длины. Для примеры выше имеем viter и witer. Оставшиеся символы будут отвечать пробелам. В связи с этим введем следующее определение (где «е» в индексе означает equal length).
Определение 1. (1) Соответствием (^,а'2)е ^ ai х а2 последовательностей а\, а2 называется произвольная пара подпоследовательностей ^ = ai..., i'f], а'2 = a2[i'{,..., г''] равной длины, I = | = |a2'I, где 0 ^ I ^ min{Li,L2} (I = 0 отвечает пуст,ом,у соответствию). Здесь символ ai [¿'J первой последовательности отвечает символу a2[i"s] второй последовательности, s = 1,1.
(2) Также положим а^ = а^[1]... а^[I], а'2 = а'2[1]... а'2[И\, где а^[s] = ai[г'<,]; а'2[s] = a2[i"].
(3) Соответствие (а\,а'2)е однозначно определяется индексам,и подпоследовательностей (i', i")e S 1,^1 х 1,L2, где i' = ii,...,i'V/ i" = i'[,...,i'' (для пуст,ого соответствия имеем пустой набор индексов).
В выравнивании, построенном по соответствию (а\, а'2)е) — I символов первой строки, не вошедших в соответствие, должны отвечать пробелам второй строки, а L2 — I символов второй строки — пробелам первой. Например, для соответствия (vine,ritе)е имеем выравнивание
/_vintner__\
.
Несмотря на то, что пробелы расставляются неоднозначно, длина предписания будет всегда равна 8. Она складывается из 2 различий в соответствии (замены), 3 пробелов вверху (вставки) и 3 пробелов внизу (удаления).
Таким образом, в общем случае длина предписания, отвечающего соответствию
(a'i,a'2)е ^ «1 х а2,
будет содержать |ai| — | позиций для вставок и |a2| — |а21 Для удалений. При этом число замен определяется числом различий в подпоследовательностях а^, а^, которое можно задать с помощью расстояния Хемминга
dH(а'1,а!2) = |{г: а[[г] = а^[г], i = М}|. (равно нулю при I = 0). В итоге, получаем следующее
Предложение 1. Расстояние Левенштейна между cm,рокам,и а\,а2 s G*n равно d(a\,a2)= min {dn(a[,a'2) + |ai| — | + |a2| — |a21} .
(a[,a^)eCai x«2
Отметим, что мы не встречали это утверждение в литературе.
Предложение 1 следующим образом обобщается на бинарное представление последовательностей.
Следствие 1. Пусть Bi,B2 s Е*. Тогда
d(Bi, В2) = min н {2 № — В2+ ^| — ^| + ^ — №||, (6)
где для, I = \В[\ = \В2,'\
I п
1В -ВЧЬ = £ \\В[\г] -ВМЦ = ££ \В[[г,— В'Цг,
И
г=1 г=1]=1
(пуль при I = 0).
Доказательство. Достаточно показать, что ||В1 — В''||1 = 2йн(В^В'^). Это легко следует из следующего равенства для строк матриц:
11 [] 2 [Л| 1 2, В[[г] = В"\г], ,
Оно верно, так как в каждом из бинарных векторов В1\^, В''[г] только по одной едипице. □
Следствие 2. (1) Формула (6) позволяет определить редакционное расстояние (1(Х1,Х2) для, произвольных весовых матриц Х1,Х2 € Ш*.
(2) Функция (1(Х1,Х2) по-прежнему удовлетворяет неравенствам (5).
Доказательство. Часть 1 следствия очевидна.
Часть 2 вытекает из (6). Пусть для определенности \Х1 ^ \Х2\. Соответствием максималь-( Х1 , Х ) | Х1 | = | Х |
ё(Х1,Х) < 2 ||Х1 — Х''|1 + \Х1\ — \Х1\ + \Х2\ — \Х''\ < 1 (|Х1|1 + ||Х''||1) + \Х2\ — \Х''\.
Имеем
|Хх| п |Хх|
||Х1|1 = ^^ \Х1 [г,;]\ = £ 1 = \Х1\,
г=1 ]=1 г=1
| Х | 1 = \ Х \
й(ХиХ2) < \Х''\ + \Х2\ —\Х''\ = \Х2\.
С другой стороны, для любых соответствий
2 ||Х1 — Х''||1 + \Х1\ — \Х1\ + \Х2\ — \Х''\ > \Х1\ + \Х2\ — 2\Х1 \ = \Х2\ — \Х1\,
откуда следует, что с1(Х1,Х2) ^ \Х2\ — \Х1\ и, аналогично, с1(Х1,Х2) ^ \Х1\ — \Х2\. □
Для доказательства теоремы 1 нам потребуется следующая лемма (где «■» обозначает скалярное произведение векторов).
Предложение 2. Пусть Х € Ш*, В € Е*п. Тогда
1^1
й(Х, В) = \Х\ + \В\ — шах _[г] Х[г], (7)
(Ц")еС1,|Х|х1,|Б| ^ ' где 1' = %!,..., 1' = г'1,..., г'' — индексы, соответствий, 0 ^ I ^ шт{ \Х\ , \ В\};
Доказательство. Используя индексы (1', \")е для задания соответствия
(Х',В")е, X' = X [1'], В" = В[М'},
по формуле (6) находим
<1(Х,В) = шт {2 УХ' - В''\\! + \Х| - \Х'| + 1В| - |Б''|}.
(1',1")е I. 2 )
Учтем, что X — весовая матрица, а В — бинарная. Поэтому X[1'3,]\ € [0,1],
Ех к ,J'] = 1
3=1
и B[i"s] = ejs для некоторого js. Отсюда следует, что 1 1 1 п
11
"X' — — \Х'| — |Б''| = 2 Е Е 1ХК,j] — В[%",j] I — 21
s=1 j=1
I / ч I
w \ 1 1 1
2 E E X ] + 1 — X Ы) — 21 = 1 E(2 — 2X [CjJ) — 21 = — E (1 + X [CjJ)
j=js ' i=1 i=1
I / n n \ In
EE x ] + E в[г" ,j ]x v.,3 ] =—EE x к,i ] + 1) =
S=1 ^ j=1 j=1 ' S=1 j=1
j=1 j=1 S=1 j=\
\X \ n \X |
EE x [', №>¡1' [*,з] =—E Bi'y w ■x w ■■
i=1 j=1 i=1
Таким образом,
M
\\X' — B''\U + | X | — |X' | + |R| — |R''| = | X | + Ш | —
1 "X' — B''\\1 + |X| — |X'| + \B\ — \B''| = |X| + \B| — EBVtV.[г] ■ X[г],
i=1
откуда следует утверждение леммы. □
Как хорошо известно максимум конечного числа линейных функций — выпуклая кусочно-линейная функция (см. [7, пример 3.5]). Отсюда и из предложения 2 немедленно вытекает
Следствие 3. Для каждого г = 1,Ь функция X[г] ^ й(Х[1] ...X[г] ...X[Ц,В) на, множестве Шп является выпуклой (вниз) кусочно-линейной функцией.
2.2. Жадный алгоритм
Жадный алгоритм поиска строковой медианы был предложен в работе [8]. Напомним, что решается дискретная задача оптимизации
с = argminF(a), F(а) = Е d(a,b), D С G*
aeG™ beD
Для известного а функция F(а) вычисляется эффективно по алгоритму Вагнера-Фишера.
Lemma 1. Пусть Lо = щу beD Н ~ средняя длина cm,роки в наборе D. Тогда, для, длины медианы справедлива оценка,
|с| < 2LD. (8)
Доказательство. Для пустой строки имеем F(А) = b<^D НЕсли |с| > |dd| beD т0 в СИЛУ (5) находим
f (с) = Y.d(с,b) > Е iici - iftii > |Е(и - ifti)| = Им - Е N > Е = F (А).
beD beD beD beD beD
Поэтому с не может быть медианой и (8) верно. □
Шаги жадного алгоритма. Начинаем с присваиваний с0 = А и Fo = F(с0) = beD М-Далее для к = 0,1,... решается задача
Cfc+i = argmin {F(а): а = скд, д е G} , Fk+1 = F(ck+i),
т.е. перебором по д е Gn к текущей строке Ск добавляется символ д и среди таких строк а = Скд вычисляются минимальное значение Fk+i функции F(а) и соответствующая оптимальная строка Ск+1- Сложность этого шага 0((к + 1)Lmax|G||D|), где Lmax = maxbeD Щ-Если Fk+i < min^=okFi, то итерации продолжаются и Ck+i новое приближение медианы.
к
Число итераций в силу леммы 1 ограничено 2Ld-
2.3. Доказательство теоремы 1
Пусть D,U С Е*, Fd(X) = YlBeDd(X,B). Напомним, что в матричной постановке медиана С = С (D, U) определяется из решения дискретной задачи
S(D,U) = minFD(X), С = С(D,U) е argS(D,U). х eu
Покажем, что
S (D,E* )=S(D,W.П). По лемме 1 с учетом следствия 2 находим
minFd(X) = min min Fd(X), (9)
хeu o^l^2ld xeu
" " d |x|=L
поэтому достаточно доказать
ДL) = S(D, W*L
S (D,EL) = S (D,WL). (10)
При L = 0 это равенство очевидно и ревизуется на пустой матрице. Поэтому пусть L ^ 1.
Очевидно, что минимум в непрерывной задаче S(D, достигается па некоторой весовой матрице X £ W,^. Этот минимум можно искать следующим образом:
S(D, Wt) = min .. min FD(X) X [1]£ff„ x [L]ewn
В силу следствия 3 функция X[г] ^ Fd(X[1] .. .X[г] .. .X[L]) выпуклая па множестве Wn. Поэтому ее минимум достигается в некоторой крайней точке (бинарном векторе) XX[г] = ej0 £ En. Это доказывает (10).
В общем случае минимум может также достигаться в выпуклой комбинации крайних точек X[i] = wjsejs> 1 ^ т ^ п. Тогда по определению округления вектор round (XT[f]) совпа-
дает с одним из ejs и также дает решение. Следовательно, round (XX) — решение S(D,EL).
Теорема 1 доказана.
3. Гладкая аппроксимация расстояния Левенштейна
Детали гладкой аппроксимации расстояния Левенштейна подробно даны в работе [14]. Здесь мы приведем основные факты.
3.1. Гладкое расстояние
По теореме 1 медиану С £ Ei** дискретного набора D С Ei** можно искать в непрерывном пространстве весовых матриц W*. Для этого с учетом (9) для каждого 0 ^ L ^ 2Ld ищется медиана Cl заданной длины L:
CL = round (Cl), Cl = argmin FD(X), (11)
x ew-ь
где, напомним, Fd(X) = ^^^d d(X, В). Глобальная медиана С определяется как лучшая из медиан Cl-
Для задачи оптимизации на непрерывном пространстве W*L естественно воспользоваться гладкими методами оптимизации. Однако функция X ^ d(X, В) в силу следствия (3) является только кусочно-дифференцируемой. Поэтому в [14] было предложено заменить в определении (6) минимум на его гладкую аппроксимацию, называемую гладким минимумом. Этот прием известен в выпуклой оптимизации [7].
Для конечного множества Y С К гладкий т-минимум определяется равенством
smin (Y; г) = ^ ^, (12)
где т < 0 — параметр аппроксимации (при т > 0 имеем гладкий максимум). Нетрудно показать, что
min Y ^ smin (Y; т) ^ min Y + СетА, (13)
где С = \Y|(mean Y — min Y) и А = miny/=y«eY W — y"V Отсюда следует, что гладкий минимум стремится к обычному экспоненциально быстро. Далее параметр т часто опускается, если это не вызывает недоразумений.
Определение 2. Пусть Xi,X2 £ W,т < 0. Гладкой аппроксимацией расстояния Левенштейна называется функция
d(X\,X2; г) = smin {2 К — Х2||i + № — |ХЦ + |Х2| — |}. (14) (x'1,x'2)eCXixX^ 2 J
Из определения гладкого минимума следует, что функция d(Xi,X2) будет дифференцируемой по аргументам, понимая под градиентом функции ||Х || i = ^ L=iYl ,=i Is, АI матричную функцию
V||X ||i = (sign X [s,t]), s = 1Д, t = T/n. (15)
Поэтому d также называется дифференцируемым редакционным расстоянием (ДРР).
ДРР обладает свойствами симметричности и неотрицательности. Однако, в отличии от оригинального редакционного расстояния, d(X, X) = 0. Поэтому также вводится следующее
Определение 3. Несмещенным ДРР называется функция
d0(Xi,X2) = d(Xi,X2) — 2 {d(Xi, Xi) + d(X2,X2)).
Из (13) следует
Предложение 3. При т ^ —те экспоненциально быстро d0(Xi,X2; т) ^ d(Xi,X2).
3.2. Рекуррентные формулы
Расчет ДРР ! напрямую по формуле (14) затруднителен из-за перебора соответствий, число которых растет экспоненциально с ростом и |Х2|. Напомним, что для классического редакционного расстояния ! данная проблема решается при помощи рекуррентной формулы (4) и основанной на ней алгоритме динамического программирования Вагнера-Фишера, который позволяет вычислить ! за полиномиальное время.
В [14] аналогичные рекуррентные формулы установлены для ! Пусть Х\,Х2 € 1Х\1 = Ь\, |Х21 = г = 0,Ь\, ] = 0,Ь2, — множество всевозможных соответствий (Х1, Х2')е € Х\[1 : г] х Х2[1 : ]} и для таких соответствий длины к = |ХЦ = |Х2'|
Пг]{х'1,х'2) = 1 ||Х1 - Х2' ||1 + 1 + 3- 2к.
Положим
ац = агз(Хъ Х2) = Е К^(Х1, Х2')^(ВД'),
) е €Пгj
ргз =^(Х1,Х2)= Е ^(Х1,Х2'),
чт0 отвечает числителю и знаменателю в гладком минимуме (12). Тогда
аЬ\Ь2
й(Х1,Х2) =
Теорема 2 ([14]). Для г = 1,ЬЪ у = 1,Ь2
а^ = (аг-1,у +а^-1 + ^-1)ет + (аг-1^-1 + Рг-1^-16^)ет6^ - (аг-1^-1 + 2^-^-1)е2т,
ргз = (¡Зг-1,3 + 13г,3-1)ет + рг-1,3-1(ет&^ - е2т),
где
ам = гетг, а0з = зет*, ро = етг, ^ = етз, 5гз = 2 ЦХ1М -Х2Ш1.
Рекуррентные формулы из теоремы 2 дают эффективный способ вычисления !(Х1,Х2) за полиномиальное время. Дифференцирование этих формул с учетом (15) позволяет получить рекуррентные соотношения для эффективного расчета градиента Vх1!(Х1, Х2). Отсюда получаем значение и градиент несмещенного ДРР.
3.3. Градиентный спуск
Рассмотрим задачу приближенного вычисления медианы Сь на множестве (см. (11)). Из определения имеем задачу условной оптимизации. Возможность вычисления градиента несмещенного ДРР позволяет воспользоваться методом проекции градиента.
В качестве начального приближения Хо берется случайная матрица из . Очередная итерация имеет вид
Хк+1/2 = Хк - Ъ -1- £ Vx(!о(Х,В))(Хк), Хк+1 = Рт(Хк+1/2), (16)
1 1 вео
где 1к > 0 — шаг метода и проекция Рг: МЬхга ^ определяется равенством
Рг (У) = aгgmin ||У -Х||2.
хеш^
Для практических приложений важен случай, когда число строк |_D| велико. В этом случае можно применить стохастический градиентный спуск, когда в (16) вместо D берется случайная подвыборка из D много меньшего размера [16].
Было проведено сравнение предложенного подхода с жадным алгоритмом. Для этого проведен вычислительный эксперимент на синтетических данных. Генерировалось 100 наборов случайных строк D над алфавитом из 4 букв длинной от 5 до 20 символов. Далее па этом множестве выполнялся поиск медианы двумя алгоритмами. Оценивалась величина
5(D) = Е (d(С, В) — d(Cg, В)),
BeD
где С — медиана, найденная предложенным алгоритмом, aCs — жадным. Было установлено, что для наборов из 10 строк mean 5 = —0.23, std 5 = 0.034 па уровне значимости р = 0.99. Это означает, что в данном примере точность алгоритма чуть выше жадного. Однако с увеличением числа строк они сравниваются по точности, поэтому требуется дальнейшее исследование проблемы.
4. Визуализация кластеров строк
Одно из интересных и важных приложений строковой медианы является визуализация строковых кластеров. Эта задача имеет прикладное значение, например, в биоинформатике. При этом часто достаточно знать медиану приближенно, например, при визуализации больших данных.
Рассмотрим данную задачу. Имеется набор строк S. Требуется разбить S па непересекающиеся кластеры S = (J^=1 Sk, так чтобы строки каждого кластера были близки по редакционному расстоянию к соответствующим медианам c(S^)• В случае векторных данных для решения этой задачи применяется популярный алгоритм k-средних (см., например, [6, гл. 3]), позволяющий находить заданное число кластеров. Он основан на нахождении центроидов. Для строк центроидом является медиана. Поэтому для адаптации метода k-средних к строкам важно эффективно находить медиану, пусть приближенно. Результаты данной работы позволяют это сделать. Отметим еще алгоритм k-медоидов, в котором среднее ищется среди точек кластера. В случае строк это будет отвечать локальной медиане. Однако этот метод имеет свои ограничения, в частности, он требует знания матрицы расстояний ( d(Bi,Bj))'^\=1-Кроме того, как уже отмечалось, локальная медиана может быть далека от глобальной, что может привести к смещению кластеров.
Чтобы проанализировать получившиеся кластеры, интересно их визуализировать. Основная проблема состоит в большой размерности данных. Стандартным подходом в этом случае является проекция данных в пространство низкой размерности (два или три) с сохранением метрических свойств (см. классические работы [18, 17]).
Популярным методом визуализации, активно используемом в машинном обучении, является метод стохастического вложения соседей с t-распределением (tSNE), предложенный в работе [13]. Данный метод позволяет визуализировать многомерные данные путем их нелинейной проекции на пространство малой размерности. В отличии от классических линейных алгоритмов снижения размерности [18] данный метод лучше отображает структуру многомерных данных. Он позволяет рассчитать координаты проекций строк на основе матрицы парных расстояний. В данном случае вычисление такой матрицы не является проблемой, поскольку для визуализации достаточно использовать небольшие подвыборки кластеров.
На рисунке приведен пример визуализации кластеров строк. Он сделан в Python с помощью класса TSNE из пакета машинного обучения Sklearn. Двумерный график (слева) получен функцией scatter из пакета Matplotlib, предназначенного для построения 2В-графиков. Для
вывода трехмержнх) графика (справа) применялась функция points3d из пакета научной 3D-визуализации Mayavi.
5. Заключение
В работе предложено теоретическое обоснование новохх) подхода для поиска строковой медианы с помощью гладких методов оптимизации. Основным преимуществом данного подхода по сравнению с существующими дискретными методами является мх) гибкость. Это позволяет применять ех'о в задачах поиска медианы с ограничениями, а также во многих других связанных проблемах строковой оптимизации. Благодаря тому, что ДРР является диффрен-цируемой функцией ее можно использовать для сравнения последовательностей в составе различных дифференцируемых моделях машинших) обучения, включая нейронные сети. Такие задачи строковой оптимизации в том числе играют важную роль при визуализации строковых данных. Планируется исследовать перечисленные направления.
1. Гасфилд Д. Строки, деревья и последовательности в алх'оритмах: Информатика и вычислительная биология. СПб.: Невский Диалект, БХВ-Петербург, 2003.
2. Офицеров Е.П. Статистическая модель периферической селекции Т-клеточных рецепторов /7 Известия ТулГУ Технические науки. 2017. Вып. 2. С. 138 143.
3. Офицеров Е.П. Глубокая модель селекции Т-клеточных рецепторов /7 Известия ТулГУ. Технические науки. 2017. Вып. 12. Ч. 2. С. 350 355.
4. Офицеров Е.П. Классификация последовательностей на основе коротких мотивов. Чсбы-шевский сборник. 2018. Том 19, №:1. С. 187 199. doi: 10.22405/2226-8383-2018-19-1-187-199
5. Офицеров Е.П. Программный комплекс для решения задач машинших) обучения на строковых данных при помощи дифференцируемо!^ редакционшнх) расстояния /7 Известия ТулГУ Технические науки. 2019. Вып. 5. С. 370 376.
6. Ту Дж., Гонсалес Р. Принципы распознавания образов. М.: Мир, 1978.
v.V.
-10 -
СПИСОК ЦИТИРОВАННОМ ЛИТЕРАТУРЫ
7. Boyd S., Vandenberghe L. Convex Optimization. Cambridge University Press, 2004.
8. Casacuberta F., de Antonio M. A greedy algorithm for computing approximate median strings // In: VII Simposium Nacional de Reconocimiento de Formas y Análisis de Imágenes. 1997. P. 193-198.
9. De la Higuera C., Casacuberta F. Topology of strings: Median string is NP-complete // Theoretical Computer Science. 2000. Vol.* 230, no. 1-2. P. 39-48. doi: 10.1016/s0304-3975(97)00240-5
10. Havashida M., Kovano H. Finding median and center strings for a probability distribution on a set of strings under Levenshtein distance based on integer linear programming //In: Fred A., Gamboa H. (eds) Biomedical Engineering Systems and Technologies. BIOSTEC 2016. Communications in Computer and Information Science. Springer, Cham. 2017. Vol 690. P. 108— 121.
11. Kruzslicz F. Improved greedy algorithm for computing approximate median strings // Acta Cvbernetica. 1999. Vol. 14, no. 2. P. 331-339.
12. Martínez-Hinarejos C.D., Juan A., Casacuberta F. Use of median string for classification // Proc. 15th Int. Conf. on Pattern Recognition. ICPR-2000, Barcelona, Spain. 2000. Vol. 2. P. 903-906. doi: 10.1109/ICPR.2000.906220
13. Van der Maaten L., Hinton G. Visualizing data using t-SNE //J. Mach. Learn. Res. 2008. Vol. 9. P. 2579-2605.
14. Ofitserov E., Tsvetkov V., Nazarov V. Soft edit distance for differentiable comparison of symbolic sequences // arXiv:1904.12562. 2019.
15. Olivares-Rodríguez C., Oncina J. A stochastic approach to median string computation // In: da Vitoria Lobo N. et al. (eds) Structural, Syntactic, and Statistical Pattern Recognition. SSPR /SPR 2008. Lecture Notes in Computer Science. Springer, Berlin, Heidelberg. 2008. Vol. 5342. P. 431-440.
16. Ruder S. An overview of gradient descent optimization algorithms // arXiv:1609.04747. 2016.
17. Sammon J.W. A nonlinear mapping for data structure analysis // IEEE Transactions on Computers, 18(5):401-409, 1969.
18. Torgerson W.S. Multidimensional scaling I: Theory and method // Psvchometrika. 1952. Vol. 17. P. 401-419.
19. Wang L., Jiang T. On the complexity of multiple sequence alignment //J. Computat. Biol. 1994. Vol. 1, no. 4. P. 337-348.
REFERENCES
1. Gusfield D. Algorithms on strings, trees, and sequences: Computer science and computational biology. New York: Cambridge University Press, 1997.
2. Ofitserov E.P. Statistical model of T-cell receptor pereferic selection // Izvestiva of Tula State University. Technical sciences. 2017. No. 2. P. 138-143.
3. Ofitserov E.P. Deep model of T-cell receptor selection // Izvestiva of Tula State University. Technical sciences. 2017. No. 12, part 2. P. 350-355.
4. Ofitserov E.P. Motif based sequence classification // Chebvshevskii Sbornik. 2018. Vol. 19, no. 1. P. 187-199. (In Russ.) https://doi.org/10.22405/2226-8383-2018-19-l-187-199
5. Ofitserov E.P. Software package for solving machine learning problems on string data using soft edit distance // Izvestiva of Tula State University. Technical sciences. 2019. No. 5. P. 370-376.
6. Tou J.T., Gonzalez R.C. Pattern recognition principles. Addison-Wesley Publishing Company, 1974.
7. Bovd S., Vandenberghe L. Convex Optimization. Cambridge University Press, 2004.
8. Casacuberta F., de Antonio M. A greedy algorithm for computing approximate median strings // In: VII Simposium Nacional de Reconocimiento de Formas y Análisis de Imágenes. 1997. P. 193-198.
9. De la Higuera C., Casacuberta F. Topology of strings: Median string is NP-complete // Theoretical Computer Science. 2000. Vol. 230, no. 1-2. P. 39-48. doi: 10.1016/s0304-3975(97)00240-5
10. Havashida M., Kovano H. Finding median and center strings for a probability distribution on a set of strings under Levenshtein distance based on integer linear programming //In: Fred A., Gamboa H. (eds) Biomedical Engineering Systems and Technologies. BIOSTEC 2016. Communications in Computer and Information Science. Springer, Cham. 2017. Vol 690. P. 108— 121.
11. Kruzslicz F. Improved greedy algorithm for computing approximate median strings // Acta Cvbernetica. 1999. Vol. 14, no. 2. P. 331-339.
12. Martínez-Hinarejos C.D., Juan A., Casacuberta F. Use of median string for classification // Proc. 15th Int. Conf. on Pattern Recognition. ICPR-2000, Barcelona, Spain. 2000. Vol. 2. P. 903-906. doi: 10.1109/ICPR.2000.906220
13. Van der Maaten L., Hinton G. Visualizing data using t-SNE //J. Mach. Learn. Res. 2008. Vol. 9. P. 2579-2605.
14. Ofitserov E., Tsvetkov V., Nazarov V. Soft edit distance for differentiable comparison of symbolic sequences // arXiv:1904.12562. 2019.
15. Olivares-Rodríguez C., Oncina J. A stochastic approach to median string computation // In: da Vitoria Lobo N. et al. (eds) Structural, Syntactic, and Statistical Pattern Recognition. SSPR /SPR 2008. Lecture Notes in Computer Science. Springer, Berlin, Heidelberg. 2008. Vol. 5342. P. 431-440.
16. Ruder S. An overview of gradient descent optimization algorithms // arXiv:1609.04747. 2016.
17. Sammon J.WT. A nonlinear mapping for data structure analysis // IEEE Transactions on Computers, 18(5):401-409, 1969.
18. Torgerson WT.S. Multidimensional scaling I: Theory and method // Psvchometrika. 1952. Vol. 17. P. 401-419.
19. Wang L., Jiang T. On the complexity of multiple sequence alignment //J. Computat. Biol. 1994. Vol. 1, no. 4. P. 337-348.
Получено 18.05.2019 г. Принято в печать 12.07.2019 г.
