Алгоритм Витерби

Приёмник GPS в машине ошибается на 5–15 метров, а соседние улицы, дублёры и съезды бывают в нескольких десятках метров друг от друга. Если каждую отметку отнести к ближайшей улице, машина на карте начнёт перескакивать на соседние улицы и сворачивать на перекрёстках туда, куда не ездила. Можно иначе: смотреть на весь трек сразу и искать маршрут, по которому машина могла проехать и который лучше всего объясняет все отметки вместе. Если перебирать улицы-кандидаты для каждой из 170 отметок поездки ниже, комбинаций — около 10⁸⁹. Алгоритм Витерби находит лучшую за миллисекунды.

Машина едет по кратчайшему пути от вокзала Алматы-2 до MEGA Alma-Ata, 10,9 км. Приёмник записывает положение раз в 5 секунд и ошибается в среднем на 15 м по каждой оси. Ближайшая улица верна для 86,5 % отметок, Витерби — для 94,7 %. Карта открывается там, где ближайшая улица ошибается чаще всего; красные кольца с пунктиром — места, куда она отнесла бы отметки. Больше половины её ошибок на этих маршрутах — параллельная дорога в нескольких метрах: встречная половина широкой улицы с разделительной полосой или дублёр. В данных OpenStreetMap это отдельные дороги, и отнести к ним машину — тоже ошибка: по встречной половине она ехала бы в другую сторону. Кнопка «к следующей ошибке» обходит ошибки выбранного метода, а касание отметки показывает улицы-кандидаты и ту, которую выбрал Витерби.

Идея

Машина в каждый момент в каком-то состоянии — на какой-то улице, — но состояния мы не видим. Видим наблюдения: отметки GPS. Модель говорит две вещи: насколько правдоподобна отметка, если машина на такой-то улице, и насколько правдоподобен переход с одной улицы на другую между двумя отметками. Ищем цепочку состояний, при которой всё вместе правдоподобнее всего. Для этого на каждом шаге достаточно помнить по одной цепочке на состояние — лучшую из тех, что в нём кончаются. Остальные уже не понадобятся.

Такую модель называют скрытой марковской моделью. В ней есть вероятности переходов a(r, s) — из состояния r в состояние s за один шаг, вероятности наблюдений b(s, o) — увидеть o, находясь в s, и вероятности начать с каждого состояния π(s). Вероятность того, что состояния были s₁, …, s_T и при этом наблюдались o₁, …, o_T, — произведение:

Наблюдения уже даны, поэтому самая вероятная при них цепочка состояний — та, у которой это произведение наибольшее. Произведения тысяч вероятностей быстро уходят за пределы чисел с плавающей точкой: 0,25^{48 502} — примерно 10^{−29 201}, а обычные числа double кончаются около 10⁻³⁰⁸. Поэтому берут логарифмы, и произведение превращается в сумму очков: очки наблюдения e_t(s) = log b(s, o_t) и очки перехода log a(r, s). Наибольшая сумма — у той же цепочки, что и наибольшее произведение.

Обозначим V_t(s) наибольший счёт цепочки из первых t состояний, которая кончается в s. Последний шаг такой цепочки пришёл из какого-то состояния r, поэтому

Таблица заполняется шаг за шагом, и в каждой клетке запоминают, из какого r пришёл максимум. Наибольшее из V_T — счёт лучшей цепочки, а саму цепочку восстанавливают с конца по запомненным стрелкам. Картину из столбцов-шагов, в каждом по клетке на состояние, со стрелками-переходами между соседними столбцами называют решёткой. В поездке выше у отметки в среднем 4,2 улицы-кандидата, у некоторых — восемь; решётка — 170 столбцов, в среднем по четыре клетки, а путей по ней — те самые 10⁸⁹.

Почему цепочка лучшая

Всё держится на одном свойстве: очки перехода и наблюдения на шаге t зависят только от состояний на шагах t − 1 и t, а не от того, как цепочка в них пришла. Это и есть марковость. Если бы штраф зависел от трёх последних улиц, состоянием пришлось бы считать пару улиц — состояний стало бы больше, но алгоритм остался бы тем же.

На языке графов это кратчайший путь по решётке: вершины — клетки, рёбра — переходы, длина ребра — очки перехода и наблюдения со знаком минус. Решётка — граф без циклов, и столбцы идут в топологическом порядке, поэтому куча Дейкстры не нужна: каждую клетку достаточно посмотреть один раз. Так же устроен шов на картинке, только там у пикселя три соседа сверху, а здесь в состояние можно прийти из любого.

Улицы по отметкам GPS

Скрытую марковскую модель для привязки трека к дорогам в таком виде описали Пол Ньюсон и Джон Крумм в 2009 году. Состояния на шаге t — улицы рядом с отметкой z_t, точнее ближайшие к отметке точки этих улиц. Очки наблюдения — нормальный закон для расстояния d от отметки до точки улицы:

где σ — ошибка приёмника; множитель перед экспонентой одинаков для всех кандидатов и на выбор не влияет. Очки перехода сравнивают два расстояния: по прямой между отметками z_t−1 и z_t и по дорогам между точками-кандидатами. Если машина и правда ехала по этим улицам, расстояния почти равны: за несколько секунд она проезжает примерно столько, на сколько сдвинулась. У неверной пары кандидатов путь по дорогам выходит длиннее — чтобы перебраться на параллельную улицу, нужно доехать до перекрёстка и вернуться. По поездкам с известным маршрутом Ньюсон и Крумм увидели, что разница расстояний распределена по показательному закону, отсюда

Ещё две детали из их статьи. Улицы дальше 200 метров от отметки в кандидаты не берут. А отметку, которая ближе 2σ к предыдущей оставленной, выбрасывают: пока машина не сдвинулась дальше, чем может ошибиться приёмник, отметка говорит больше о шуме, чем о движении. В их данных так ушло 38,9 % отметок. Параметры оценили по поездке, привязанной вручную, — около 80 километров за два часа, 7531 отметка раз в секунду. Ошибка приёмника вышла σ = 4,07 м; её считали через медиану отклонений, чтобы редкие большие выбросы не испортили оценку.

На этой поездке алгоритм дал ровно тот же маршрут, что и ручная привязка. Потом данные портили. Когда оставляли одну отметку из тридцати, доля неверно привязанного пути выросла всего до 0,11 %. А к шуму алгоритм оказался чувствительнее при записи раз в секунду, чем при редкой: частые шумные точки, как пишут авторы, тянут маршрут то в одну, то в другую сторону. При редких отметках он выдерживал шум до 50 метров.

В сцене то же самое на графе дорог центра Алматы. Трек смоделирован: машина едет по кратчайшему пути с разрешённой на улице скоростью, а к её положению добавлен случайный шум — поэтому для каждой отметки известна настоящая улица. Кандидаты — до восьми улиц в радиусе 4σ, но не меньше 50 м, β = 15 м. Путь по дорогам между кандидатами ищет Дейкстра из обоих концов улицы-кандидата и обрывает поиск, когда путь длиннее трёх расстояний по прямой плюс 500 м: такой переход всё равно почти невероятен. На маршруте от вокзала это 170 шагов, 703 поиска и 3382 пары кандидатов.

Переключайте ошибку и частоту записи. При ошибке 5 м ближайшая улица и так почти всегда верна, а Витерби доводит долю до 98–100 %. При 30 м и записи раз в 15 секунд ближайшая улица верна для 71,7 % отметок, Витерби — для 83,3 %. А при 30 м и записи раз в секунду он даже чуть хуже ближайшей улицы: 56,5 % против 58,2 %. Это тот же эффект, что у Ньюсона и Крумма. После прореживания на 2σ = 60 м между оставленными отметками машина проезжает ненамного больше, чем ошибается приёмник, и сравнение расстояний почти ничего не говорит. К тому же первой дальше 60 м от предыдущей часто оказывается отметка, которую шум отнёс особенно далеко.

Привязку по этой модели делают и в открытых маршрутизаторах: сервис match в OSRM ссылается на статью Ньюсона и Крумма и считает по ней путь Витерби, так же устроен модуль Meili в Valhalla.

Участки генома

Тот же алгоритм размечает ДНК. В геноме фага λ 48 502 буквы, и доля G и C меняется по его длине: где-то больше 55 %, где-то около 40 %. Модель с двумя скрытыми состояниями описывает это так. В состоянии «богато G и C» каждая из букв G и C встречается с вероятностью 29 %, A и T — по 21 %; в состоянии «богато A и T» — 20 % и 30 %. Очки буквы — логарифм отношения её вероятности в состоянии к 1/4, в сотых долях бита: G в первом состоянии приносит 100 · log₂(0,29 / 0,25) ≈ 21, A — 100 · log₂(0,21 / 0,25) ≈ −25. Смена состояния стоит штраф — логарифм вероятности перехода с обратным знаком.

При цене смены 2000 — это 20 битов, то есть вероятность перехода около одной миллионной на букву, — путь делит геном на шесть участков. Левые 21 923 буквы богаты G и C, 56,8 %, и на них лежат все гены головки и хвоста. Дальше — участок с 40,1 % G и C, где гены встраивания в хромосому int и xis и гены с неизвестной ролью. Потом короткий участок с 52 % ровно на трёх генах рекомбинации exo, bet и gam, потом гены регуляции — репрессор cI, cro, N, cII — на участке, богатом A и T. Снова G и C — на гене репликации P, и правый конец с генами выхода из клетки S и R. Алгоритм ничего не знает о генах — только буквы.

Сколько участков получится, решает цена смены. В участке, где G и C ровно 58 %, буква в среднем приносит верному состоянию примерно на 9 сотых бита больше, чем неверному: это расстояние Кульбака — Лейблера между двумя составами. Две смены по 2000 окупаются, если такой участок длиннее примерно 400 букв. Двигайте ползунок: при цене 200 участков 401, при 2000 — шесть, при 2900 — четыре, а начиная примерно с 3600 остаются два: левая часть генома с головкой и хвостом и правая. Под лупой — 48 букв решётки. У каждой буквы два столбика, лучшие счета путей, которые кончаются в каждом из двух состояний, и стрелки — откуда пришёл каждый путь. На границе участков один из путей приходит наискосок, из другого ряда, заплатив за смену.

Так размечал геном Гэри Черчилл в 1989 году. Он описал ДНК скрытой марковской моделью, в которой состояния задают состав букв, оценил её параметры по самой последовательности и применил к митохондриальной ДНК дрожжей, человека и мыши и к полному геному фага λ. Скрытые состояния он восстанавливал не путём Витерби, а сглаживанием — вероятностью каждого состояния у каждой буквы; о разнице — в вариантах. Через восемь лет скрытые модели стали искать гены: программа VEIL Хендерсона, Зальцберга и Фасмана, «Viterbi Exon-Intron Locator», размечала ДНК на экзоны и интроны путём Витерби.

Сообщение через шум

Сам Витерби придумал алгоритм для другой задачи — декодирования свёрточных кодов. Канал связи иногда переворачивает биты. Чтобы ошибки можно было исправить, передатчик посылает больше битов, чем в сообщении. В кодере регистр из K битов. Очередной бит сообщения вдвигается в него сверху, самый старый выпадает, и в канал уходят два бита — чётности двух наборов битов регистра. У кода с K = 7 наборы задают образующие 171 и 133, записанные восьмеричными числами: 171 = 1111001₂ берёт новый бит и биты, вдвинутые на 1, 2, 3 и 6 шагов раньше, 133 = 1011011₂ — новый и на 2, 3, 5 и 6 шагов раньше. После сообщения в регистр вдвигают шесть нулей, чтобы кодер вернулся в ноль.

Здесь скрытое состояние — шесть последних битов сообщения, их 64. Наблюдения — пары принятых битов. Если канал переворачивает каждый бит независимо с вероятностью p < 1/2, то вероятность принять пару — (1 − p)^2−kp^k, где k — число расхождений с тем, что выдал кодер, и самое вероятное сообщение — то, чей код расходится с принятым в наименьшем числе мест. Очки — число расхождений, лучший путь — с наименьшей суммой. Из каждого состояния выходят две стрелки, по биту 0 и по биту 1, в каждое приходят две, и шаг декодера — 64 одинаковые операции «сложить — сравнить — выбрать»: прибавить расхождения к счёту каждого из двух предшественников, сравнить, оставить лучший.

Листинг ниже передаёт фразу из «Войны и мира» «Пока есть жизнь, есть и счастье.» — 456 битов в UTF-8 — через канал, который переворачивает 4 % битов. Без кода перевёрнуто 19 битов, и половина букв превращается в мусор. С кодом в канал уходит 924 бита, перевёрнуто 34, а Витерби находит сообщение, чей код отличается от принятого ровно в 34 местах, — отправленное, без единой ошибки.

Почему код справляется, видно из его свободного расстояния: коды любых двух разных сообщений расходятся хотя бы в 10 битах. Поэтому, если перевёрнуто не больше четырёх битов, отправленное сообщение наверняка ближе к принятому, чем любое другое: до любого другого не меньше 10 − 4 = 6 расхождений. Когда ошибки редкие и разбросаны, то же рассуждение работает на каждом коротком отрезке. Когда их много и они стоят кучно, декодер находит ближайшее сообщение, но оно может оказаться чужим, — это видно во второй задаче ниже.

Код с K = 7 и образующими 171 и 133 с 1970-х годов — стандарт связи с космическими аппаратами: он летал на «Вояджерах» и до сих пор записан в рекомендациях CCSDS по телеметрии. Тот же код — образующие там перечислены в другом порядке, 133 и 171, — работает в Wi-Fi стандартов 802.11a и 802.11g.

Код

Листинг делает обе вещи: размечает геном фага λ при трёх ценах смены и передаёт фразу через шумный канал. Функция viterbi — общий случай с любым числом состояний и одинаковой ценой любой смены; decode — та же решётка, где из каждого состояния ровно две стрелки, а очки — расхождения битов.

При цене смены 1000 участков одиннадцать, при 2000 — шесть, при 4000 — два, и счёт лучшего пути падает: дорогие смены ограничивают выбор. Список участков при цене 2000 — тот же, что в сцене. Во второй части видно, что Витерби не знает, какие биты перевернул канал: он находит сообщение, ближайшее к принятому, а то, что оно совпало с отправленным, — заслуга кода.

Сколько это стоит

На каждом из T шагов для каждого из S состояний перебираются все предшественники — время O(T · S²). Если из каждого состояния ведёт мало переходов, считают только их: O(T · E), где E — число переходов за шаг. У свёрточного кода их 2 · 64 = 128, у генома — четыре на букву. Для одного счёта лучшего пути хватает двух столбцов таблицы, O(S) памяти, а для обратного хода нужны стрелки — O(T · S). Их можно хранить экономно: у двух состояний генома стрелка — один бит, у свёрточного кода тоже один бит, тот, что выпал из регистра.

Декодер в приёмнике не может ждать конца передачи. Он пользуется тем, что лучшие пути ко всем 64 состояниям, если пройти по ним назад, почти наверняка сходятся в одну цепочку за несколько K шагов, и решает о бите с таким опозданием. Окно обычно берут около пяти K; оно ограничивает и память, и задержку.

В привязке трека дорого не само сравнение, а очки переходов: для каждой пары кандидатов нужен путь по дорогам. Поэтому поиск обрывают на разумной длине, и одного запуска Дейкстры от кандидата хватает на всех кандидатов следующей отметки. На маршруте от вокзала это 703 поиска на 3382 пары, и вся привязка занимает миллисекунды — сколько именно у вас, сцена показывает под картой.

Аппаратные декодеры делают все «сложить — сравнить — выбрать» шага одновременно. Уже в 1971 году декодер на 64 состояния работал со скоростью 2 Мбит/с. В начале 1990-х для связи с дальним космосом собрали «Большой декодер Витерби» на 2¹⁴ состояний: 8192 параллельных узла «сложить — сравнить — выбрать», около 1 Мбит/с.

Варианты

Прямой алгоритм. Если в формуле складывать вероятности всех путей в клетку, а не выбирать лучший, получится вероятность наблюдений по всем цепочкам сразу. Ею сравнивают модели: какое семейство белков, какое слово из словаря лучше объясняет наблюдения.

Сглаживание. Прямой проход с начала и такой же с конца дают для каждого шага вероятность каждого состояния при всех наблюдениях. Самое вероятное состояние на каждом шаге отдельно — не то же самое, что самая вероятная цепочка: отдельные ответы могут сложиться в цепочку с невозможным переходом. Зато сглаживание говорит, насколько модель уверена в каждом месте. Так размечал геномы Черчилл.

Мягкие решения. Приёмник видит не биты, а уровень сигнала. Если вместо расхождений битов считать очки по самому уровню — насколько уверенно пришёл ноль или единица, — декодер ошибается заметно реже при том же шуме. Витерби с мягкими решениями работает без изменений: меняются только очки.

Лучи. Когда состояний миллионы, как в распознавании речи, на каждом шаге оставляют только лучшие — те, чей счёт не слишком отстаёт от максимума. Ответ уже не гарантирован, зато время ограничено.

Обучение. Вероятности модели подбирают по данным. Алгоритм Баума — Велша — частный случай EM — уточняет их по сглаженным вероятностям. Проще «обучение Витерби»: разметить данные лучшим путём, пересчитать вероятности по разметке и повторить.

История

Эндрю Витерби вспоминал, что алгоритм вырос из преподавания: теория информации давалась студентам трудно, и он искал инструменты попроще. Статья написана в марте 1966 года и вышла в апреле 1967-го в IEEE Transactions on Information Theory. Подробности алгоритма его отговаривали публиковать, и включить их убедил один из рецензентов, Джим Месси. В названии статьи алгоритм назван асимптотически оптимальным, в тексте — «явно неоптимальным», а в выводах сказано, что из-за огромной памяти он непрактичен, но помогает понять свёрточные коды.

Вскоре Дэвид Форни нарисовал, как сливаются пути свёрточного кода, — и назвал картинку решёткой. Оттуда оставался небольшой шаг: алгоритм Витерби — точный способ найти кратчайший путь по решётке, то есть оптимальный декодер. Форни вспоминает, что позвонил Витерби и сказал, что тот был слишком скромен. Джим Омура в статье, отправленной в мае 1968 года, показал, что это обычное прямое динамическое программирование. Позже Форни писал, что любопытно гадать, сложилась бы история иначе, если бы алгоритм с самого начала назывался динамическим программированием.

Практичность увидел Джерри Хеллер: в 1968–1969 годах он смоделировал коды с короткой длиной регистра и нашёл, что код всего на 64 состояния даёт выигрыш порядка 6 дБ. В 1971 году готов был декодер на 2 Мбит/с, а в 1970-х алгоритм вошёл в стандарт связи с дальним космосом — позже в паре с внешним кодом Рида — Соломона. В 1973 году Форни написал обзор, главная мысль которого — алгоритм годится для любой задачи, где нужно восстановить выход конечного автомата по наблюдениям с шумом, то есть для любой скрытой марковской модели.

Сами скрытые марковские модели описали Леонард Баум и Тед Петри в 1966 году. В 1970–1980-х на них и на алгоритме Витерби построили распознавание речи, а в 1989 году Лоуренс Рабинер собрал этот опыт в подробный обзор. Потом пришли биология — Черчилл в 1989 году, поиск генов в 1990-х — и привязка треков: Ньюсон и Крумм опубликовали свою модель в 2009 году вместе с данными поездки, чтобы другие могли сравнивать свои алгоритмы на тех же данных.

Задачи