Wi-Fiでストリーミング動画を視聴しているとき、電波環境が多少悪くても映像が途切れにくいのはなぜでしょうか。SSD（ソリッドステートドライブ）に大量のデータを書き込んでも、読み出し時に1ビットも化けないのはなぜでしょうか。その裏側で活躍しているのがLDPC符号（Low-Density Parity-Check Code）です。
LDPC符号は、パリティ検査行列がほとんどゼロで埋め尽くされた「疎行列」であることが特徴の誤り訂正符号です。この「疎」であるという性質が、信念伝搬（Belief Propagation）という効率的な復号アルゴリズムを可能にし、理論限界であるシャノン限界に極めて近い性能を実現します。
LDPC符号を理解すると、以下のような広い分野への見通しが開けます。
- 5G NR（New Radio）: 3GPPの5Gでは、データチャネル（PDSCH/PUSCH）の誤り訂正符号としてLDPC符号が採用されています。高速・大容量の5G通信を支える基盤技術です
- Wi-Fi（IEEE 802.11n/ac/ax/be）: Wi-Fi 4以降の全ての規格でLDPC符号がオプションまたは必須として採用されており、無線LANの信頼性向上に貢献しています
- ストレージ（SSD/HDD）: NANDフラッシュメモリの微細化に伴う信頼性低下を補うため、最新のSSDコントローラにはLDPC符号が搭載されています
- 確率的グラフィカルモデル: LDPC符号のBP復号は、ベイズネットワークやマルコフ確率場での推論アルゴリズムと密接に関連しており、機械学習や統計物理学との橋渡しとなります
本記事の内容
- LDPC符号の歴史 — Gallagerの先見性と再発見の経緯
- パリティ検査行列 $\bm{H}$ の構造と正則・非正則LDPC符号
- タナーグラフ（二部グラフ表現）の直感的理解
- 信念伝搬（BP）復号アルゴリズム：対数尤度比メッセージの更新則
- 密度進化による復号閾値解析
- Pythonでの正則LDPC符号化・BP復号シミュレーションとBERカーブの描画
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
LDPC符号の歴史
Gallagerの先見的な研究（1960-1963）
LDPC符号の歴史は、Robert G. Gallagerの1960年のMIT博士論文に遡ります。Gallagerは、パリティ検査行列の各行と各列に含まれる1の数が少ない（低密度な）符号を提案し、確率的な繰り返し復号アルゴリズムの性能を理論的に解析しました。
1962年には、この成果をIEEE Transactions on Information Theoryに発表し、1963年にはMIT Pressからモノグラフ “Low-Density Parity-Check Codes” を出版しています。Gallagerは、この符号が大きな符号長に対してシャノン限界に近い性能を発揮することを示しました。
しかし、1960年代のコンピュータ技術では、大きな符号長のLDPC符号を実装することは事実上不可能でした。また、同時期にReed-Solomon符号や畳み込み符号+ビタビ復号といった、当時の計算能力で実装可能な符号が次々と開発されたこともあり、LDPC符号は約30年間ほとんど忘れ去られることになります。
再発見と復興（1990年代）
1993年にターボ符号がシャノン限界に0.5 dB以内の性能を示したことは、符号理論コミュニティに衝撃を与えました。この衝撃がきっかけとなり、研究者たちはシャノン限界に近い性能を持つ他の符号の可能性を再探索し始めます。
1996年、David MacKayとRadford Nealは、Gallagerの忘れられていたLDPC符号を再発見し、現代の計算機でシミュレーションを行いました。その結果、適切に設計されたLDPC符号がターボ符号に匹敵する（場合によってはそれを上回る）性能を示すことが確認されたのです。
MacKayとNealの研究は、LDPC符号の再評価を引き起こしました。その後、Richardson、Shokrollahi、Urbankeらによる密度進化（Density Evolution）の理論、Luby、Mitzenmacher、Shokrollahi、Spielmanらによる非正則LDPC符号の設計、Richardsonらによるシャノン限界まで0.0045 dBという驚異的な性能を持つLDPC符号の設計など、急速な理論的発展が続きました。
現在の位置づけ
LDPC符号は現在、最も広く実用化されている近シャノン限界符号の一つです。ターボ符号と比較した場合のLDPC符号の主な利点は以下の通りです。
- 復号器の並列化: タナーグラフ上のメッセージ更新は局所的であり、高度に並列化できる。これにより、高スループットの復号器が実現できる
- エラーフロア: 適切に設計されたLDPC符号は、ターボ符号よりも低いエラーフロアを達成できる
- 柔軟な符号設計: 次数分布の最適化により、特定の通信路に合わせた符号設計が可能
これらの利点から、5G NRではデータチャネルにLDPC符号が選択されました。
それでは、LDPC符号の数学的な構造を見ていきましょう。まず、全ての線形ブロック符号の基盤であるパリティ検査行列から始めます。
パリティ検査行列 $\bm{H}$ の構造
線形ブロック符号の復習
LDPC符号の構造を理解するために、まず線形ブロック符号の基本を確認しておきましょう。
図書館の本に例えてみます。各本にはISBNコード（10桁または13桁）がありますが、最後の1桁はチェックディジットです。他の桁から計算で求まるこの冗長な桁があることで、番号の入力ミスを検出できます。LDPC符号はこれと同じ発想の延長線上にあります — ただし、1つのチェック式ではなく、数百から数千のチェック式を同時に使い、しかもそれぞれのチェック式が少数のビットだけを参照するというのがポイントです。
$(N, K)$ 線形ブロック符号は、$K$ ビットの情報語を $N$ ビットの符号語に変換する符号です（$N > K$）。符号化率は $R = K/N$ です。
符号語 $\bm{c} = (c_1, c_2, \dots, c_N)$ は、$M \times N$ のパリティ検査行列 $\bm{H}$（$M = N – K$）に対して次の条件を満たします。
$$ \bm{H} \bm{c}^T = \bm{0} $$
ここで演算は $\text{GF}(2)$（二元体、すなわち排他的論理和）上で行います。つまり、$\bm{H}$ の各行は1つのパリティ検査方程式を表し、全ての検査方程式を同時に満たすベクトル $\bm{c}$ が有効な符号語です。
「低密度」の意味
LDPC符号のパリティ検査行列 $\bm{H}$ は、ほとんどの要素が0で、わずかな要素だけが1である疎行列（sparse matrix）です。これが “Low-Density” の意味です。
一般的な線形ブロック符号のパリティ検査行列では、1の密度は行列サイズに対して一定の割合を占めます。例えば、BCH符号やReed-Solomon符号の検査行列は密行列です。密行列に対する復号は、符号長が大きくなると計算量が急激に増大し、実用的ではなくなります。
一方、LDPC符号では $\bm{H}$ の各行に含まれる1の数（行重み）と各列に含まれる1の数（列重み）が、符号長 $N$ に対して一定（$N$ に依存しない固定値）です。このため、行列全体に含まれる1の数は $O(N)$ であり、$N \times (N-K)$ の行列サイズに対する1の密度は $N$ が大きくなるにつれて0に近づきます。
正則LDPC符号と非正則LDPC符号
LDPC符号は、パリティ検査行列の1の分布によって2つのカテゴリに分類されます。
正則（regular）LDPC符号: $\bm{H}$ の全ての列の重み（1の数）が $d_v$ で均一、全ての行の重みが $d_c$ で均一な符号です。$(d_v, d_c)$-正則LDPC符号と呼ばれます。$M$ 行 $N$ 列の行列なので、1の総数は $N d_v = M d_c$ を満たします。符号化率は次のように表されます。
$$ R = 1 – \frac{M}{N} = 1 – \frac{d_v}{d_c} $$
例えば、$(3, 6)$-正則LDPC符号は $R = 1 – 3/6 = 1/2$ の符号化率を持ちます。Gallagerが元々提案したのはこの正則LDPC符号です。
非正則（irregular）LDPC符号: 列重みや行重みが均一ではなく、確率分布に従う符号です。次数分布多項式で特徴づけられます。
$$ \lambda(x) = \sum_{i=2}^{d_{v,\max}} \lambda_i x^{i-1}, \quad \rho(x) = \sum_{j=2}^{d_{c,\max}} \rho_j x^{j-1} $$
ここで、$\lambda_i$ は辺の観点から見て列重み $i$ の割合、$\rho_j$ は辺の観点から見て行重み $j$ の割合です。$\lambda(x)$ を変数ノード次数分布、$\rho(x)$ をチェックノード次数分布と呼びます。
非正則LDPC符号は正則LDPC符号よりも優れた性能を持つことが知られています。密度進化を用いた次数分布の最適化により、Richardson、Shokrollahi、Urbankeらはシャノン限界まで0.0045 dBという性能を持つLDPC符号を設計しました。5G NRで採用されているLDPC符号も非正則構造です。
パリティ検査行列の構造がわかったところで、次にこの行列をグラフとして表現する方法を見ていきましょう。グラフ表現は、BP復号アルゴリズムの理解と設計に不可欠です。
タナーグラフ
二部グラフとしての表現
タナーグラフ（Tanner Graph）は、1981年にR. Michael Tannerが提案した、LDPC符号のパリティ検査行列を二部グラフとして表現する方法です。
二部グラフとは、頂点を2つのグループに分け、同じグループ内の頂点同士には辺がなく、異なるグループの頂点間にのみ辺があるグラフです。タナーグラフでは、次の2種類のノードを持ちます。
変数ノード（Variable Node, VN）: 符号語の各ビット $c_n$（$n = 1, 2, \dots, N$）に対応する $N$ 個のノードです。図では通常、丸（$\bigcirc$）で描かれます。
チェックノード（Check Node, CN）: パリティ検査方程式の各行 $m$（$m = 1, 2, \dots, M$）に対応する $M$ 個のノードです。図では通常、四角（$\square$）で描かれます。
パリティ検査行列 $\bm{H}$ の $(m, n)$ 要素が1であるとき、チェックノード $m$ と変数ノード $n$ の間に辺を引きます。つまり、$\bm{H}$ の1の位置がタナーグラフの辺に対応します。
具体例
例えば、次のパリティ検査行列を考えます。
$$ \bm{H} = \begin{pmatrix} 1 & 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 & 0 & 1 \end{pmatrix} $$
これは $(N, K) = (6, 3)$ の符号で、3つのパリティ検査方程式があります。
- チェック1: $c_1 \oplus c_2 \oplus c_4 = 0$
- チェック2: $c_2 \oplus c_3 \oplus c_5 = 0$
- チェック3: $c_1 \oplus c_3 \oplus c_6 = 0$
タナーグラフでは、6個の変数ノード（$v_1$ 〜 $v_6$）と3個のチェックノード（$f_1$ 〜 $f_3$）があり、$\bm{H}$ の1の位置に対応する辺で接続されます。各変数ノードの次数（列重み）は2、各チェックノードの次数（行重み）は3なので、これは $(2, 3)$-正則LDPC符号です。
サイクルとガースの重要性
タナーグラフにおいて、BP復号の性能に大きく影響するのがサイクル（cycle）の存在です。サイクルとは、あるノードから出発して辺を辿り、同じノードに戻ってくる経路のことです。
タナーグラフは二部グラフなので、サイクルの長さは必ず偶数になります。最短のサイクルの長さをガース（girth）と呼びます。二部グラフの最短サイクルは長さ4です。
なぜサイクルが問題なのでしょうか。BP復号は、メッセージがグラフ上を伝搬することで動作しますが、短いサイクルがあると、あるノードから送り出されたメッセージがすぐに自分自身に戻ってきてしまいます。これは、独立であるべき情報が循環的に再利用されることを意味し、BP復号の前提（メッセージの独立性）が崩れます。
したがって、LDPC符号の設計ではガースをできるだけ大きく（最低でも6以上）することが重要です。ガースが大きいほど、BP復号はより正確に動作します。ただし、ガースを大きくするのと符号の距離特性を良くするのは必ずしも両立しないため、実際の設計ではトレードオフが存在します。
タナーグラフの構造が理解できたところで、いよいよこのグラフ上でメッセージを伝搬させるBP復号アルゴリズムの詳細に入ります。
信念伝搬（BP）復号アルゴリズム
メッセージパッシングの直感
BP復号を直感的に理解するために、次のような状況を想像してみてください。
6人のメンバーがいるチームで、3つの「ルール」があるとします。
- ルール1: メンバーA、B、Dの「合計」は偶数でなければならない
- ルール2: メンバーB、C、Eの「合計」は偶数でなければならない
- ルール3: メンバーA、C、Fの「合計」は偶数でなければならない
各メンバーは自分の値（0か1）についてある程度の確信を持っていますが、雑音のせいで完全には確定していません。ここで、各メンバーがルールを通じて互いに情報を交換すると、確信度が高まっていきます。
例えば、メンバーBがルール1を確認すると、AとDの情報から自分の値に関する追加の手がかりが得られます。同時に、ルール2からもCとEの情報を通じて手がかりが得られます。複数のルールから独立な情報が集まることで、Bは自分の値にますます自信を持てるようになります。
BP復号はまさにこのプロセスを数学的に定式化したものです。
対数尤度比（LLR）メッセージの定義
BP復号では、タナーグラフの辺に沿ってLLRメッセージを伝搬させます。以下の2種類のメッセージを定義します。
変数ノード $n$ からチェックノード $m$ へのメッセージ $L_{n \to m}$:
変数ノード $n$ が、チェックノード $m$ を除く全ての隣接チェックノードから集めた情報を基に、ビット $c_n$ の推定値を $L_{n \to m}$ として送ります。
チェックノード $m$ から変数ノード $n$ へのメッセージ $L_{m \to n}$:
チェックノード $m$ が、変数ノード $n$ を除く全ての隣接変数ノードから集めた情報を基に、パリティ検査方程式を使って $c_n$ に関する情報を $L_{m \to n}$ として送ります。
これらのメッセージは全てLLR形式です。正の値は $c_n = 0$ を示唆し、負の値は $c_n = 1$ を示唆します。絶対値が大きいほど確信度が高いことを意味します。
初期化
AWGN通信路でBPSK変調（$c_n = 0 \to x_n = +1$、$c_n = 1 \to x_n = -1$）を考えます。受信信号 $y_n = x_n + w_n$（$w_n \sim \mathcal{N}(0, \sigma^2)$）に対するチャネルLLRは次の通りです。
$$ L_{\text{ch}}(n) = \frac{2}{\sigma^2} y_n $$
初期化として、変数ノードからチェックノードへの最初のメッセージをチャネルLLRに設定します。
$$ L_{n \to m}^{(0)} = L_{\text{ch}}(n) $$
チェックノード更新（Check Node Update）
チェックノード $m$ は、パリティ検査方程式「接続された変数ノードの値のXOR合計が0」を使って情報を処理します。チェックノード $m$ に接続された変数ノードの集合を $\mathcal{N}(m)$ とします。
チェックノード $m$ から変数ノード $n$ へのメッセージ更新則は次のように書けます。
$$ L_{m \to n} = 2 \tanh^{-1} \left( \prod_{n’ \in \mathcal{N}(m) \setminus n} \tanh\left(\frac{L_{n’ \to m}}{2}\right) \right) $$
この式の意味を直感的に説明します。$\tanh(L/2)$ は確率領域での変換であり、$L > 0$ のとき正、$L < 0$ のとき負の値をとります。チェックノードはXOR演算（偶数パリティ）を表すので、接続された全ての変数ノード（$n$ を除く）の「符号の積」が、$c_n$ に対する情報を与えます。
もう少し具体的に言うと、パリティ検査方程式 $c_{n_1} \oplus c_{n_2} \oplus \dots \oplus c_{n_{d_c}} = 0$ から $c_n$ の値を推定するには、残りの $d_c – 1$ 個のビットの値がわかればよいのです。残りのビットの推定値に基づいて、パリティ検査を満たすような $c_n$ の値を推測する — これがチェックノード更新の本質です。
チェックノード更新には重要な性質があります。入力メッセージの絶対値が小さい（確信度が低い）ものがあると、出力メッセージの絶対値も小さくなります。つまり、最も不確実なビットが全体の品質を制限するのです。これは直感に合います — 1つでも確信のないビットがあると、パリティ検査からの情報は限られます。
変数ノード更新（Variable Node Update）
変数ノード $n$ に接続されたチェックノードの集合を $\mathcal{M}(n)$ とします。変数ノード $n$ からチェックノード $m$ へのメッセージ更新則は次の通りです。
$$ L_{n \to m} = L_{\text{ch}}(n) + \sum_{m’ \in \mathcal{M}(n) \setminus m} L_{m’ \to n} $$
変数ノード更新はチェックノード更新よりもシンプルです。チャネルからの直接情報（$L_{\text{ch}}(n)$）と、他の全てのチェックノード（$m$ を除く）から受け取った情報の合計です。
メッセージの宛先であるチェックノード $m$ からの情報を除いている点が重要です。これは、$m$ に送るメッセージに $m$ 自身からの情報を含めると、情報が循環してしまい（エコー効果）、BP復号の前提が崩れるためです。
事後LLRと硬判定
各反復の最後に、変数ノード $n$ の事後LLRを全ての受信情報の合計として計算します。
$$ L_{\text{total}}(n) = L_{\text{ch}}(n) + \sum_{m \in \mathcal{M}(n)} L_{m \to n} $$
ここでは、全てのチェックノードからのメッセージを含めます（どこかに送るメッセージではなく、最終推定なので）。硬判定は次の通りです。
$$ \hat{c}_n = \begin{cases} 0 & \text{if } L_{\text{total}}(n) \geq 0 \\ 1 & \text{if } L_{\text{total}}(n) < 0 \end{cases} $$
停止条件
反復は以下のいずれかの条件で停止します。
- パリティ検査の充足: $\bm{H} \hat{\bm{c}}^T = \bm{0}$ が成立すれば、有効な符号語が見つかったので停止
- 最大反復回数: 所定の反復回数（通常20〜50回）に達したら停止
パリティ検査による早期停止は、低SNR領域（多くの反復が必要）では計算量を節約できませんが、高SNR領域では数回の反復で停止できるため、平均的な復号遅延を大幅に削減できます。
Min-Sum近似
チェックノード更新における $\tanh$ と $\tanh^{-1}$ の計算は計算コストが高いため、実装ではしばしばMin-Sum近似が使われます。
$$ L_{m \to n} \approx \left(\prod_{n’ \in \mathcal{N}(m) \setminus n} \text{sign}(L_{n’ \to m})\right) \cdot \min_{n’ \in \mathcal{N}(m) \setminus n} |L_{n’ \to m}| $$
この近似は、$\tanh$ 演算を符号（sign）の積と絶対値の最小値で置き換えたものです。直感的には、「パリティ検査の結果の符号は入力符号の積で決まり、確信度は最も不確実な入力に律速される」という性質を利用しています。
Min-Sum近似は0.2〜0.5 dBの性能劣化がありますが、補正係数 $\alpha$（通常0.75〜0.9）を絶対値に乗じるNormalized Min-Sumや、固定のオフセット $\beta$ を絶対値から引くOffset Min-Sumにより、性能劣化を大幅に軽減できます。5G NRの復号器では、これらの近似アルゴリズムが広く使われています。
ここまでで、BP復号アルゴリズムの具体的な動作を理解しました。次のセクションでは、このアルゴリズムの理論的な性能限界を解析するための強力な手法 — 密度進化 — を紹介します。
密度進化による閾値解析
密度進化の概念
BP復号の性能は、符号長 $N$ とSNRに依存しますが、符号長 $N \to \infty$ の極限では、その振る舞いが鋭い閾値現象を示します。つまり、あるSNR（または雑音分散 $\sigma^2$）の閾値 $\sigma^*$ が存在し、$\sigma < \sigma^*$ ではBP復号のビット誤り率が0に収束し、$\sigma > \sigma^*$ ではBERが正の値に留まります。
この閾値を計算する手法が密度進化（Density Evolution, DE）であり、Richardson、Shokrollahi、Urbankeらによって2001年に提出されました。
密度進化の基本的なアイデアは次の通りです。符号長 $N \to \infty$ の極限では、タナーグラフは局所的に木（サイクルなしのグラフ）に見えます。木構造であれば、BP復号の各反復でのメッセージは互いに独立であり、その確率分布を正確に追跡できます。
密度進化の手順
密度進化は、各反復におけるLLRメッセージの確率密度関数の進化を追跡します。
初期化: チャネルLLR $L_{\text{ch}}$ の分布はチャネルモデルで決まります。AWGN通信路・BPSKの場合、$L_{\text{ch}} \sim \mathcal{N}(2/\sigma^2, 4/\sigma^2)$、すなわち平均 $\mu = 2/\sigma^2$、分散 $\sigma_L^2 = 4/\sigma^2$ のガウス分布です。
反復ステップ: 各反復 $\ell$ において、以下を計算します。
- 変数ノード → チェックノードメッセージの分布: 変数ノード更新は加算なので、入力メッセージの分布を畳み込みます
- チェックノード → 変数ノードメッセージの分布: チェックノード更新は $\tanh$ 変換を含むため、分布の変換が必要です
正確な密度進化は分布全体を追跡するため計算量が大きくなります。実用的には、ガウス近似（Gaussian Approximation）を用いて、メッセージの分布がガウス分布（対称条件：平均 = 分散/2）であると仮定することで、1つのパラメータ（平均値）のみを追跡する簡略化が可能です。
ガウス近似による閾値計算
ガウス近似では、変数ノードからチェックノードへのメッセージの平均を $\mu_v^{(\ell)}$ で追跡します。
変数ノード更新: 次数 $d_v$ の変数ノードからのメッセージの平均は、チャネルLLRの平均 $\mu_0 = 2/\sigma^2$ とチェックノードからのメッセージの平均 $\mu_c^{(\ell)}$ を使って次のように表されます。
$$ \mu_v^{(\ell+1)} = \mu_0 + (d_v – 1) \mu_c^{(\ell)} $$
チェックノード更新: 次数 $d_c$ のチェックノードからのメッセージの平均は、近似的に次の関数 $\phi$ を用いて計算されます。
$$ \mu_c^{(\ell)} = \phi^{-1}\left(1 – \left[1 – \phi(\mu_v^{(\ell)})\right]^{d_c – 1}\right) $$
ここで $\phi(\mu) = 1 – \frac{1}{\sqrt{4\pi\mu}} \int_{-\infty}^{\infty} \tanh\left(\frac{x}{2}\right) e^{-\frac{(x-\mu)^2}{4\mu}} dx$ です。この関数は数値的に事前計算しテーブル化できます。
閾値の決定: 初期値 $\mu_v^{(0)} = \mu_0 = 2/\sigma^2$ から反復を開始し、$\mu_v^{(\ell)} \to \infty$ に発散する（BER → 0）か、有限値に留まる（BER > 0）かを調べます。発散する最大の $\sigma^2$（最小のSNR）が閾値 $\sigma^*$ です。
例えば、$(3, 6)$-正則LDPC符号の閾値は $\sigma^* \approx 0.881$（$E_b/N_0 \approx 1.11$ dB）であることが知られています。AWGN通信路で $R = 1/2$ のシャノン限界は $E_b/N_0 \approx 0.19$ dBなので、正則LDPC符号はシャノン限界から約0.9 dBの位置にあります。
非正則LDPC符号の次数分布を最適化することで、この差をさらに縮めることができます。Richardson、Shokrollahi、Urbankeらは、最適化された非正則LDPC符号でシャノン限界まで0.0045 dBという驚異的な結果を報告しています。
密度進化の理論が示す閾値が、実際のシミュレーションでどの程度正確に現れるかを確認してみましょう。次のセクションでは、正則LDPC符号のPython実装とBERシミュレーションを行います。
Pythonでの実装
パリティ検査行列の生成
まず、$(d_v, d_c)$-正則LDPC符号のパリティ検査行列を生成する関数を実装します。Gallagerの方法に基づく構成を使います。
import numpy as np
def make_regular_ldpc(N, dv, dc):
"""
(dv, dc)-正則LDPC符号のパリティ検査行列を生成（Gallager法）
N: 符号長（dc の倍数であること）
dv: 列重み（変数ノード次数）
dc: 行重み（チェックノード次数）
戻り値: パリティ検査行列 H (M x N, M = N*dv/dc)
"""
M = N * dv // dc # パリティ検査数
# 基本サブ行列: 各列に1が1つだけある M x N 行列
# 最初のサブ行列: 1が連続ブロックに配置
H_sub0 = np.zeros((M, N), dtype=int)
for j in range(N):
row = j * M // N
H_sub0[row, j] = 1
# dv 個のサブ行列を重ねる（2番目以降は列をランダム置換）
H = H_sub0.copy()
for _ in range(dv - 1):
perm = np.random.permutation(N)
H += H_sub0[:, perm]
# 値を 0/1 に変換（重複を無視して二元行列にする）
H = (H % 2).astype(int)
return H
このコードでは、Gallagerの構成法に基づいてパリティ検査行列を生成しています。まず「各列にちょうど1つの1がある」基本サブ行列を作り、列の置換を変えながら $d_v$ 枚のサブ行列を重ね合わせます。これにより、各列の重みが $d_v$、各行の重みが約 $d_c$ の正則構造が得られます。重複（同じ位置に複数回1が加わる）はmod 2演算で処理しています。
組織符号化
パリティ検査行列 $\bm{H}$ から組織形式の生成行列 $\bm{G}$ を構成して符号化します。
def ldpc_encode(u, H):
"""
LDPC符号化（ガウス消去法で組織形式を構成）
u: 情報ビット列
H: パリティ検査行列
戻り値: 符号語 c
"""
M, N = H.shape
K = N - M
# GF(2) 上のガウス消去法でHを [P | I_M] 形式に変換
H_work = H.copy()
# 列入れ替えの追跡
col_order = np.arange(N)
for i in range(M):
# ピボット探索
pivot_found = False
for j in range(i, N):
if H_work[i, j] == 1:
# 列入れ替え
if j != i + K:
target = i + K
H_work[:, [j, target]] = H_work[:, [target, j]]
col_order[[j, target]] = col_order[[target, j]]
pivot_found = True
break
if not pivot_found:
# ピボットが見つからない場合はi行目以降を探索
for j in range(K, N):
if j != i + K and H_work[i, j] == 1:
H_work[:, [j, i + K]] = H_work[:, [i + K, j]]
col_order[[j, i + K]] = col_order[[i + K, j]]
pivot_found = True
break
if not pivot_found:
continue
# 行操作
for k in range(M):
if k != i and H_work[k, i + K] == 1:
H_work[k, :] = (H_work[k, :] + H_work[i, :]) % 2
# H_work = [P | I_M] の形になったはず
P = H_work[:, :K]
# 生成行列 G = [I_K | P^T]
c_reordered = np.zeros(N, dtype=int)
c_reordered[:K] = u
c_reordered[K:] = P.T @ u % 2 # パリティビット
# 元の列順序に戻す
c = np.zeros(N, dtype=int)
for i in range(N):
c[col_order[i]] = c_reordered[i]
return c
この符号化関数は、パリティ検査行列 $\bm{H}$ をガウス消去法で $[\bm{P} | \bm{I}_M]$ の形に変換し、生成行列 $\bm{G} = [\bm{I}_K | \bm{P}^T]$ を使って組織的な符号化を行います。GF(2) 上の演算なので、加算はXOR、乗算はANDです。
BP復号器の実装
タナーグラフ上でのBP復号を実装します。
def bp_decode(y, H, sigma2, max_iter=50):
"""
信念伝搬（BP）復号器
y: 受信信号（実数値、BPSK: 0->+1, 1->-1）
H: パリティ検査行列
sigma2: 雑音分散
max_iter: 最大反復回数
戻り値: 推定符号語
"""
M, N = H.shape
# チャネルLLR
Lch = 2.0 * y / sigma2
# エッジの接続情報を構築
# check_neighbors[m] = チェックノードmに接続された変数ノードのリスト
# var_neighbors[n] = 変数ノードnに接続されたチェックノードのリスト
check_neighbors = [np.where(H[m, :] == 1)[0] for m in range(M)]
var_neighbors = [np.where(H[:, n] == 1)[0] for n in range(N)]
# メッセージの初期化
# L_v2c[m][n]: 変数ノードn → チェックノードm へのメッセージ
L_v2c = {}
for m in range(M):
for n in check_neighbors[m]:
L_v2c[(m, n)] = Lch[n]
# L_c2v[m][n]: チェックノードm → 変数ノードn へのメッセージ
L_c2v = {}
for m in range(M):
for n in check_neighbors[m]:
L_c2v[(m, n)] = 0.0
for iteration in range(max_iter):
# チェックノード更新
for m in range(M):
neighbors = check_neighbors[m]
for n in neighbors:
# n を除く隣接変数ノードのメッセージの積
product = 1.0
for n_prime in neighbors:
if n_prime != n:
val = L_v2c[(m, n_prime)]
product *= np.tanh(val / 2.0)
# クリッピング（数値安定性）
product = np.clip(product, -1 + 1e-15, 1 - 1e-15)
L_c2v[(m, n)] = 2.0 * np.arctanh(product)
# 変数ノード更新
for n in range(N):
neighbors = var_neighbors[n]
for m in neighbors:
# m を除く隣接チェックノードからのメッセージの和
total = Lch[n]
for m_prime in neighbors:
if m_prime != m:
total += L_c2v[(m_prime, n)]
L_v2c[(m, n)] = total
# 事後LLR（仮判定用）
L_total = Lch.copy()
for n in range(N):
for m in var_neighbors[n]:
L_total[n] += L_c2v[(m, n)]
# 硬判定
c_hat = (L_total < 0).astype(int)
# パリティ検査
syndrome = H @ c_hat % 2
if np.sum(syndrome) == 0:
return c_hat
return c_hat
このBP復号器は、前のセクションで解説した変数ノード更新とチェックノード更新を交互に実行します。辞書 L_v2c と L_c2v でエッジごとのメッセージを管理し、各反復後にパリティ検査を行って早期停止を判定しています。np.tanh と np.arctanh はチェックノード更新のtanh則をそのまま実装したもので、数値安定性のためにクリッピングを入れています。
BERシミュレーション
LDPC符号の性能をシミュレーションします。
import numpy as np
import matplotlib.pyplot as plt
def simulate_ldpc_ber(N=504, dv=3, dc=6, EbN0_dB_list=None,
num_frames=50, max_iter=50):
"""
LDPC符号のBERシミュレーション
N: 符号長
dv, dc: 正則LDPC符号のパラメータ
EbN0_dB_list: Eb/N0 [dB] のリスト
num_frames: フレーム数
max_iter: BP復号最大反復回数
"""
if EbN0_dB_list is None:
EbN0_dB_list = np.arange(1.0, 5.1, 0.5)
R = 1.0 - dv / dc # 符号化率
M = N * dv // dc
K = N - M
# パリティ検査行列の生成（全フレーム共通）
H = make_regular_ldpc(N, dv, dc)
ber_list = []
for EbN0_dB in EbN0_dB_list:
EbN0 = 10 ** (EbN0_dB / 10.0)
sigma2 = 1.0 / (2.0 * R * EbN0)
sigma = np.sqrt(sigma2)
total_bits = 0
total_errors = 0
for frame in range(num_frames):
# 情報ビット生成
u = np.random.randint(0, 2, K)
# LDPC符号化
c = ldpc_encode(u, H)
# BPSK変調: 0 -> +1, 1 -> -1
x = 1.0 - 2.0 * c
# AWGN通信路
noise = sigma * np.random.randn(N)
y = x + noise
# BP復号
c_hat = bp_decode(y, H, sigma2, max_iter=max_iter)
# 情報ビットの誤りカウント（組織符号の先頭Kビット）
# 注: 列入れ替えがあるため、全ビットで評価
errors = np.sum(c != c_hat)
total_errors += errors
total_bits += N
ber = total_errors / total_bits if total_bits > 0 else 0
ber_list.append(ber)
print(f"Eb/N0 = {EbN0_dB:.1f} dB: BER = {ber:.2e}")
return np.array(EbN0_dB_list), np.array(ber_list)
この関数は各 $E_b/N_0$ の値に対して複数フレームの送受信を行い、BERを計算します。LDPC符号の符号化率は $R = 1 – d_v/d_c$ であり、雑音分散をこの符号化率を考慮して設定しています。
# シミュレーション実行
np.random.seed(42)
EbN0_dBs = np.arange(1.0, 5.1, 0.5)
EbN0_dBs_result, ber_result = simulate_ldpc_ber(
N=504, dv=3, dc=6,
EbN0_dB_list=EbN0_dBs,
num_frames=30,
max_iter=50
)
# 符号化なしBPSKの理論BER
from scipy.special import erfc
EbN0_lin = 10 ** (EbN0_dBs / 10.0)
ber_uncoded = 0.5 * erfc(np.sqrt(EbN0_lin))
# BERカーブ描画
plt.figure(figsize=(10, 7))
plt.semilogy(EbN0_dBs_result, ber_result, 'bs-', linewidth=2,
markersize=8, label='LDPC (3,6)-regular, N=504, BP')
plt.semilogy(EbN0_dBs, ber_uncoded, 'r--', linewidth=1.5,
label='Uncoded BPSK')
plt.axvline(x=0.19, color='green', linestyle=':', linewidth=1.5,
label='Shannon Limit (R=1/2)')
plt.axvline(x=1.11, color='orange', linestyle='--', linewidth=1.5,
label='DE Threshold (3,6)-regular')
plt.xlabel('$E_b/N_0$ [dB]')
plt.ylabel('BER')
plt.title('LDPC Code BER Performance (AWGN, BPSK)')
plt.legend(fontsize=11)
plt.grid(True, which='both', alpha=0.3)
plt.xlim([0, 6])
plt.ylim([1e-6, 1])
plt.tight_layout()
plt.show()
上のグラフから、LDPC符号の性能について重要な特徴が読み取れます。
- ウォーターフォール特性: $E_b/N_0 \approx 2 \sim 3$ dBの領域でBERが急激に低下するウォーターフォール領域が確認できます。密度進化の閾値（$(3,6)$-正則LDPC符号で $E_b/N_0 \approx 1.11$ dB）は $N \to \infty$ の理論値ですが、有限の符号長 $N = 504$ でも閾値の近くでBERが急速に低下し始めています
- 符号化利得: 符号化なしBPSKと比較して、数dBの符号化利得が得られています。この利得は符号長 $N$ を大きくするとさらに増大し、ウォーターフォール領域がシャノン限界に近づいていきます
- 有限長効果: $N = 504$ は比較的短い符号長です。$N$ を数千〜数万に増やすと、ウォーターフォール領域がより急峻になり、密度進化の閾値に近づきます。5G NRでは $N$ が最大数千ビットのLDPC符号が使用されています
BP復号の収束過程の可視化
BP復号の反復過程で、LLRメッセージがどのように変化するかを可視化します。
import numpy as np
import matplotlib.pyplot as plt
def visualize_bp_convergence(N=504, dv=3, dc=6, EbN0_dB=3.0,
max_iter=30):
"""
BP復号の反復ごとのBER変化を可視化
"""
R = 1.0 - dv / dc
M = N * dv // dc
K = N - M
EbN0 = 10 ** (EbN0_dB / 10.0)
sigma2 = 1.0 / (2.0 * R * EbN0)
sigma = np.sqrt(sigma2)
H = make_regular_ldpc(N, dv, dc)
# データ生成・符号化・変調・通信路
u = np.random.randint(0, 2, K)
c = ldpc_encode(u, H)
x = 1.0 - 2.0 * c
y = x + sigma * np.random.randn(N)
# BP復号（各反復のBERを記録）
Lch = 2.0 * y / sigma2
check_neighbors = [np.where(H[m, :] == 1)[0] for m in range(M)]
var_neighbors = [np.where(H[:, n] == 1)[0] for n in range(N)]
L_v2c = {}
L_c2v = {}
for m in range(M):
for n in check_neighbors[m]:
L_v2c[(m, n)] = Lch[n]
L_c2v[(m, n)] = 0.0
ber_per_iter = []
syndrome_per_iter = []
for iteration in range(max_iter):
# チェックノード更新
for m in range(M):
neighbors = check_neighbors[m]
for n in neighbors:
product = 1.0
for n_prime in neighbors:
if n_prime != n:
product *= np.tanh(L_v2c[(m, n_prime)] / 2.0)
product = np.clip(product, -1 + 1e-15, 1 - 1e-15)
L_c2v[(m, n)] = 2.0 * np.arctanh(product)
# 変数ノード更新
for n in range(N):
neighbors = var_neighbors[n]
for m in neighbors:
total = Lch[n]
for m_prime in neighbors:
if m_prime != m:
total += L_c2v[(m_prime, n)]
L_v2c[(m, n)] = total
# 事後LLRと判定
L_total = Lch.copy()
for n in range(N):
for m in var_neighbors[n]:
L_total[n] += L_c2v[(m, n)]
c_hat = (L_total < 0).astype(int)
ber = np.mean(c != c_hat)
ber_per_iter.append(ber)
syndrome = H @ c_hat % 2
syndrome_per_iter.append(np.sum(syndrome))
if np.sum(syndrome) == 0:
# 収束後の残りの反復は同じ値
for _ in range(max_iter - iteration - 1):
ber_per_iter.append(ber)
syndrome_per_iter.append(0)
break
# 可視化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
iterations = np.arange(1, len(ber_per_iter) + 1)
ax1.plot(iterations, ber_per_iter, 'bo-', linewidth=2, markersize=6)
ax1.set_xlabel('Iteration')
ax1.set_ylabel('BER')
ax1.set_title(f'BER vs Iteration ($E_b/N_0$ = {EbN0_dB} dB)')
ax1.grid(True, alpha=0.3)
ax1.set_ylim(bottom=-0.001)
ax2.plot(iterations, syndrome_per_iter, 'rs-', linewidth=2, markersize=6)
ax2.set_xlabel('Iteration')
ax2.set_ylabel('Unsatisfied Checks')
ax2.set_title(f'Syndrome Weight vs Iteration ($E_b/N_0$ = {EbN0_dB} dB)')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 実行
np.random.seed(42)
visualize_bp_convergence(N=504, dv=3, dc=6, EbN0_dB=3.0, max_iter=30)
上の2つのグラフから、BP復号の収束過程を明確に観察できます。
- BERの急速な減少（左図）: 最初の数反復でBERが急速に低下し、その後ゼロに収束しています。これは、BP復号の初期段階で外部情報の交換が効果的に機能し、大半のビットが正しく推定されることを示しています
- シンドローム重みの減少（右図）: 充足されていないパリティ検査の数も反復とともに急速に減少します。シンドローム重みがゼロになった時点で、推定された符号語が全てのパリティ検査を満たし、有効な符号語として受け入れられます
- 早期停止の効果: 十分なSNRがあれば、最大反復回数（50回）に達する前にパリティ検査が充足され復号が終了します。これにより、高SNR領域での平均復号遅延が大幅に短縮されます
LDPC符号の実用と設計
5G NRのLDPC符号
5G NR（New Radio）では、データチャネル（PDSCH/PUSCH）にLDPC符号が採用されています。5G NRのLDPC符号は以下の特徴を持ちます。
準巡回（QC: Quasi-Cyclic）構造: パリティ検査行列は、$Z \times Z$ の巡回シフト行列を要素とするブロック行列として構成されます。$Z$ はリフティングサイズと呼ばれ、符号長の柔軟な調整を可能にします。QC構造により、符号化と復号の並列処理が効率的に行えます。
ベースグラフ: 5G NRには2つのベースグラフ（BG1, BG2）が定義されています。BG1は最大情報ビット長8,448ビット（$R$ は最大8/9）に対応し、大容量データ伝送に使われます。BG2は最大情報ビット長3,840ビット（$R$ は最大2/3）に対応し、制御情報や小パケット伝送に使われます。
レートマッチング: 符号化率の調整は、パリティビットのパンクチャリング（間引き）や反復（リピティション）により行われます。幅広い符号化率に対応できるため、通信路の品質に応じた適応変調・符号化（AMC）が可能です。
ストレージでの応用
SSD（ソリッドステートドライブ）のNANDフラッシュメモリは、微細化が進むにつれてセルあたりの信頼性が低下しています。特にTLC（Triple Level Cell）やQLC（Quad Level Cell）では、1セルに3ビットや4ビットの情報を格納するため、閾値電圧の分布が密になり、誤り率が高くなります。
この問題に対して、LDPC符号は強力な誤り訂正能力で対処しています。SSDコントローラに搭載されるLDPC復号器は、ソフト判定（チャネルから複数回読み出して軟情報を取得）とBP復号を組み合わせることで、NANDフラッシュの信頼性を大幅に向上させています。
LDPC符号の設計手法
実用的なLDPC符号の設計では、以下の手法が用いられます。
PEG（Progressive Edge Growth）アルゴリズム: ガースを最大化するようにタナーグラフの辺を1本ずつ追加していく貪欲アルゴリズムです。短い符号長でも良好な距離特性を持つLDPC符号を構成できます。
プロトグラフ（Protograph）ベース設計: 小さなプロトグラフ（テンプレート）を定義し、それをリフティング（コピーと辺の並べ替え）により大きな符号に拡張する方法です。5G NRのLDPC符号もこの方式で設計されています。
EXIT chart最適化: ターボ符号と同様に、LDPC符号にもEXIT chartによる次数分布の最適化が適用できます。2つの復号器の代わりに、変数ノードとチェックノードのEXIT曲線の間にトンネルが開くように次数分布を最適化します。
まとめ
本記事では、LDPC符号の理論について解説しました。
- 疎パリティ検査行列: LDPC符号のパリティ検査行列 $\bm{H}$ は低密度（疎）であり、各行・各列の1の数が符号長に依存しない固定値である。この疎性がBP復号の効率性を保証する
- タナーグラフ: パリティ検査行列を変数ノードとチェックノードからなる二部グラフとして表現する。ガース（最短サイクル長）が大きいほどBP復号の精度が向上する
- 信念伝搬（BP）復号: タナーグラフの辺に沿ってLLRメッセージを伝搬させ、チェックノード更新（tanh則）と変数ノード更新（和則）を交互に反復する。パリティ検査の充足を早期停止条件として利用できる
- 密度進化: 符号長 $N \to \infty$ の極限でのBP復号の性能を解析する理論的ツールであり、復号閾値を計算できる。非正則LDPC符号の次数分布最適化により、シャノン限界に極めて近い性能が達成される
- 広範な実用化: 5G NR、Wi-Fi、SSDなど、現代の通信・ストレージシステムで広く採用されている
LDPC符号は、ターボ符号とともに近シャノン限界符号の双璧をなし、信念伝搬という確率的推論アルゴリズムの工学的応用の代表例です。LDPC符号のBP復号は、ベイズネットワークでの確率推論やグラフニューラルネットワークのメッセージパッシングとも深い関係があり、情報理論と機械学習を結ぶ架け橋となっています。
次のステップとして、以下の記事も参考にしてください。