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

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

Версия от 10:54, 16 марта 2013; Kropotov (Обсуждение | вклад)
(разн.) ← Предыдущая | Текущая версия (разн.) | Следующая → (разн.)
Перейти к: навигация, поиск

Содержание

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

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

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

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

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

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

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

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

Hx=0,

где H\in\{0,1\}^{(N-K){\times}N} — матрица ранга N-K. Такая матрица называется проверочной матрицей кода, т.к. с её помощью можно проверить, является ли слово x кодовым словом путём проверки соотношения Hx=0 (здесь и далее все операции проводятся по модулю 2).

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

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

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

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

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

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

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

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

p(y|x) = \prod_{n=1}^Np(y_n|x_n) = \prod_{n=1}^Nq^{y_n+x_n}(1-q)^{y_n+x_n+1}.

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

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

p(y,x)\propto p(y|x)I[Hx=0] = \prod_{n=1}^Np(y_n|x_n)\prod_{m=1}^MI[h_m^Tx = 0].

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

Процесс декодирования, т.е. восстановление кодового слова x по полученному слову y, предлагается осуществлять как x_n^* = \arg\max p(x_n|y), а маргинальные распределения p(x_n|y) оценивать в помощью циклического алгоритма передачи сообщений (sum-product loopy BP) на фактор-графе. При этом для упрощения алгоритма декодирования предлагается избавиться от факторов-унарных потенциалов (путем их включения в сообщения от переменных x_n к факторам f_m), а в качестве расписания пересчёта сообщений выбрать параллельное расписание, при котором сначала все переменные одновременно посылают сообщения во все факторы, а затем все факторы одновременно посылают сообщения во все вершины.

Введём обозначения N(n)=\{m:H_{mn}=1\} — множество факторов, в которых участвует переменная x_n, и N(m)=\{n:H_{mn}=1\} — множество переменных, которые входят в фактор f_m. Тогда общая схема алгоритма декодирования выглядит следующим образом:

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

\mu_{x_n\rightarrow f_m}(x_n) = p(y_n|x_n);

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

\mu_{f_m\rightarrow x_n}(x_n) = \sum_{x_{n^'}:n^'\in N(m)\backslash n}I[x_n+\sum_{n^'}x_{n^'}=0]\prod_{n^'\in N(m)\backslash n}\mu_{x_{n^'}\rightarrow f_m}(x_{n^'});

3. Пересчет сообщений от переменных:

\mu_{x_n\rightarrow f_m}(x_n) \propto p(y_n|x_n)\prod_{m^'\in N(n)\backslash m}\mu_{f_{m^'}\rightarrow x_n}(x_n);
\hat{p}_n(x_n|y)\propto p(y_n|x_n)\prod_{m^'\in N(n)}\mu_{f_{m^'}\rightarrow x_n}(x_n);
Символом \propto обозначается пропорциональность. Таким образом, при пересчете все сообщения от переменных и оценки на маргинальные распределения должны нормироваться так, чтобы \mu_{x_n\rightarrow f_m}(0)+\mu_{x_n\rightarrow f_m}(1)=1 и \hat{p}_n(0|y)+\hat{p}_n(1|y)=1;

4. Оценка кодового слова:

\hat{x}_n=\arg\max\hat{p}_n(x_n|y);

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

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

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

Рассмотрим более эффективную схему реализации шага 2 путем его сведения к задаче вывода в графической модели с графом-цепочкой. Пусть нам необходимо вычислить сообщение \mu_{f_m\rightarrow x_n}(x_n). Перенумеруем все переменные, входящие в m-ый фактор (кроме переменной x_n), как x_1,x_2,\dots,x_{N(m)-1}. Рассмотрим графическую модель с N(m)-1 бинарными переменными s_i и графом

.

Положим, что s_i=\sum_{j=1}^ix_j. Определим априорное распределение и вероятность перехода как

p(s_1) = \mu_{x_1\rightarrow f_m}(s_1); p(s_i|s_{i-1}) = \begin{cases}\mu_{x_i\rightarrow f_m}(0),& s_i=s_{i-1},\\ \mu_{x_i\rightarrow f_m}(1),& s_i\neq s_{i-1}.\end{cases}

Тогда требуемое сообщение \mu_{f_m\rightarrow x_n}(x_n) поэлементно равно маргинальному распределению p(s_{N(m)-1}), которое может быть эффективно вычислено путем однократной пересылки сообщений вдоль цепи от s_1 до s_{N(m)-1}.

Дальнейшее повышение эффективности реализации шага 2 связано с рассмотрением разностей значений сообщений \delta\mu_{x_n\rightarrow f_m} = \mu_{x_n\rightarrow f_m}(0) - \mu_{x_n\rightarrow f_m}(1) и \delta\mu_{f_m\rightarrow x_n} = \mu_{f_m\rightarrow x_n}(0) - \mu_{f_m\rightarrow x_n}(1). В терминах разностей можно показать, что новое значение \delta\mu_{f_m\rightarrow x_n} может быть вычислено как

\delta\mu_{f_m\rightarrow x_n} = \prod_{n^'\in N(m)\backslash n}\delta\mu_{x_{n^'}\rightarrow f_m}.   (*)

Зная значение \delta\mu_{f_m\rightarrow x_n} и учитывая условие нормировки на сообщения, сами сообщения от факторов могут быть вычислены как

\mu_{f_m\rightarrow x_n}(0) = \frac{1}{2}(1+\delta\mu_{f_m\rightarrow x_n}); \mu_{f_m\rightarrow x_n}(1) = \frac{1}{2}(1-\delta\mu_{f_m\rightarrow x_n}).

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

Покажем теперь справедливость условия (*). Для этого рассмотрим две произвольные бинарные независимые случайные величины u_1 и u_2. Обозначим через p_i^u величину \mathbb{P}\{u_i=u\}. Тогда очевидно, что

\mathbb{P}\{u_1+u_2=0\} = p_1^0p_2^0 + p_1^1p_2^1; \mathbb{P}\{u_1+u_2=1\} = p_1^0p_2^1 + p_1^1p_2^0; \mathbb{P}\{u_1+u_2=0\}-\mathbb{P}\{u_1+u_2=1\} = (p_1^0-p_1^1)(p_2^0-p_2^1) = \delta p_1\delta p_2.

Аналогичные рассуждения справедливы для произвольного количества случайных величин u. В частности,

\mathbb{P}\{u_1+u_2+u_3=0\} - \mathbb{P}\{u_1+u_2+u_3=1\} = \delta p_1\delta p_2\delta p_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.
Личные инструменты