Графические модели (курс лекций)/2014/Задание 2

Материал из MachineLearning.

(Различия между версиями)
Перейти к: навигация, поиск
(Оптимизация в процедуре декодирования)
(Оптимизации в процедуре декодирования)
Строка 88: Строка 88:
==== Расписание пересчёта сообщений и дэмпфирование ====
==== Расписание пересчёта сообщений и дэмпфирование ====
 +
 +
Общая схема циклического алгоритма передачи сообщений оставляет определённый произвол в выборе расписания пересчёта сообщений. Обычно здесь рассматриваются следующие подходы:
 +
* ''Параллельное расписание''. В данном случае сначала все вершины посылают сообщения во все факторы, а затем все факторы посылают сообщения во все вершины.
 +
* ''Последовательное расписание''. Здесь выбирается некоторый (например, случайный) порядок последовательного пересчёта всех сообщений (и от вершин к факторам, и от факторов к вершинам). При этом данный порядок может меняться от итерации к итерации.
 +
 +
При использовании дэмпфирования с параметром <tex>\lambda\in(0,1]</tex> сообщения на итерации <tex>t</tex> пересчитываются как
 +
:<tex>\mu^t_{h_j\rightarrow e_i}(e_i) = \lambda\mu_{h_j\rightarrow e_i}(e_i) + (1-\lambda)\mu^{t-1}_{h_j\rightarrow e_i}(e_i)</tex>;
 +
:<tex>\mu^t_{e_i\rightarrow h_j}(e_i) = \lambda\mu_{e_i\rightarrow h_j}(e_i) + (1-\lambda)\mu^{t-1}_{e_i\rightarrow h_j}(e_i)</tex>.
 +
 +
Здесь <tex>\mu_{h_j\rightarrow e_i}(e_i)</tex> и <tex>\mu_{e_i\rightarrow h_j}(e_i)</tex> — сообщения, вычисляемые на шагах 2 и 3 общего алгоритма декодирования.
== Формулировка задания ==
== Формулировка задания ==

Версия 18:51, 3 марта 2014

Внимание! Текст задания находится в стадии разработки. Просьба не приступать к выполнению задания до тех пор, пока это предупреждение не будет удалено.


Содержание

Демонстрация процесса кодирования/декодирования с использованием низкоплотностного кода. Кодируемое сообщение (изображение цифры) содержит значительную избыточность для возможности визуальной оценки процесса декодирования. Алгоритм декодирования не использует свойство избыточности входного сигнала.
Демонстрация процесса кодирования/декодирования с использованием низкоплотностного кода. Кодируемое сообщение (изображение цифры) содержит значительную избыточность для возможности визуальной оценки процесса декодирования. Алгоритм декодирования не использует свойство избыточности входного сигнала.

Начало выполнения задания: 4 марта 2014 г.
Срок сдачи: 18 марта 2014 г., 23:59.

Среда для выполнения задания — MATLAB.

Низкоплотностные коды

Задача помехоустойчивого кодирования

Рассмотрим решение задачи безошибочной передачи потока битовой информации по каналу с шумом с помощью кодов, исправляющих ошибки. При блоковом кодировании входящий поток информации разбивается на блоки фиксированной длины k, и каждый блок кодируется/декодируется независимо. Обозначим один такой блок через u\in\{0,1\}^k. Предположим, что во входном потоке данных, вообще говоря, нет избыточности. Поэтому для реализации схемы, способной исправлять ошибки, необходимо закодировать блок u в некоторое кодовое слово большей длины путем добавления избыточности в передаваемые данные. Обозначим кодовое слово через v\in\{0,1\}^n, n>k. Для кодирования всевозможных блоков u необходимо использовать 2^k кодовых слов длины n. Назовём множество 2^k кодовых слов длины n (n,k)-блоковым кодом, а величину r=k/nскоростью кода. При передаче по каналу с шумом кодовое слово v превращается в принятое слово w, которое, вообще говоря, отличается от v. Предположим, что при передаче по каналу длина сообщения не изменяется, т.е. w\in\{0,1\}^n, а происходит лишь инверсия некоторых бит. Задача алгоритма декодирования состоит в восстановлении по w переданного слова v (например, путем поиска среди всевозможных кодовых слов ближайшего к w). Обозначим результат работы алгоритма декодирования через \hat{v}. На последнем этапе декодированное слово \hat{v} переводится в декодированное слово исходного сообщения \hat{u}.

Кодирование с помощью (n,k)-линейного блокового кода

Множество \{0,1\}^n с операциями суммы и произведения по модулю 2 образует линейное пространство над конечным полем из двух элементов \{0,1\}. (n,k)-блоковый код называется линейным, если множество его кодовых слов образует линейное подпространство размерности k общего линейного пространства \{0,1\}^n. Одним из способов задания k-мерного линейного подпространства является рассмотрение множества решений следующей системы линейных уравнений:

Hv=0,

где H\in\{0,1\}^{(n-k){\times}n} — матрица ранга m=n-k (эта матрица задаёт базис линейного подпространства, ортогонального к рассматриваемому (n,k)-коду). Матрица H называется проверочной матрицей кода, т.к. с её помощью можно проверить, является ли слово v кодовым словом путём проверки соотношения Hv=0 (здесь и далее все операции проводятся по модулю 2).

Рассмотрим задачу кодирования слов исходного сообщения u в кодовые слова v (n,k)-линейного блокового кода, заданного своей проверочной матрицей H. Для этого можно найти базис k-мерного линейного подпространства g_1,\dots,g_k\in\{0,1\}^n. Тогда, рассматривая базисные вектора как столбцы общей матрицы G\in\{0,1\}^{n{\times}k}, операция кодирования может быть представлена как v=Gu. Матрица G называется порождающей матрицей кода. Кодирование называется систематическим, если все биты слова u копируются в некоторые биты кодового слова v, т.е. в матрице G некоторое подмножество строк образует единичную матрицу размера k{\times}k. При систематическом кодировании обратный процесс преобразования из декодированного кодового слова \hat{v} в декодированное сообщение \hat{u} становится тривиальным.

Одним из способов построения порождающей матрицы кода по заданной проверочной матрице является преобразование проверочной матрицы к каноническому ступенчатому виду. Такое преобразование всегда может быть сделано с помощью гауссовских исключений. С точностью до перестановки столбцов канонический ступенчатый вид матрицы H эквивалентен её представлению в виде \begin{bmatrix} I_{m} & P\end{bmatrix}, где I_{m} — единичная матрица размера m{\times}m. Тогда в качестве порождающей матрицы, обеспечивающей систематическое кодирование, можно выбрать матрицу

G = \begin{bmatrix}P\\ I_k\end{bmatrix}.

Действительно, в этом случае HGu = (P+P)u = 0.

Декодирование низкоплотностного кода

Фактор-граф для (16,4)-низкоплотностного кода, в проверочной матрице которого в каждом столбце по 3 единицы, а в каждой строке - по 4 единицы.
Фактор-граф для (16,4)-низкоплотностного кода, в проверочной матрице которого в каждом столбце по 3 единицы, а в каждой строке - по 4 единицы.

Низкоплотностным кодом (или кодом с малой плотностью проверок на чётность) называется бинарный (n,k)-линейный блоковый код, в котором проверочная матрица H\in\{0,1\}^{m{\times}n} является сильно разреженной.

Рассмотрим бинарный симметричный канал для передачи данных. Здесь при передаче каждый бит независимо инвертируется с некоторой вероятностью q. В результате бинарный симметричный канал задает распределение p(w|v) для передаваемого кодового слова v\in\{0,1\}^n и полученного на выходе слова w\in\{0,1\}^n как

p(w|v) = \prod_{i=1}^np(w_i|v_i) = \prod_{i=1}^nq^{w_i+v_i}(1-q)^{w_i+v_i+1}.

Пропускная способность данного канала определяется величиной 1+q\log_2q+(1-q)\log_2(1-q).

Объединяя низкоплотностный код с бинарным симметричным каналом, получаем следующую вероятностную модель для пары w,v:

p(v|w)\propto p(w|v)I[Hv=0] \propto \prod_{i=1}^np(w_i|v_i)\prod_{j=1}^mI[h_j^Tv = 0].

Здесь m=n-k — количество проверок на чётность, h_j^T — j-ая строка матрицы H, а I[\cdot] — индикаторная функция. Фактор-граф введённой модели показан на рис. справа.

Назовём синдромом принятого слова w вектор s\in\{0,1\}^m, определяемый как s=Hw. Процесс передачи кодового слова по бинарному симметричному каналу можно представить как w=v+e, где e\in\{0,1\}^n — вектор ошибок (e_i=1, если в позиции i произошла ошибка). Тогда s = Hw = H(v+e) = Hv+He = He. Далее можно перейти от вероятностной модели для переменных v,w к аналогичной для переменных e,s:

p(e|s)\propto p(e)p(s|e) \propto\prod_{i=1}^np(e_i)\prod_{j=1}^mI[h_j^Te=s_j].

Здесь p(e_i)=q^{e_i}(1-q)^{1-e_i}. Зная значение вектора ошибок \hat{e}, результат декодирования вычисляется как \hat{v}=w+\hat{e}. При использовании вероятностной модели для e,s тестирование алгоритма декодирования можно проводить без предварительной реализации процедуры кодирования.

При использовании побитовой функции потерь \lambda(e,\tilde{e}) = \sum_{i=1}^nI[e_i\neq \tilde{e}_i] оптимальная процедура декодирования связана с максимизацией маргиналов отдельных переменных, т.е. \hat{e}_n = \arg\max p(e_n|s). Для поиска маргинальных распределений p(e_n|s) воспользуемся циклическим алгоритмом передачи сообщений (sum-product loopy BP) на фактор-графе. Введём обозначения N(i)=\{j:H_{ji}=1\} — множество факторов, в которых участвует переменная e_i, и N(j)=\{i:H_{ji}=1\} — множество переменных, которые входят в фактор h_j. Тогда общая схема алгоритма декодирования выглядит следующим образом:

1. Инициализация:

\mu_{e_i\rightarrow h_j}(e_i) = p(e_i);

2. Пересчет сообщений от факторов:

\mu_{h_j\rightarrow e_i}(e_i) = \sum_{e_k:k\in N(j)\backslash i}I[e_i+\sum_{k}e_{k}=s_j]\prod_{k\in N(j)\backslash i}\mu_{e_{k}\rightarrow h_j}(e_{k});

3. Пересчет сообщений от переменных и вычисление beliefs (оценок на маргинальные распределения):

\mu_{e_i\rightarrow h_j}(e_i) \propto p(e_i)\prod_{k\in N(i)\backslash j}\mu_{h_{k}\rightarrow e_i}(e_i);
b_i(e_i)\propto p(e_i)\prod_{k\in N(i)}\mu_{h_{k}\rightarrow e_i}(e_i);

4. Оценка вектора ошибок:

\hat{e}_i=\arg\max b_i(e_i);

5. Критерий остановки:

Если H\hat{e}=s, то выход алгоритма со статусом 0;
Если достигнуто максимальное число итераций или произошла стабилизация по всем b_i, то выход алгоритма со статусом -1;
Переход к шагу 2.

Оптимизации в процедуре декодирования

Пересчёт сообщений от факторов

При прямой реализации шага 2 общего алгоритма декодирования требуется рассмотрение для каждого фактора 2^{|N(j)|-1} различных конфигураций значений переменных e_k. Это может приводить как к низкой скорости пересчета сообщений, так и к большим требованиям по памяти.

Рассмотрим более эффективную схему реализации шага 2. Пусть нам необходимо вычислить сообщение \mu_{h_j\rightarrow e_i}(e_i). Перенумеруем все переменные, входящие в j-ый фактор (кроме переменной e_i), как e_{(1)},e_{(2)},\dots,e_{(l)}, где l=|N(j)|-1. Тогда вычисление сообщения \mu_{h_j\rightarrow e_i}(e_i) можно записать как

\mu_{h_j\rightarrow e_i}(e_i)\propto\sum_{e_{(k)}}I[\sum_{k=1}^le_{(k)} = s_j+e_i]\prod_{k=1}^l\mu_{e_{(k)}\rightarrow h_j}(e_{(k)}).

Данный результат можно интерпретировать как вычисление вероятности p(\sum_{k=1}^le_{(k)}=s_j+e_i) для набора независимых бинарных переменных e_{(k)} с распределением p_{(k)}(e_{(k)}) = \mu_{e_{(k)}\rightarrow h_j}(e_{(k)}). Обозначим через e^k сумму первых k переменных e_{(k)}, т.е. e^k=e_{(1)}+e_{(2)}+\dots+e_{(k)}. Тогда распределение на e^k можно итерационно пересчитывать по формуле

p_k(e^k) = \sum_{e^{k-1}}p_{k-1}(e^{k-1})p_{(k)}(e^k-e^{k-1}).

В результате вычисление p_l(e^l) имеет линейную по l сложность.

Дальнейшее повышение эффективности реализации шага 2 связано с рассмотрением разностей вероятностей \delta p_{(k)} = p_{(k)}(0) - p_{(k)}(1) и \delta p_k = p_k(0)-p_k(1). Нетрудно показать, что

\delta p_l = \prod_{k=1}^l\delta p_{(k)}. (*)

Зная значение \delta p_l и учитывая условие нормировки p_l(0)+p_l(1)=1, сами вероятности могут быть вычислены как p_l(0)=\frac{1}{2}(1+\delta p_l), p_l(1)=\frac{1}{2}(1-\delta p_l). Таким образом, \mu_{h_j\rightarrow e_i}(e_i) = p_l(s_j+e_i). Основное преимущество условия (*) по сравнению с последовательным пересчётом распределений p_k(e^k) связано с тем, что формула (*) может быть реализована с помощью векторных операций в MATLAB.

Расписание пересчёта сообщений и дэмпфирование

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

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

При использовании дэмпфирования с параметром \lambda\in(0,1] сообщения на итерации t пересчитываются как

\mu^t_{h_j\rightarrow e_i}(e_i) = \lambda\mu_{h_j\rightarrow e_i}(e_i) + (1-\lambda)\mu^{t-1}_{h_j\rightarrow e_i}(e_i);
\mu^t_{e_i\rightarrow h_j}(e_i) = \lambda\mu_{e_i\rightarrow h_j}(e_i) + (1-\lambda)\mu^{t-1}_{e_i\rightarrow h_j}(e_i).

Здесь \mu_{h_j\rightarrow e_i}(e_i) и \mu_{e_i\rightarrow h_j}(e_i) — сообщения, вычисляемые на шагах 2 и 3 общего алгоритма декодирования.

Формулировка задания

  1. Реализовать алгоритм построения по заданной проверочной матрице чётности H порождающей матрицы кода G для систематического кодирования;
  2. Реализовать алгоритм декодирования низкоплотностного кода на основе loopy BP; при реализации шага 2 пересчета сообщений от факторов к переменным необходимо использовать эффективные схемы, обозначенные выше; при реализации на MATLAB одной итерации схемы передачи сообщений использование вложенных циклов является нежелательным; провести временные замеры реализованного алгоритма для различных значений входных параметров;
  3. Рассмотрим две характеристики качества кода — вероятность совершить ошибку хотя бы в одном бите при декодировании блока (p(\hat{x}\neq x)) и среднюю вероятность совершить ошибку при декодировании в одном бите (\frac{1}{N}\sum_{n=1}^Np(\hat{x}_n\neq x_n)). Требуется реализовать алгоритм оценки вероятности битовой и блоковой ошибки кода с помощью метода стат. испытаний (многократная случайная генерация слова t, его преобразование к кодовому слову x, передача по каналу с независимым инвертированием каждого бита с заданной вероятностью q, восстановление кодового слова \hat{x} с помощью алгоритма декодирования и подсчет необходимых характеристик);
  4. Провести эксперименты по оцениванию битовой и блоковой ошибки низкоплотностного кода для различных значений длины кодового слова N, скорости кода R, вероятности инвертирования бита при передаче по каналу связи q и среднего количества единиц в столбце проверочной матрицы j. В частности, необходимо проанализировать следующие ситуации:
    • Теорема Шеннона определяет пропускную способность канала как максимально допустимую скорость кода, при которой возможно осуществление надежной коммуникации. Требуется проверить, как меняются характеристики кода при изменении скорости R от минимального значения до пропускной способности канала.
    • Теорема Шеннона предполагает, что качество кода растет при увеличении длины кодового слова N. Требуется проверить это предположение.
    • Одно из следствий теоремы Шеннона утверждает, что хорошими кодами являются коды со случайной проверочной матрицей H. В частности, здесь предполагается, что качество кода должно расти при увеличении среднего количества единиц в столбце проверочной матрицы j. Требуется проверить это утверждение для низкоплотностных кодов (путём рассмотрения нескольких значений j, начиная от 3).
  5. Составить отчёт в формате PDF с описанием всех проведённых исследований. Данный отчёт обязательно должен содержать описание отладочных тестов, которые проводились для проверки корректности реализованных алгоритмов. Также в этом отчёте должны быть графики поведения характеристик кода в зависимости от значений различных параметров.

Рекомендации по выполнению задания

1. Разреженную проверочную матрицу кода заданных размеров можно строить с помощью случайной генерации (с соблюдением условия полноранговости). Однако, здесь рекомендуется воспользоваться готовой реализацией, которая генерирует проверочную матрицу заданного размера с заданным количеством единиц в каждом столбце. При этом реализованный алгоритм старается сократить количество циклов длины 4 в генерируемой матрице, т.к. наличие коротких циклов в графе, как правило, значительно усложняет работу алгоритма loopy BP.

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

3. Алгоритм декодирования удобно реализовывать в т.н. синдромном представлении. Полученное на выходе канала слово y можно представить как x+e, где x — посылаемое кодовое слово, а e\in\{0,1\}^N — вектор ошибок, которые вносит канал. Назовём величину z = Hy синдромом полученного сообщения. Тогда z=Hy=H(x+e)=He, и на этапе декодирования можно вместо x оценивать вектор ошибок e путем оценки аргмаксимумов маргинальных распределений p(e_n|z) в вероятностной модели

p(e,z)\propto p(e)I[He=z]=\prod_{n=1}^Np(e_n)\prod_{m=1}^MI[h_m^Te=z_m].

Для бинарного симметричного канала p(e_n) = q^{e_n}(1-q)^{e_n+1}. Зная e, искомое кодовое слово можно найти как x=y+e. Заметим, что реализация алгоритма декодирования в синдромном представлении избавляет от необходимости предварительной реализации алгоритма кодирования.

4. При тестировании алгоритма декодирования рекомендуется пробовать запускать итерационный процесс из различных начальных приближений для сообщений (не только из значений унарных потенциалов). При корректной реализации метода должна наблюдаться сходимость к одним и тем же финальным сообщениям независимо от начальных приближений.

Оформление задания

Выполненное задание следует отправить письмом по адресу bayesml@gmail.com с заголовком письма «[ГМ13] Задание 2 <ФИО>». Убедительная просьба присылать выполненное задание только один раз с окончательным вариантом. Также убедительная просьба строго придерживаться заданных ниже прототипов реализуемых функций.

Присланный вариант задания должен содержать в себе:

  • Файл отчёта в формате PDF с указанием ФИО.
  • Все исходные коды с необходимыми комментариями.

 

Построение порождающей матрицы для систематического кодирования
[G, ind] = ldpc_gen_matrix(H)
ВХОД
H — проверочная матрица чётности, бинарная матрица размера MxN;
ВЫХОД
G — порождающая матрица кода, бинарная матрица размера Nx(N-M);
ind — номера позиций кодового слова, в которые копируются биты исходного сообщения, т.е. G(ind, :) является единичной матрицей.

 

Алгоритм декодирования LDPC-кода в синдромном представлении
[e, status] = ldpc_decoding(z, H, q, param_name1, param_value1, ...)
ВХОД
z — наблюдаемый синдром, бинарный вектор-столбец длины M;
H — проверочная матрица чётности, бинарная матрица размера MxN;
q — вероятность инверсии бита при передаче по каналу связи, число от 0 до 0.5;
(param_name, param_value) — набор необязательных параметров алгоритма, следующие имена и значения возможны:
'max_iter' — максимальное число итераций алгоритма декодирования, число, по умолчанию = 200;
'eps' — порог стабилизации для сообщений, число, по умолчанию = 1e-4;
'display' — режим отображения, true или false, если true, то отображается промежуточная информация на итерациях, например, номер итерации, текущее число ошибок декодирования, невязка для сообщений и т.д.
ВЫХОД
e — восстановленный вектор ошибок, бинарный вектор-столбец длины N;
status — результат декодирования, равен 0, если вектор e восстановлен без ошибок, равен -1, если произошел выход по максимальному числу итераций или стабилизации значений сообщений.

 

Оценка характеристик LDPC-кода с помощью метода Монте Карло
[err_bit, err_block, diver] = ldpc_mc(H, G, q, num_points)
ВХОД
H — проверочная матрица чётности, бинарная матрица размера MxN;
G — порождающая матрица кода, бинарная матрица размера Nx(N-M);
q — вероятность инверсии бита при передаче по каналу связи, число от 0 до 0.5;
num_points — общее количество экспериментов, число;
ВЫХОД
err_bit — вероятность битовой ошибки декодирования (относительно N бит кодового слова), число от 0 до 1;
err_block — вероятность блоковой ошибки декодирования, число от 0 до 1;
diver — доля ситуаций расходимости алгоритма декодирования, число от 0 до 1.
Личные инструменты