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

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

(Различия между версиями)
Перейти к: навигация, поиск
(Задача помехоустойчивого кодирования)
м
 
(11 промежуточных версий не показаны.)
Строка 1: Строка 1:
-
{{stop|Внимание! Текст задания находится в стадии разработки. Просьба не приступать к выполнению задания до тех пор, пока это предупреждение не будет удалено.}}
 
-
 
{{Main|Графические модели (курс лекций)}}
{{Main|Графические модели (курс лекций)}}
Строка 8: Строка 6:
|}
|}
-
'''Начало выполнения задания''': 4 марта 2014 г.<br>
+
'''Начало выполнения задания''': 5 марта 2014 г.<br>
'''Срок сдачи''': {{важно|18 марта 2014 г., 23:59.}}
'''Срок сдачи''': {{важно|18 марта 2014 г., 23:59.}}
Строка 19: Строка 17:
Рассмотрим решение задачи безошибочной передачи потока битовой информации по каналу с шумом с помощью кодов, исправляющих ошибки. При ''блоковом'' кодировании входящий поток информации разбивается на блоки фиксированной длины <tex>k</tex>, и каждый блок кодируется/декодируется независимо. Обозначим один такой блок через <tex>u\in\{0,1\}^k</tex>. Предположим, что во входном потоке данных, вообще говоря, нет избыточности. Поэтому для реализации схемы, способной исправлять ошибки, необходимо закодировать блок <tex>u</tex> в некоторое кодовое слово большей длины путем добавления избыточности в передаваемые данные. Обозначим кодовое слово через <tex>v\in\{0,1\}^n</tex>, <tex>n>k</tex>. Для кодирования всевозможных блоков <tex>u</tex> необходимо использовать <tex>2^k</tex> кодовых слов длины <tex>n</tex>. Назовём множество <tex>2^k</tex> кодовых слов длины <tex>n</tex> ''(n,k)-блоковым кодом'', а величину <tex>r=k/n</tex> — ''скоростью кода''. При передаче по каналу с шумом кодовое слово <tex>v</tex> превращается в принятое слово <tex>w</tex>, которое, вообще говоря, отличается от <tex>v</tex>. Предположим, что при передаче по каналу длина сообщения не изменяется, т.е. <tex>w\in\{0,1\}^n</tex>, а происходит лишь инверсия некоторых бит. Задача алгоритма декодирования состоит в восстановлении по <tex>w</tex> переданного слова <tex>v</tex> (например, путем поиска среди всевозможных кодовых слов ближайшего к <tex>w</tex>). Обозначим результат работы алгоритма декодирования через <tex>\hat{v}</tex>. На последнем этапе декодированное слово <tex>\hat{v}</tex> переводится в декодированное слово исходного сообщения <tex>\hat{u}</tex>.
Рассмотрим решение задачи безошибочной передачи потока битовой информации по каналу с шумом с помощью кодов, исправляющих ошибки. При ''блоковом'' кодировании входящий поток информации разбивается на блоки фиксированной длины <tex>k</tex>, и каждый блок кодируется/декодируется независимо. Обозначим один такой блок через <tex>u\in\{0,1\}^k</tex>. Предположим, что во входном потоке данных, вообще говоря, нет избыточности. Поэтому для реализации схемы, способной исправлять ошибки, необходимо закодировать блок <tex>u</tex> в некоторое кодовое слово большей длины путем добавления избыточности в передаваемые данные. Обозначим кодовое слово через <tex>v\in\{0,1\}^n</tex>, <tex>n>k</tex>. Для кодирования всевозможных блоков <tex>u</tex> необходимо использовать <tex>2^k</tex> кодовых слов длины <tex>n</tex>. Назовём множество <tex>2^k</tex> кодовых слов длины <tex>n</tex> ''(n,k)-блоковым кодом'', а величину <tex>r=k/n</tex> — ''скоростью кода''. При передаче по каналу с шумом кодовое слово <tex>v</tex> превращается в принятое слово <tex>w</tex>, которое, вообще говоря, отличается от <tex>v</tex>. Предположим, что при передаче по каналу длина сообщения не изменяется, т.е. <tex>w\in\{0,1\}^n</tex>, а происходит лишь инверсия некоторых бит. Задача алгоритма декодирования состоит в восстановлении по <tex>w</tex> переданного слова <tex>v</tex> (например, путем поиска среди всевозможных кодовых слов ближайшего к <tex>w</tex>). Обозначим результат работы алгоритма декодирования через <tex>\hat{v}</tex>. На последнем этапе декодированное слово <tex>\hat{v}</tex> переводится в декодированное слово исходного сообщения <tex>\hat{u}</tex>.
-
=== Кодирование с помощью (N,K)-линейного блокового кода ===
+
=== Кодирование с помощью (n,k)-линейного блокового кода ===
-
Множество <tex>\{0,1\}^N</tex> с операциями суммы и произведения по модулю 2 образует линейное пространство над конечным полем из двух элементов <tex>\{0,1\}</tex>. (N,K)-блоковый код называется ''линейным'', если множество его кодовых слов образует линейное подпространство размерности <tex>K</tex> общего линейного пространства <tex>\{0,1\}^N</tex>. Одним из способов задания <tex>K</tex>-мерного линейного подпространства является рассмотрение множества решений следующей системы линейных уравнений:
+
Множество <tex>\{0,1\}^n</tex> с операциями суммы и произведения по модулю 2 образует линейное пространство над конечным полем из двух элементов <tex>\{0,1\}</tex>. (n,k)-блоковый код называется ''линейным'', если множество его кодовых слов образует линейное подпространство размерности <tex>k</tex> общего линейного пространства <tex>\{0,1\}^n</tex>. Одним из способов задания <tex>k</tex>-мерного линейного подпространства является рассмотрение множества решений следующей системы линейных уравнений:
-
:<tex>Hx=0</tex>,
+
:<tex>Hv=0</tex>,
-
где <tex>H\in\{0,1\}^{(N-K){\times}N}</tex> — матрица ранга <tex>N-K</tex>. Такая матрица называется ''проверочной матрицей'' кода, т.к. с её помощью можно проверить, является ли слово <tex>x</tex> кодовым словом путём проверки соотношения <tex>Hx=0</tex> (здесь и далее все операции проводятся по модулю 2).
+
где <tex>H\in\{0,1\}^{(n-k){\times}n}</tex> — матрица ранга <tex>m=n-k</tex> (эта матрица задаёт базис линейного подпространства, ортогонального к рассматриваемому (n,k)-коду). Матрица <tex>H</tex> называется ''проверочной матрицей'' кода, т.к. с её помощью можно проверить, является ли слово <tex>v</tex> кодовым словом путём проверки соотношения <tex>Hv=0</tex> (здесь и далее все операции проводятся по модулю 2).
-
Рассмотрим задачу кодирования слов исходного сообщения <tex>t</tex> в кодовые слова <tex>x</tex> (N,K)-линейного блокового кода, задаваемого проверочной матрицей <tex>H</tex>. Для этого необходимо найти базис <tex>K</tex>-мерного линейного подпространства <tex>g_1,\dots,g_K\in\{0,1\}^N</tex>. Тогда, рассматривая базисные вектора как столбцы общей матрицы <tex>G\in\{0,1\}^{N{\times}K}</tex>, операция кодирования может быть представлена как <tex>x=Gt</tex>. Матрица <tex>G</tex> называется ''порождающей матрицей'' кода. Кодирование называется ''систематическим'', если все биты слова <tex>t</tex> копируются в некоторые биты кодового слова <tex>x</tex>, т.е. в матрице G некоторое подмножество строк образует единичную матрицу размера <tex>K{\times}K</tex>. При систематическом кодировании обратный процесс преобразования из декодированного кодового слова <tex>\hat{x}</tex> в декодированное слово сообщения <tex>\hat{t}</tex> становится тривиальным.
+
Рассмотрим задачу кодирования слов исходного сообщения <tex>u</tex> в кодовые слова <tex>v</tex> (n,k)-линейного блокового кода, заданного своей проверочной матрицей <tex>H</tex>. Для этого можно найти базис <tex>k</tex>-мерного линейного подпространства <tex>g_1,\dots,g_k\in\{0,1\}^n</tex>. Тогда, рассматривая базисные вектора как столбцы общей матрицы <tex>G\in\{0,1\}^{n{\times}k}</tex>, операция кодирования может быть представлена как <tex>v=Gu</tex>. Матрица <tex>G</tex> называется ''порождающей матрицей'' кода. Кодирование называется ''систематическим'', если все биты слова <tex>u</tex> копируются в некоторые биты кодового слова <tex>v</tex>, т.е. в матрице <tex>G</tex> некоторое подмножество строк образует единичную матрицу размера <tex>k{\times}k</tex>. При систематическом кодировании обратный процесс преобразования из декодированного кодового слова <tex>\hat{v}</tex> в декодированное сообщение <tex>\hat{u}</tex> становится тривиальным.
-
Одним из способов построения порождающей матрицы кода по заданной проверочной матрице является преобразование проверочной матрицы к [http://ru.wikipedia.org/wiki/Ступенчатый_вид_по_строкам каноническому ступенчатому виду]. Такое преобразование всегда может быть сделано с помощью [http://en.wikipedia.org/wiki/Gaussian_elimination гауссовских исключений]. С точностью до перестановки столбцов канонический ступенчатый вид матрицы <tex>H</tex> эквивалентен её представлению в виде <tex>\begin{bmatrix} I_{N-K} & P\end{bmatrix}</tex>, где <tex>I_{N-K}</tex> — единичная матрица размера <tex>(N-K)\times(N-K)</tex>. Тогда в качестве порождающей матрицы, обеспечивающей систематическое кодирование, можно выбрать матрицу
+
Одним из способов построения порождающей матрицы кода по заданной проверочной матрице является преобразование проверочной матрицы к [http://ru.wikipedia.org/wiki/Ступенчатый_вид_по_строкам каноническому ступенчатому виду]. Такое преобразование всегда может быть сделано с помощью [http://en.wikipedia.org/wiki/Gaussian_elimination гауссовских исключений]. С точностью до перестановки столбцов канонический ступенчатый вид матрицы <tex>H</tex> эквивалентен её представлению в виде <tex>\begin{bmatrix} I_{m} & P\end{bmatrix}</tex>, где <tex>I_{m}</tex> — единичная матрица размера <tex>m{\times}m</tex>. Тогда в качестве порождающей матрицы, обеспечивающей систематическое кодирование, можно выбрать матрицу
-
:<tex>G = \begin{bmatrix}P\\ I_K\end{bmatrix}</tex>.
+
:<tex>G = \begin{bmatrix}P\\ I_k\end{bmatrix}</tex>.
-
Действительно, в этом случае <tex>HGt = (P+P)t = 0</tex>.
+
Действительно, в этом случае <tex>HGu = (P+P)u = 0</tex>.
=== Декодирование низкоплотностного кода ===
=== Декодирование низкоплотностного кода ===
Строка 35: Строка 33:
[[Image:GM13_task2_fg.png|thumb|Фактор-граф для (16,4)-низкоплотностного кода, в проверочной матрице которого в каждом столбце по 3 единицы, а в каждой строке - по 4 единицы.]]
[[Image:GM13_task2_fg.png|thumb|Фактор-граф для (16,4)-низкоплотностного кода, в проверочной матрице которого в каждом столбце по 3 единицы, а в каждой строке - по 4 единицы.]]
-
''Низкоплотностным кодом'' (или кодом с малой плотностью проверок на чётность) называется бинарный (N,K)-линейный блоковый код, в котором проверочная матрица <tex>H\in\{0,1\}^{(N-K)\times N}</tex> является сильно разреженной.
+
''Низкоплотностным кодом'' (или кодом с малой плотностью проверок на чётность) называется бинарный (n,k)-линейный блоковый код, в котором проверочная матрица <tex>H\in\{0,1\}^{m{\times}n}</tex> является сильно разреженной.
-
Рассмотрим бинарный симметричный канал для передачи данных. Здесь при передаче каждый бит независимо инвертируется с некоторой вероятностью q. В результате бинарный симметричный канал задает распределение <tex>p(y|x)</tex> для передаваемого кодового слова <tex>x\in\{0,1\}^N</tex> и полученного на выходе слова <tex>y\in\{0,1\}^N</tex> как
+
Рассмотрим бинарный симметричный канал для передачи данных. Здесь при передаче каждый бит независимо инвертируется с некоторой вероятностью <tex>q</tex>. В результате бинарный симметричный канал задает распределение <tex>p(w|v)</tex> для передаваемого кодового слова <tex>v\in\{0,1\}^n</tex> и полученного на выходе слова <tex>w\in\{0,1\}^n</tex> как
-
:<tex>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}</tex>.
+
:<tex>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}</tex>.
[http://en.wikipedia.org/wiki/Channel_capacity Пропускная способность] данного канала определяется величиной <tex>1+q\log_2q+(1-q)\log_2(1-q)</tex>.
[http://en.wikipedia.org/wiki/Channel_capacity Пропускная способность] данного канала определяется величиной <tex>1+q\log_2q+(1-q)\log_2(1-q)</tex>.
-
Объединяя низкоплотностный код с бинарным симметричным каналом, получаем следующую вероятностную модель для пары <tex>x,y</tex>:
+
Объединяя низкоплотностный код с бинарным симметричным каналом, получаем следующую вероятностную модель для пары <tex>w,v</tex>:
-
:<tex>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]</tex>.
+
:<tex>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]</tex>.
-
Здесь M=N-K — количество проверок на чётность, <tex>h_m^T</tex> — m-ая строка матрицы H, а <tex>I[\cdot]</tex> — индикаторная функция. Фактор-граф введённой модели показан на рис. справа.
+
Здесь <tex>m=n-k</tex> — количество проверок на чётность, <tex>h_j^T</tex> — j-ая строка матрицы H, а <tex>I[\cdot]</tex> — индикаторная функция. Фактор-граф введённой модели показан на рис. справа.
-
Процесс декодирования, т.е. восстановление кодового слова <tex>x</tex> по полученному слову <tex>y</tex>, предлагается осуществлять как
+
Назовём ''синдромом'' принятого слова <tex>w</tex> вектор <tex>s\in\{0,1\}^m</tex>, определяемый как <tex>s=Hw</tex>. Процесс передачи кодового слова по бинарному симметричному каналу можно представить как <tex>w=v+e</tex>, где <tex>e\in\{0,1\}^n</tex> — вектор ошибок (<tex>e_i=1</tex>, если в позиции <tex>i</tex> произошла ошибка). Тогда <tex>s = Hw = H(v+e) = Hv+He = He</tex>. Далее можно перейти от вероятностной модели для переменных <tex>v,w</tex> к аналогичной для переменных <tex>e,s</tex>:
-
<tex>x_n^* = \arg\max p(x_n|y)</tex>, а маргинальные распределения <tex>p(x_n|y)</tex> оценивать в помощью циклического алгоритма передачи сообщений (sum-product loopy BP) на фактор-графе. При этом для упрощения алгоритма декодирования предлагается избавиться от факторов-унарных потенциалов (путем их включения в сообщения от переменных <tex>x_n</tex> к факторам <tex>f_m</tex>), а в качестве расписания пересчёта сообщений выбрать параллельное расписание, при котором сначала все переменные одновременно посылают сообщения во все факторы, а затем все факторы одновременно посылают сообщения во все вершины.
+
:<tex>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]</tex>.
-
Введём обозначения <tex>N(n)=\{m:H_{mn}=1\}</tex> — множество факторов, в которых участвует переменная <tex>x_n</tex>, и <tex>N(m)=\{n:H_{mn}=1\}</tex> — множество переменных, которые входят в фактор <tex>f_m</tex>. Тогда общая схема алгоритма декодирования выглядит следующим образом:
+
Здесь <tex>p(e_i)=q^{e_i}(1-q)^{1-e_i}</tex>. Зная значение вектора ошибок <tex>\hat{e}</tex>, результат декодирования вычисляется как <tex>\hat{v}=w+\hat{e}</tex>. При использовании вероятностной модели для <tex>e,s</tex> тестирование алгоритма декодирования можно проводить без предварительной реализации процедуры кодирования.
 +
 
 +
При использовании побитовой функции потерь <tex>\lambda(e,\tilde{e}) = \sum_{i=1}^nI[e_i\neq \tilde{e}_i]</tex> оптимальная процедура декодирования связана с максимизацией маргиналов отдельных переменных, т.е. <tex>\hat{e}_n = \arg\max p(e_n|s)</tex>. Для поиска маргинальных распределений <tex>p(e_n|s)</tex> воспользуемся циклическим алгоритмом передачи сообщений (sum-product loopy BP) на фактор-графе. Введём обозначения <tex>N(i)=\{j:H_{ji}=1\}</tex> — множество факторов, в которых участвует переменная <tex>e_i</tex>, и <tex>N(j)=\{i:H_{ji}=1\}</tex> — множество переменных, которые входят в фактор <tex>h_j</tex>. Тогда общая схема алгоритма декодирования выглядит следующим образом:
1. Инициализация:
1. Инициализация:
-
: <tex>\mu_{x_n\rightarrow f_m}(x_n) = p(y_n|x_n)</tex>;
+
: <tex>\mu_{e_i\rightarrow h_j}(e_i) = p(e_i)</tex>;
2. Пересчет сообщений от факторов:
2. Пересчет сообщений от факторов:
-
: <tex>\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^'})</tex>;
+
: <tex>\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})</tex>;
-
3. Пересчет сообщений от переменных:
+
3. Пересчет сообщений от переменных и вычисление beliefs (оценок на маргинальные распределения):
-
: <tex>\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)</tex>;<br>
+
: <tex>\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)</tex>;<br>
-
: <tex>\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)</tex>;
+
: <tex>b_i(e_i)\propto p(e_i)\prod_{k\in N(i)}\mu_{h_{k}\rightarrow e_i}(e_i)</tex>;
-
: Символом <tex>\propto</tex> обозначается пропорциональность. Таким образом, при пересчете все сообщения от переменных и оценки на маргинальные распределения должны нормироваться так, чтобы <tex>\mu_{x_n\rightarrow f_m}(0)+\mu_{x_n\rightarrow f_m}(1)=1</tex> и <tex>\hat{p}_n(0|y)+\hat{p}_n(1|y)=1</tex>;
+
4. Оценка вектора ошибок:
-
4. Оценка кодового слова:
+
: <tex>\hat{e}_i=\arg\max b_i(e_i)</tex>;
-
: <tex>\hat{x}_n=\arg\max\hat{p}_n(x_n|y)</tex>;
+
5. Критерий остановки:
5. Критерий остановки:
-
: Если <tex>H\hat{x}=0</tex>, то выход алгоритма со статусом 0;
+
: Если <tex>H\hat{e}=s</tex>, то выход алгоритма со статусом 0;
-
: Если достигнуто максимальное число итераций или суммарная норма разности между сообщениями на текущей и предыдущей итерациях меньше определенного порога <tex>\varepsilon</tex>, то выход алгоритма со статусом -1;
+
: Если произошла стабилизация по всем <tex>b_i</tex>, то выход алгоритма со статусом 1;
 +
: Если достигнуто максимальное число итераций, то выход алгоритма со статусом 2;
: Переход к шагу 2.
: Переход к шагу 2.
-
При прямой реализации данного алгоритма на шаге 2 требуется рассмотрение для каждого фактора <tex>2^{N(m)-1}</tex> различных конфигураций переменных. Это может приводить как к низкой скорости пересчета сообщений, так и к большим требованиям по памяти.
+
=== Оптимизации в процедуре декодирования ===
-
Рассмотрим более эффективную схему реализации шага 2 путем его сведения к задаче вывода в графической модели с графом-цепочкой. Пусть нам необходимо вычислить сообщение <tex>\mu_{f_m\rightarrow x_n}(x_n)</tex>. Перенумеруем все переменные, входящие в m-ый фактор (кроме переменной <tex>x_n</tex>), как <tex>x_1,x_2,\dots,x_{N(m)-1}</tex>. Рассмотрим графическую модель с <tex>N(m)-1</tex> бинарными переменными <tex>s_i</tex> и графом
+
==== Пересчёт сообщений от факторов ====
-
[[Image:GM13_task2_chain.png|500px]].
+
При прямой реализации шага 2 общего алгоритма декодирования требуется рассмотрение для каждого фактора <tex>2^{|N(j)|-1}</tex> различных конфигураций значений переменных <tex>e_k</tex>. Это может приводить как к низкой скорости пересчета сообщений, так и к большим требованиям по памяти.
-
Положим, что <tex>s_i=\sum_{j=1}^ix_j</tex>. Определим априорное распределение и вероятность перехода как
+
Рассмотрим более эффективную схему реализации шага 2. Пусть нам необходимо вычислить сообщение <tex>\mu_{h_j\rightarrow e_i}(e_i)</tex>. Перенумеруем все переменные, входящие в j-ый фактор (кроме переменной <tex>e_i</tex>), как <tex>e_{(1)},e_{(2)},\dots,e_{(l)}</tex>, где <tex>l=|N(j)|-1</tex>. Тогда вычисление сообщения <tex>\mu_{h_j\rightarrow e_i}(e_i)</tex> можно записать как
-
:<tex>p(s_1) = \mu_{x_1\rightarrow f_m}(s_1)</tex>; <tex>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}</tex>
+
:<tex>\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)})</tex>.
-
Тогда требуемое сообщение <tex>\mu_{f_m\rightarrow x_n}(x_n)</tex> поэлементно равно маргинальному распределению <tex>p(s_{N(m)-1})</tex>, которое может быть эффективно вычислено путем однократной пересылки сообщений вдоль цепи от <tex>s_1</tex> до <tex>s_{N(m)-1}</tex>.
+
Данный результат можно интерпретировать как вычисление вероятности <tex>p(\sum_{k=1}^le_{(k)}=s_j+e_i)</tex> для набора независимых бинарных переменных <tex>e_{(k)}</tex> с распределением <tex>p_{(k)}(e_{(k)}) = \mu_{e_{(k)}\rightarrow h_j}(e_{(k)}).</tex> Обозначим через <tex>e^k</tex> сумму первых <tex>k</tex> переменных <tex>e_{(k)}</tex>, т.е. <tex>e^k=e_{(1)}+e_{(2)}+\dots+e_{(k)}</tex>. Тогда распределение на <tex>e^k</tex> можно итерационно пересчитывать по формуле
 +
:<tex>p_k(e^k) = \sum_{e^{k-1}}p_{k-1}(e^{k-1})p_{(k)}(e^k-e^{k-1})</tex>.
-
Дальнейшее повышение эффективности реализации шага 2 связано с рассмотрением разностей значений сообщений <tex>\delta\mu_{x_n\rightarrow f_m} = \mu_{x_n\rightarrow f_m}(0) - \mu_{x_n\rightarrow f_m}(1)</tex> и <tex>\delta\mu_{f_m\rightarrow x_n} = \mu_{f_m\rightarrow x_n}(0) - \mu_{f_m\rightarrow x_n}(1)</tex>. В терминах разностей можно показать, что новое значение <tex>\delta\mu_{f_m\rightarrow x_n}</tex> может быть вычислено как
+
В результате вычисление <tex>p_l(e^l)</tex> имеет линейную по <tex>l</tex> сложность.
-
:<tex>\delta\mu_{f_m\rightarrow x_n} = \prod_{n^'\in N(m)\backslash n}\delta\mu_{x_{n^'}\rightarrow f_m}</tex>.&nbsp;&nbsp;&nbsp;(*)
+
-
Зная значение <tex>\delta\mu_{f_m\rightarrow x_n}</tex> и учитывая условие нормировки на сообщения, сами сообщения от факторов могут быть вычислены как
+
-
:<tex>\mu_{f_m\rightarrow x_n}(0) = \frac{1}{2}(1+\delta\mu_{f_m\rightarrow x_n})</tex>; <tex>\mu_{f_m\rightarrow x_n}(1) = \frac{1}{2}(1-\delta\mu_{f_m\rightarrow x_n})</tex>.
+
-
Основное преимущество условия (*) по сравнению с пересчетом сообщений для <tex>s_i</tex> связано с тем, что формула (*) может быть реализована с помощью векторных операций в MATLAB.
+
Дальнейшее повышение эффективности реализации шага 2 связано с рассмотрением разностей вероятностей <tex>\delta p_{(k)} = p_{(k)}(0) - p_{(k)}(1)</tex> и <tex>\delta p_k = p_k(0)-p_k(1)</tex>. Нетрудно показать, что
 +
:<tex>\delta p_l = \prod_{k=1}^l\delta p_{(k)}</tex>. (*)
-
Покажем теперь справедливость условия (*). Для этого рассмотрим две произвольные бинарные независимые случайные величины <tex>u_1</tex> и <tex>u_2</tex>. Обозначим через <tex>p_i^u</tex> величину <tex>\mathbb{P}\{u_i=u\}</tex>. Тогда очевидно, что
+
Зная значение <tex>\delta p_l</tex> и учитывая условие нормировки <tex>p_l(0)+p_l(1)=1</tex>, сами вероятности могут быть вычислены как <tex>p_l(0)=\frac{1}{2}(1+\delta p_l)</tex>, <tex>p_l(1)=\frac{1}{2}(1-\delta p_l)</tex>. Таким образом, <tex>\mu_{h_j\rightarrow e_i}(e_i) = p_l(s_j+e_i)</tex>. Основное преимущество условия (*) по сравнению с последовательным пересчётом распределений <tex>p_k(e^k)</tex> связано с тем, что формула (*) может быть реализована с помощью векторных операций в MATLAB.
-
:<tex>\mathbb{P}\{u_1+u_2=0\} = p_1^0p_2^0 + p_1^1p_2^1</tex>; <tex>\mathbb{P}\{u_1+u_2=1\} = p_1^0p_2^1 + p_1^1p_2^0</tex>; <tex>\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</tex>.
+
 
-
Аналогичные рассуждения справедливы для произвольного количества случайных величин <tex>u</tex>. В частности,
+
==== Расписание пересчёта сообщений и дэмпфирование ====
-
:<tex>\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</tex>.
+
 
-
Отсюда немедленно следует справедливость (*).
+
Общая схема циклического алгоритма передачи сообщений оставляет определённый произвол в выборе расписания пересчёта сообщений. Обычно здесь рассматриваются следующие подходы:
 +
* ''Параллельное расписание''. В данном случае сначала все вершины посылают сообщения во все факторы, а затем все факторы посылают сообщения во все вершины.
 +
* ''Последовательное расписание''. Здесь выбирается некоторый (например, случайный) порядок последовательного пересчёта всех сообщений (и от вершин к факторам, и от факторов к вершинам). При этом данный порядок может меняться от итерации к итерации.
 +
 
 +
При использовании дэмпфирования с параметром <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 общего алгоритма декодирования.
== Формулировка задания ==
== Формулировка задания ==
# Реализовать алгоритм построения по заданной проверочной матрице чётности H порождающей матрицы кода G для систематического кодирования;
# Реализовать алгоритм построения по заданной проверочной матрице чётности H порождающей матрицы кода G для систематического кодирования;
-
# Реализовать алгоритм декодирования низкоплотностного кода на основе loopy BP; при реализации шага 2 пересчета сообщений от факторов к переменным необходимо использовать эффективные схемы, обозначенные выше; при реализации на MATLAB одной итерации схемы передачи сообщений использование вложенных циклов является нежелательным; провести временные замеры реализованного алгоритма для различных значений входных параметров;
+
# Реализовать алгоритм декодирования низкоплотностного кода на основе loopy BP; при реализации шага 2 пересчета сообщений от факторов к переменным необходимо использовать эффективные схемы, обозначенные выше; реализовать последовательное и параллельное расписание пересчёта сообщений, а также дэмпфирование сообщений; при реализации на MATLAB одной итерации схемы передачи сообщений использование вложенных циклов является нежелательным;
-
# Рассмотрим две характеристики качества кода — вероятность совершить ошибку хотя бы в одном бите при декодировании блока (<tex>p(\hat{x}\neq x)</tex>) и среднюю вероятность совершить ошибку при декодировании в одном бите (<tex>\frac{1}{N}\sum_{n=1}^Np(\hat{x}_n\neq x_n)</tex>). Требуется реализовать алгоритм оценки вероятности битовой и блоковой ошибки кода с помощью метода стат. испытаний (многократная случайная генерация слова <tex>t</tex>, его преобразование к кодовому слову <tex>x</tex>, передача по каналу с независимым инвертированием каждого бита с заданной вероятностью <tex>q</tex>, восстановление кодового слова <tex>\hat{x}</tex> с помощью алгоритма декодирования и подсчет необходимых характеристик);
+
# Провести эксперименты с различными расписаниями пересчёта сообщений и коэффициентами дэмпфирования; в частности, оценить долю стабилизировавшихся beliefs в зависимости от номера итерации алгоритма декодирования (усреднённую по различным запускам); оценить время работы алгоритма декодирования в зависимости от выбранного расписания и коэффициента дэмпфирования;
-
# Провести эксперименты по оцениванию битовой и блоковой ошибки низкоплотностного кода для различных значений длины кодового слова N, скорости кода R, вероятности инвертирования бита при передаче по каналу связи q и среднего количества единиц в столбце проверочной матрицы j. В частности, необходимо проанализировать следующие ситуации:
+
# Рассмотрим две характеристики качества кода — вероятность совершить ошибку хотя бы в одном бите при декодировании блока (<tex>p(\hat{e}\neq e)</tex>) и среднюю вероятность совершить ошибку при декодировании в одном бите (<tex>\frac{1}{n}\sum_{i=1}^np(\hat{e}_i\neq e_i)</tex>). Требуется реализовать алгоритм оценки вероятности битовой и блоковой ошибки кода с помощью метода стат. испытаний (многократная случайная генерация вектора ошибок <tex>e</tex> (<tex>e_i=1</tex> с вероятностью <tex>q</tex>), вычисление по нему синдрома <tex>s=He</tex>, восстановление вектора ошибок <tex>\hat{e}</tex> с помощью алгоритма декодирования и подсчет необходимых характеристик);
-
#* [http://ru.wikipedia.org/wiki/Теоремы_Шеннона_для_канала_с_шумами Теорема Шеннона] определяет пропускную способность канала как максимально допустимую скорость кода, при которой возможно осуществление надежной коммуникации. Требуется проверить, как меняются характеристики кода при изменении скорости R от минимального значения до пропускной способности канала.
+
# Провести эксперименты по оцениванию битовой и блоковой ошибки низкоплотностного кода для различных значений длины кодового слова <tex>n</tex>, скорости кода <tex>r</tex>, вероятности инвертирования бита при передаче по каналу связи <tex>q</tex> и среднего количества единиц в столбце проверочной матрицы <tex>j</tex>. В частности, необходимо проанализировать следующие ситуации:
-
#* Теорема Шеннона предполагает, что качество кода растет при увеличении длины кодового слова N. Требуется проверить это предположение.
+
#* [http://ru.wikipedia.org/wiki/Теоремы_Шеннона_для_канала_с_шумами Теорема Шеннона] определяет пропускную способность канала как максимально допустимую скорость кода, при которой возможно осуществление надежной коммуникации. Требуется проверить, как меняются характеристики кода при изменении скорости <tex>r</tex> от минимального значения до пропускной способности канала.
-
#* Одно из следствий теоремы Шеннона утверждает, что хорошими кодами являются коды со случайной проверочной матрицей H. В частности, здесь предполагается, что качество кода должно расти при увеличении среднего количества единиц в столбце проверочной матрицы j. Требуется проверить это утверждение для низкоплотностных кодов (путём рассмотрения нескольких значений j, начиная от 3).
+
#* Теорема Шеннона предполагает, что качество кода растет при увеличении длины кодового слова <tex>n</tex>. Требуется проверить это предположение.
-
# Составить отчёт в формате PDF с описанием всех проведённых исследований. Данный отчёт обязательно должен содержать описание отладочных тестов, которые проводились для проверки корректности реализованных алгоритмов. Также в этом отчёте должны быть графики поведения характеристик кода в зависимости от значений различных параметров.
+
#* Одно из следствий теоремы Шеннона утверждает, что хорошими кодами являются коды со случайной проверочной матрицей H. В частности, здесь предполагается, что качество кода должно расти при увеличении среднего количества единиц в столбце проверочной матрицы <tex>j</tex>. Требуется проверить это утверждение для низкоплотностных кодов (путём рассмотрения нескольких значений <tex>j</tex>, начиная от 3).
 +
# Провести эксперименты по сравнению низкоплотностного кода [[Практикум на ЭВМ (317)/2013/Коды БЧХ|с кодами БЧХ]]; в частности, провести эксперименты для различных <tex>q</tex> при фиксированной скорости кода <tex>r</tex>, а также для различных <tex>r</tex> при фиксированной <tex>q</tex>;
 +
# Составить отчёт в формате PDF с описанием всех проведённых исследований.
== Рекомендации по выполнению задания ==
== Рекомендации по выполнению задания ==
-
1. Разреженную проверочную матрицу кода заданных размеров можно строить с помощью случайной генерации (с соблюдением условия полноранговости). Однако, здесь рекомендуется воспользоваться [[Media:make_ldpc_mex.zip|готовой реализацией]], которая генерирует проверочную матрицу заданного размера с заданным количеством единиц в каждом столбце. При этом реализованный алгоритм старается сократить количество циклов длины 4 в генерируемой матрице, т.к. наличие коротких циклов в графе, как правило, значительно усложняет работу алгоритма loopy BP.
+
1. Разреженную проверочную матрицу кода заданных размеров можно строить с помощью случайной генерации (с соблюдением условия полноранговости). Однако, здесь рекомендуется воспользоваться [[Media:make_ldpc_mex.zip|готовой реализацией]], которая генерирует проверочную матрицу заданного размера с заданным количеством единиц в каждом столбце. При этом реализованный алгоритм старается сократить количество циклов длины 4 в генерируемой матрице, т.к. наличие коротких циклов в графе, как правило, усложняет работу алгоритма loopy BP.
-
2. При построении порождающей матрицы кода по заданной проверочной матрице с помощью алгоритма гауссовских исключений разрешается пользоваться сторонними кодами по реализации стандартных алгоритмов линейной алгебры (с соответствующими указаниями в отчёте). Однако, здесь рекомендуется реализовывать алгоритм гауссовских исключений в логике по модулю 2 самостоятельно, т.к. такая реализация может быть сделана значительно эффективнее по сравнению с общим случаем логики по модулю произвольного простого числа. Действительно, в логике по модулю два нет необходимости использовать операцию деления, а также искать ведущий элемент при очередном вычитании ввиду того, что все ненулевые элементы равны единице.
+
2. При тестировании алгоритма декодирования рекомендуется пробовать запускать итерационный процесс из различных начальных приближений для сообщений (не только из значений унарных потенциалов). При корректной реализации метода должна наблюдаться сходимость к одним и тем же финальным сообщениям независимо от начальных приближений.
-
 
+
-
3. Алгоритм декодирования удобно реализовывать в т.н. синдромном представлении. Полученное на выходе канала слово <tex>y</tex> можно представить как <tex>x+e</tex>, где <tex>x</tex> — посылаемое кодовое слово, а <tex>e\in\{0,1\}^N</tex> — вектор ошибок, которые вносит канал. Назовём величину <tex>z = Hy</tex> ''синдромом'' полученного сообщения. Тогда <tex>z=Hy=H(x+e)=He</tex>, и на этапе декодирования можно вместо <tex>x</tex> оценивать вектор ошибок <tex>e</tex> путем оценки аргмаксимумов маргинальных распределений <tex>p(e_n|z)</tex> в вероятностной модели
+
-
::<tex>p(e,z)\propto p(e)I[He=z]=\prod_{n=1}^Np(e_n)\prod_{m=1}^MI[h_m^Te=z_m]</tex>.
+
-
Для бинарного симметричного канала <tex>p(e_n) = q^{e_n}(1-q)^{e_n+1}</tex>. Зная <tex>e</tex>, искомое кодовое слово можно найти как <tex>x=y+e</tex>. Заметим, что реализация алгоритма декодирования в синдромном представлении избавляет от необходимости предварительной реализации алгоритма кодирования.
+
-
 
+
-
4. При тестировании алгоритма декодирования рекомендуется пробовать запускать итерационный процесс из различных начальных приближений для сообщений (не только из значений унарных потенциалов). При корректной реализации метода должна наблюдаться сходимость к одним и тем же финальным сообщениям независимо от начальных приближений.
+
== Оформление задания ==
== Оформление задания ==
-
Выполненное задание следует отправить письмом по адресу ''bayesml@gmail.com'' с заголовком письма «[ГМ13] Задание 2 <ФИО>». Убедительная просьба присылать выполненное задание '''только один раз''' с окончательным вариантом. Также убедительная просьба строго придерживаться заданных ниже прототипов реализуемых функций.
+
Выполненное задание следует отправить письмом по адресу ''bayesml@gmail.com'' с заголовком письма «[ГМ14] Задание 2 <ФИО>». Убедительная просьба присылать выполненное задание '''только один раз''' с окончательным вариантом. Также убедительная просьба строго придерживаться заданных ниже прототипов реализуемых функций.
Присланный вариант задания должен содержать в себе:
Присланный вариант задания должен содержать в себе:
Строка 133: Строка 136:
|
|
{|border="0"
{|border="0"
-
|H — проверочная матрица чётности, бинарная матрица размера MxN;
+
|H — проверочная матрица чётности, бинарная матрица размера m x n;
|}
|}
|-
|-
Строка 140: Строка 143:
|
|
{|
{|
-
|G — порождающая матрица кода, бинарная матрица размера Nx(N-M);
+
|G — порождающая матрица кода, бинарная матрица размера n x k;
|-
|-
|ind — номера позиций кодового слова, в которые копируются биты исходного сообщения, т.е. G(ind, :) является единичной матрицей.
|ind — номера позиций кодового слова, в которые копируются биты исходного сообщения, т.е. G(ind, :) является единичной матрицей.
Строка 151: Строка 154:
!''Алгоритм декодирования LDPC-кода в синдромном представлении''
!''Алгоритм декодирования LDPC-кода в синдромном представлении''
|-
|-
-
|[e, status] = '''ldpc_decoding'''(z, H, q, param_name1, param_value1, ...)
+
|[e, status] = '''ldpc_decoding'''(s, H, q, param_name1, param_value1, ...)
|-
|-
|ВХОД
|ВХОД
Строка 157: Строка 160:
|
|
{|border="0"
{|border="0"
-
|z — наблюдаемый синдром, бинарный вектор-столбец длины M;
+
|s — наблюдаемый синдром, бинарный вектор-столбец длины m;
|-
|-
-
|H — проверочная матрица чётности, бинарная матрица размера MxN;
+
|H — проверочная матрица чётности, бинарная матрица размера m x n;
|-
|-
|q — вероятность инверсии бита при передаче по каналу связи, число от 0 до 0.5;
|q — вероятность инверсии бита при передаче по каналу связи, число от 0 до 0.5;
Строка 167: Строка 170:
|
|
{|border="0"
{|border="0"
 +
|'schedule' — расписание пересчёта сообщений, возможные значения 'parallel' и 'sequential', по умолчанию = 'parallel';
 +
|-
 +
|'damping' — коэффициент дэмпфирования при пересчёте сообщений, число от 0 до 1, по умолчанию = 1;
 +
|-
|'max_iter' — максимальное число итераций алгоритма декодирования, число, по умолчанию = 200;
|'max_iter' — максимальное число итераций алгоритма декодирования, число, по умолчанию = 200;
|-
|-
-
|'eps' — порог стабилизации для сообщений, число, по умолчанию = 1e-4;
+
|'eps' — порог стабилизации для beliefs, число, по умолчанию = 1e-4;
|-
|-
|'display' — режим отображения, true или false, если true, то отображается промежуточная информация на итерациях, например, номер итерации, текущее число ошибок декодирования, невязка для сообщений и т.д.
|'display' — режим отображения, true или false, если true, то отображается промежуточная информация на итерациях, например, номер итерации, текущее число ошибок декодирования, невязка для сообщений и т.д.
Строка 180: Строка 187:
|
|
{|
{|
-
|e — восстановленный вектор ошибок, бинарный вектор-столбец длины N;
+
|e — восстановленный вектор ошибок, бинарный вектор-столбец длины n;
|-
|-
-
|status — результат декодирования, равен 0, если вектор e восстановлен без ошибок, равен -1, если произошел выход по максимальному числу итераций или стабилизации значений сообщений.
+
|status — результат декодирования, равен 0, если найден вектор e, соответствующий входному синдрому s, равен 1, если произошла стабилизация значений beliefs, равен 2, если произошел выход по максимальному числу итераций.
|}
|}
|}
|}
Строка 191: Строка 198:
!''Оценка характеристик LDPC-кода с помощью метода Монте Карло''
!''Оценка характеристик LDPC-кода с помощью метода Монте Карло''
|-
|-
-
|[err_bit, err_block, diver] = '''ldpc_mc'''(H, G, q, num_points)
+
|[err_bit, err_block, diver] = '''ldpc_mc'''(H, q, num_points)
|-
|-
|ВХОД
|ВХОД
Строка 197: Строка 204:
|
|
{|border="0"
{|border="0"
-
|H — проверочная матрица чётности, бинарная матрица размера MxN;
+
|H — проверочная матрица чётности, бинарная матрица размера m x n;
-
|-
+
-
|G — порождающая матрица кода, бинарная матрица размера Nx(N-M);
+
|-
|-
|q — вероятность инверсии бита при передаче по каналу связи, число от 0 до 0.5;
|q — вероятность инверсии бита при передаче по каналу связи, число от 0 до 0.5;
Строка 210: Строка 215:
|
|
{|
{|
-
|err_bit — вероятность битовой ошибки декодирования (относительно N бит кодового слова), число от 0 до 1;
+
|err_bit — вероятность битовой ошибки декодирования (относительно n бит кодового слова), число от 0 до 1, вычисляется по тем ситуациям, когда алгоритм декодирования сошёлся (status < 2);
|-
|-
-
|err_block — вероятность блоковой ошибки декодирования, число от 0 до 1;
+
|err_block — вероятность блоковой ошибки декодирования, число от 0 до 1, вычисляется по тем ситуациям, когда алгоритм декодирования сошёлся (status < 2)
|-
|-
|diver — доля ситуаций расходимости алгоритма декодирования, число от 0 до 1.
|diver — доля ситуаций расходимости алгоритма декодирования, число от 0 до 1.
|}
|}
|}
|}
 +
 +
[[Категория:Учебные задания]]
 +
[[Категория:Байесовские методы]]

Текущая версия

Содержание

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

Начало выполнения задания: 5 марта 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 общего алгоритма декодирования требуется рассмотрение для каждого фактора 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. Провести эксперименты с различными расписаниями пересчёта сообщений и коэффициентами дэмпфирования; в частности, оценить долю стабилизировавшихся beliefs в зависимости от номера итерации алгоритма декодирования (усреднённую по различным запускам); оценить время работы алгоритма декодирования в зависимости от выбранного расписания и коэффициента дэмпфирования;
  4. Рассмотрим две характеристики качества кода — вероятность совершить ошибку хотя бы в одном бите при декодировании блока (p(\hat{e}\neq e)) и среднюю вероятность совершить ошибку при декодировании в одном бите (\frac{1}{n}\sum_{i=1}^np(\hat{e}_i\neq e_i)). Требуется реализовать алгоритм оценки вероятности битовой и блоковой ошибки кода с помощью метода стат. испытаний (многократная случайная генерация вектора ошибок e (e_i=1 с вероятностью q), вычисление по нему синдрома s=He, восстановление вектора ошибок \hat{e} с помощью алгоритма декодирования и подсчет необходимых характеристик);
  5. Провести эксперименты по оцениванию битовой и блоковой ошибки низкоплотностного кода для различных значений длины кодового слова n, скорости кода r, вероятности инвертирования бита при передаче по каналу связи q и среднего количества единиц в столбце проверочной матрицы j. В частности, необходимо проанализировать следующие ситуации:
    • Теорема Шеннона определяет пропускную способность канала как максимально допустимую скорость кода, при которой возможно осуществление надежной коммуникации. Требуется проверить, как меняются характеристики кода при изменении скорости r от минимального значения до пропускной способности канала.
    • Теорема Шеннона предполагает, что качество кода растет при увеличении длины кодового слова n. Требуется проверить это предположение.
    • Одно из следствий теоремы Шеннона утверждает, что хорошими кодами являются коды со случайной проверочной матрицей H. В частности, здесь предполагается, что качество кода должно расти при увеличении среднего количества единиц в столбце проверочной матрицы j. Требуется проверить это утверждение для низкоплотностных кодов (путём рассмотрения нескольких значений j, начиная от 3).
  6. Провести эксперименты по сравнению низкоплотностного кода с кодами БЧХ; в частности, провести эксперименты для различных q при фиксированной скорости кода r, а также для различных r при фиксированной q;
  7. Составить отчёт в формате PDF с описанием всех проведённых исследований.

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

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

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

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

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

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

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

 

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

 

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

 

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