量子情報理論の基本：量子誤り訂正（古典線形符号）

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

## はじめに

[前回の記事](https://qiita.com/SamN/items/650549224eb1f9b511fc)で量子誤り訂正がどんな枠組みで成り立っているのかがわかりました。Shorの符号以外にも、いろんな量子誤り訂正があり得るという展望が見えてきたところで、その他にどんな手法があるのかが気になってきます。が、ちょっとその前に「古典線形符号」について理解しておく必要がありそうです。量子誤り訂正の各種方式の開発には、古典線形符号の知見がベースとして使われているとのことです。というわけで、今回は、古典線形符号について勉強します。一通り理解できた後、古典線形符号の動作をPythonプログラムによって確認したいと思います。

この後の節で何を説明しようとしているか、最初に言っておきます。本記事は、大きく分けて「理論の説明」と「シミュレーション」のパートから成ります。「理論の説明」では、古典線形符号の基本概念である「生成行列」「パリティ検査行列」「符号の距離」について説明した後、符号化方式の具体例として「Hamming符号」と「水平垂直パリティ符号」について見ていきます。さらに、「Gilbert-Varshamovの限界」と「双対符号」という若干発展的な話題を取り上げます。「シミュレーション」では、Hamming符号の動作をシミュレーションします。Pythonでの実装例と処理結果を見ながら、確かに誤りが訂正できていることを確認します。

参考にさせていただいたのは、以下の文献です。

1. [ニールセン、チャン「量子コンピュータと量子通信（３）」オーム社（2005年）](https://www.ohmsha.co.jp/book/9784274200090/)
2. [wikipedia - 線形符号](https://ja.m.wikipedia.org/wiki/%E7%B7%9A%E5%9E%8B%E7%AC%A6%E5%8F%B7)
3. [wikipedia - ハミング符号](https://ja.m.wikipedia.org/wiki/%E3%83%8F%E3%83%9F%E3%83%B3%E3%82%B0%E7%AC%A6%E5%8F%B7)
4. [Quantum Native Dojo - 第9章 量子誤り訂正](https://dojo.qulacs.org/ja/latest/notebooks/9_quantum_error_correction.html)
5. [ハミング符号の最小距離は一般に3であることの証明](http://www.umekkii.jp/college/syllabus/03w_digital_media/note05.pdf)
6. [Gilbert-Varshamov限界](http://www-infosec.ist.osaka-u.ac.jp/~yasunaga/appmath5/gvbound.pdf)


## 理論の説明

### 生成行列

まずは、簡単な話からはじめます。古典線形符号の非常に原始的なものとして、もとの情報を繰り返すだけの符号化方式があります（以下、「繰り返し符号」と呼ぶことにします）。例えば、2進数のビット系列の各ビットを、

```math
\begin{align}
0 \rightarrow 000 \\
1 \rightarrow 111 \tag{1}
\end{align}
```

のように3倍に冗長化して伝送します。もともと0だった情報は000になるので、例えば、この中のどれかのビットが反転して010になったとしても、これは000の一つにビット反転があったせいであり、もともとは0だったのだろうと推測できる（つまり誤りを訂正できる）わけです。2ビット反転した場合は、間違って訂正されてしまうのですが、その確率が低ければ、何もしないよりもいいですよね、という方式です。一般に誤り訂正符号というのは、もとのビット情報を冗長化して膨らますというのが基本です。冗長化して膨らますことを「符号化」と呼んでいます。この符号化する方式がいろいろあるわけです。

少し数学的に抽象化して言うと、$k$個の$\\{ 0,1 \\}$からなるビット系列というのは、集合$\\{ 0,1 \\}^{k}$の要素であると言えます。同様に$n$個のビット系列は、集合$\\{ 0,1 \\}^{n}$の要素です。したがって、符号化というのは、$k$次元のベクトル$x \in \\{ 0,1 \\}^{k}$から$n$次元のベクトル$y \in \\{ 0,1 \\}^{n}$への写像に相当します。この写像を線形変換に限定した符号のことを「線形符号(linear code)」と呼びます[^1]。また、線形変換した先の符号によって構成される空間のことを「符号空間」と呼びます。

[^1]: 量子との対比を明示するため、本記事のタイトルは「古典線形符号」としましたが、本文では以後面倒なので単に「線形符号」と言うことにします。

例えば、先程の1ビットを3ビットに冗長化する符号化は、

```math
G =
\begin{pmatrix}
1\\
1\\
1
\end{pmatrix}  \tag{2}
```

という行列で定義される線形変換と言えます。$x=(0), x=(1)$に対して変換行列$G$を適用すると、

```math
G
\begin{pmatrix}
0
\end{pmatrix} =
\begin{pmatrix}
1\\
1\\
1
\end{pmatrix}
\begin{pmatrix}
0
\end{pmatrix}
=
\begin{pmatrix}
0\\
0\\
0
\end{pmatrix}, \space
G
\begin{pmatrix}
1
\end{pmatrix}
 =
\begin{pmatrix}
1\\
1\\
1
\end{pmatrix}
\begin{pmatrix}
1
\end{pmatrix}
=
\begin{pmatrix}
1\\
1\\
1
\end{pmatrix}  \tag{3}
```

となり、確かに、式(1)と同じことをやっていることがわかります。というわけで、繰り返し符号は線形符号です。この$G$のことを、もとの情報から符号を生成する行列という意味で「生成行列(generator matrix)」と呼びます。これ以降、いろんな線形演算が登場しますが、すべての代数演算は2を法として行うこととします。そうしないと、各ビットが$\\{ 0,1 \\}$に収まらなくなるので[^2]。

[^2]: ちなみに、今回の記事では2進数をベースとした符号しか扱いませんが、一般の線形符号理論では何進数をベースとしても良いように理論構築がなされていて、q進数をベースにした符号のことを、特に「q元線形符号」と言うようです。

では、ここで例題をひとつ。以下の生成行列$G$はどんな符号化操作を表しているでしょうか？

```math
\begin{pmatrix}
1 & 0 \\
1 & 0 \\
1 & 0 \\
0 & 1 \\
0 & 1 \\
0 & 1
\end{pmatrix} \tag{4}
```

簡単ですね。

```math
\begin{align}
(0,0)^T \rightarrow (0,0,0,0,0,0)^T  \\
(0,1)^T \rightarrow (0,0,0,1,1,1)^T  \\
(1,0)^T \rightarrow (1,1,1,0,0,0)^T  \\
(1,1)^T \rightarrow (1,1,1,1,1,1)^T  \tag{5}
\end{align}
```

という符号化です。先程は1ビットを3ビットに冗長化しましたが、今回は2ビットを6ビットに冗長化しています。一般に、kビットをnビットに符号化する符号のことを[n,k]符号と言います。最初の例は[3,1]符号、今回の例は[6,2]符号ということになります。容易にわかると思いますが、[n,k]符号の生成行列は$n \times k$行列になります[^3]。

[^3]: 生成行列$G$を$k \times n$行列、符号$x$を横ベクトルにして、符号化操作を$xG$のように記述している文献・テキストもありますが、本記事の「理論の説明」のパートでは、[ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)に従い、符号$x$を縦ベクトル、生成行列を$n \times k$行列として、符号化操作を$Gx$のように記述することにします。あとで説明しますが、「シミュレーション」のパートでは都合により逆にします。が、いましばらくは、もとの情報と符号は縦ベクトルと思っておいてください。

こんな風に適当に生成行列を定めれば、どんな線形符号も自由自在に作れるのでは、と思ったかもしれませんが、それは違います。行列$G$が生成行列であるためには、生成行列を構成する各列ベクトル$\\{ g_1, g_2, \cdots ,g_k\\}$が線形独立になっている必要があります。線形従属だとすると、異なる情報$x$と$x^{\prime}$が同じ符号$y$にマッピングされることがあります。

どういうことか？以下で説明します。

符号$y$はもとの情報$x$から、

```math
y = Gx = \sum_{i=0}^{k} g_{i} x_{i}  \tag{6}
```

のように、生成行列$G$を適用することで生成されます。いま、$x$とは別の$x^{\prime}$が、

```math
y^{\prime} = Gx^{\prime} = \sum_{i=0}^{k} g_{i} x_{i}^{\prime}  \tag{7}
```

のように生成されたとします。もし$\\{ g_i \\}$が線形独立ならば、式(7)から式(6)を引いた式、

```math
y^{\prime} - y = \sum_{i=0}^{k} g_{i} (x_{i}^{\prime} - x_{i})  \tag{8}
```

において、左辺が$0$になるのは、すべての$i$について$x_{i}^{\prime} - x_{i} = 0$のときのみです（線形独立の定義より）。つまり、$y$と$y^{\prime}$が等しくなるのは、$x$と$x^{\prime}$が等しいときのみです。逆に言うと、ベクトル$x$と$x^{\prime}$が異なれば、必ず$y$と$y^{\prime}$は異なります。

一方、もし、$\\{ g_i \\}$が線形従属であるとするとどうなるかを考えてみます。例えば、簡単のため、

```math
g_1 = g_2 + g_3  \tag{9}
```

という関係が成り立っていたとします。これ、線形従属ってことですよね[^4]。とすると、式(6)は、

[^4]: いま、たまたま2番目と3番目の$g_i$を持ってきて、それを加算したものがちょうど1番目の$g_i$に等しいと置きましたが、別にそうでなくても（任意の組み合わせをもってきても）、以降の議論は成り立ちます。

```math
\begin{align}
y &= G x = \sum_{i=0}^{k} g_{i} x_{i}  \\
&= (g_2 + g_3) x_1 + g_2 x_2 + g_3 x_3 + \cdots + g_k x_k \\
&= g_2 (x_1 + x_2) + g_3 (x_1 + x_3) + g_4 x_4 + \cdots + g_k x_k \tag{10}
\end{align}
```

となり、$x=(x_1,x_2,x_3,x_4,\cdots ,x_k)^T$と$x^{\prime}=(\bar{x_1},\bar{x_2},\bar{x_3},x_4,\cdots ,x_k)^T$とが同じ符号を生成することになります（ここで、$\bar{x_i}$は$x_i$のビット反転とします）。というわけで、生成行列を構成する各列ベクトルは線形独立でなければなりません[^5]。

[^5]: これで、[ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の「演習問題10.15」ができたことになります（多分）。

### パリティ検査行列

もとの情報に生成行列を適用することで符号化できたので、次に、その符号に雑音が加わったとして、それを検知することを考えます。まず、具体例で示します。最初に示した[3,1]繰り返し符号の生成行列は式(2)でした。生成された符号に対して、

```math
H =
\begin{pmatrix}
1 & 1 & 0 \\
0 & 1 & 1
\end{pmatrix}  \tag{11}
```

という行列を適用すると、どうなるでしょうか。符号$x$が$(0,0,0)^T$と$(1,1,1)^{T}$のときのみ（すなわち誤りがない場合のみ）、

```math
H x = 0  \tag{12}
```

になります。それ以外は$0$になりません。したがって、これをもって誤りがあったと判定できます。どのビットに誤りがあったかは、符号方式によって異なるアルゴリズムで求めることになります。いまの例の場合、$Hx$が$(1,0)^{T}$の場合、1番目のビット、$(1,1)^{T}$の場合、2番目のビット、$(0,1)^{T}$の場合、3番目のビットに誤りがあったと判定できます。このように、誤りを検査するために構成された行列$H$のことを「パリティ検査行列(parity check matrix)」と呼び、式(12)の関係を満たすものとして定義します。この例の場合、繰り返し符号なので、3ビットごとに多数決をとって判定するやり方でも良いのですが、いろんな線形符号に対応するためには、線形演算で結果をチェックする、とした方が汎用的な言い方になります。

次に例としてあげた[6,2]繰り返し符号の生成行列は式(4)でしたが、これのパリティ検査行列は何でしょうか？、、、1分待ちます。

どうでしょう？わかりましたでしょうか？

では、答えを言います（見えていたと思いますが）。

```math
\begin{pmatrix}
1 & 1 & 0 & 0 & 0 & 0 \\
0 & 1 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 1 & 0 \\
0 & 0 & 0 & 0 & 1 & 1
\end{pmatrix}  \tag{13}
```

ですね[^6]。確かにそうなっていることは、容易に確認できると思います。式(5)の右側の6次元ベクトルを上の行列で変換すると全部$0$になりますが、どこか1つのビットに誤りがあると、$0$にはなりません。どのビットに誤りがあったかは、先程と同じようなルールで判断することが可能です。そして、[n,k]符号の生成行列は$n \times k$行列でしたが、それに対応したパリティ検査行列は$(n-k) \times n$行列です。

少し一般的な言い方をしてみます。符号化された$y$に雑音が加わり、$y^{\prime} = y + e$のように変化したとします。これにパリティ検査行列をかけると、$H y^{\prime} = H y + H e = H e$となり、雑音により加わった誤り情報に関する寄与のみが残ることになります。この$H y^{\prime}$のことを「誤りシンドローム」と呼びます。そして、この誤りシンドロームの計算結果に基づき、具体的に誤り訂正を行うわけです。もし誤りが$j$番目のビットに生じたとするなら、$e$は$j$番目の成分が$1$、残りがすべて$0$のベクトルです。なので、$H y^{\prime}$の計算結果は、ちょうど$H$の$j$列目の列ベクトルに等しくなります。したがって、誤りシンドロームが、$H$の中のどの列に等しいかを判別すれば、何番目のビットに誤りがあるのかがわかります。そのビットを反転すれば、これで誤り訂正ができたことになります。

[^6]: [ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の「演習問題10.17」の答えです。

さて、ここで、何か自分で考えついた生成行列からパリティ検査行列を求めたり、あるいは、逆に、与えられたパリティ検査行列から生成行列を求められると良いです、と思いませんか？（思ってください！）実は、以下のような考え方でできます。

生成行列からパリティ検査行列を得るには、$G$の列と直交する$n-k$個の線形独立ベクトル$y_1, y_2, \cdots , y_{n-k}$を選び、$H$の行が$y_{1}, y_{2}, \cdots , y_{n-k}$となるように設定します。逆に、パリティ検査行列から生成行列を求めるには、$H$の核を張る$k$個の線形独立なベクトル$y_1, y_2, \cdots , y_k$を選び、$G$が$y_1, y_2, \cdots , y_k$という列を持つように設定します。

ん？という感じかもしれないので、少し噛み砕きます。

もとの情報$x$に$G$を適用して符号化した後、そのまま（誤りなしで）$H$を適用すると$0$になります。任意の$x$に対して$0$なので、

```math
HG = 0  \tag{14}
```

という関係が成り立ちます[^7]。与えらた$G$から$H$を求める場合は、式(14)を満たすようにしないといけないので、$G$を構成するすべての列ベクトルと$H$を構成するすべての行ベクトルの内積が全部$0$でないといけません。つまり、直交していないといけません。全体の符号空間は$n$次元で、$G$を構成する列ベクトルは$k$個です。これに直交する線形独立なベクトルを選ぶとすると、必然的にそれは$n-k$個になります。というわけで、これで$G$から$H$を得るイメージができたと思います。次に、与えられた$H$から$G$を求める場合です。まず、$H$の「核」を張る空間を考えます。線形演算子$A$の核というのは、$Ax=0$となるベクトル$x$全体からなる部分空間のことです。いまの場合、$H$を演算して$0$になる空間というのは、$G$によって生成された符号空間そのものなので、$H$の核を張る線形独立なベクトルを求めれば、それが符号空間を張る線形独立ベクトルになっています。全体が$n$次元の空間に対して、$H$は$n-k$行なので$k$個の線形独立ベクトルを得ることができます[^8]。そしてそれを$G$の列ベクトルとして設定してやれば、これで$G$が完成です。先程、パリティ検査行列は$(n-k) \times n$行列と言いましたが、そうなる理由がこれでわかったかと思います。

[^7]: [ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の「演習問題10.18」です。パリティ検査行列の定義から自明だと思います。

[^8]: $Hx=0$を$n$変数で$n-k$個の方程式からなる連立方程式と思えばわかりやすいかもしれません。$n$次元の自由度に$n-k$個の束縛条件がついているので、結局、正味の自由度は$k$次元になります。

これで、この節を終えても良いのですが、パリティ検査行列について、もう一つ大事なことがあります。それは、パリティ検査行列は、いつでも、その効果を温存しながら「標準形」に書き直せるということです。標準形というのは、

```math
H = (A|I_{n-k})  \tag{15}
```

という形の行列です。ここで、$A$は何らかの$(n-k) \times k$行列、$I_{n-k}$は$n-k$次元の単位行列を表しています。それを横に並べてマージしたような行列が式(15)です。$H$がパリティ検査行列なら、符号空間上のベクトル$x$に対して、$H x = 0$が成り立ちます。これを$n$変数からなる$n-k$個の連立方程式と見ると、なぜ式(15)の形に書き直せるのかがわかります。中学の数学を思い出してください。連立方程式を解くときに、各方程式を睨みながら両辺を何かで掛けて足したり引いたりしながら（つまり、加減法です）、次第に簡単な形にもっていったと思いますが、それと同じ操作をイメージしてください。ただし、変数の数に対して方程式の数が少ないので、ばっちり解は求まりません。$n-k$個の変数各々を、他の$k$個の変数の線形和で表したものが答えとして求まります。式(15)を$Hx=0$に代入したものが、まさにそれです。

このようにパリティ検査行列を標準形にすると、生成行列がここからダイレクトに求まります。$I_k$を$k$次元の単位行列として、

```math
G =
\begin{pmatrix}
I_k \\
-A
\end{pmatrix}  \tag{16}
```

とするだけです。$HG = 0$になっているので、確かにこれは式(15)のパリティ検査行列に対する生成行列になっています[^9]。どうでしょう。簡単ですよね。先程、演算子の核とか線形独立とかゴニョゴニョ説明しましたが、こうやれば一瞬です！

[^9]: 式(16)の$-A$は、今回のように2進数を前提にするなら、つまり2元線形符号を考えるなら、マイナスを取っ払って$A$としても良いです。計算はすべて2を法としているので。

それから、このように作成した生成行列を使えば、もとの情報は符号化された$n$次元ベクトルの最初の$k$個そのものです。先の手続きで誤り訂正ができたとすれば、その最初の$k$個を取り出すだけで、もとの情報を復元することができます（簡単！）。

### 符号の距離

さて、生成行列とパリティ検査行列の説明ができたところで、次に、もう一つ重要な概念である「符号の距離」について説明します。

符号空間$\\{ 0,1 \\}^{n}$上の2つのベクトル$y,y^{\prime}$に距離$d(y,y^{\prime})$を定義することができます。$y$と$y^{\prime}$の各要素を比べて、異なるものの総計を距離とします（ただし、2を法とする計算ではなく、普通に和を計算します）。このように定義された距離を「Hamming距離」と呼びます。また、符号$y$の「Hamming重み」$wt(y)$は、$0$ベクトルからの距離として、$wt(y) \equiv d(y,0)$として定義されます。

$y = G x$で規定される符号化で得られる符号空間$C$の「符号の距離」$d(C)$を、このHamming距離を使って、以下のように定義します。

```math
d(C) = \min_{x,y \in C, x \neq y} d(x,y)  \tag{17}
```

つまり、任意に選んだ2つの符号のHamming距離の最小値が符号の距離です（なので、符号の「最小距離」とも呼ばれたりします）。$d(x,y) = wt(x+y)$が成り立つので[^10]、式(17)は、

[^10]: $d(x,y)$は$x,y$の各要素を比べて異なるものの総計でした。これは、$x+y （mod \space 2）$をやってから、要素として含まれる$1$の個数をカウントしたものに等しいです。一方、$wt(x+y)$は定義より、まさにその要素として含まれる$1$の個数です。したがって、$d(x+y) = wx(x+y)$が成り立ちます。

```math
d(C) = \min_{x \in C, x \neq 0} wt(x)  \tag{18}
```

と書き直しても、同じことです。

符号の距離は符号空間の特徴を表すとても重要な特徴量なので、符号の距離が$d=d(C)$である[n,k]符号のことを特に[n,k,d]という具合に記載することも多いです。

なぜ、そんなに重要なのかというと、1つには、それによって何ビットまでの誤り訂正が可能であるかということがわかるからです。[3,1]繰り返し符号を思い出していただくと、この符号の距離は、いまの定義からすると3です。3倍に冗長化したので、多数決で誤り訂正すると1ビットの誤りまでであれば間違いなく対応できます。[5,1]繰り返し符号の場合であれば、その距離は5であり、多数決方式で誤り訂正するなら2ビットの誤りまでであれば対応できます。ということを延長すると、符号の距離が$2t+1$の符号であれば、$t$ビットの誤りを間違いなく訂正できます。ということで、符号の距離は誤り訂正可能なビット数と密接に関係しています。一般に、$2t+1$の距離にある符号によって、損なわれた符号$y^{\prime}$を$d(y,y^{\prime}) \leq t$を満たす一義の符号語$y$として単に復号化すると、$t$ビットまでの誤りが訂正できます。

符号の距離に関連して、言っておくべき大事なことが、あと2つあるので、それについて順に説明します。

まず、一つ目です。

与えられた[n,k]符号から符号の距離を推定する方法があります。パリティ検査行列を構成する列ベクトルの中から任意の$d-1$個を選ぶと必ず線形独立になっていて、$d$個を選ぶと線形従属になる場合があるとします。このとき符号の距離は$d$に等しいです[^11]。では、証明してみます。

[^11]: [ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の「演習問題10.20」です。

【証明】

$H$を[n,k]符号のパリティ検査行列とし、これを構成する列ベクトルを$\\{ h_1, h_2, \cdots , h_n\\}$とします。符号空間C上の任意の符号を$x$とすると、$H x = 0$なので、

```math
h_1 x_1 + h_2 x_2 + \cdots + h_n x_n = 0  \tag{19}
```

が成り立ちます。いま、ある正の整数$s (\leq n)$に対して、

```math
\begin{align}
& x_{i_1}, x_{i_2}, \cdots , x_{i_s} \neq 0 \\
& x_{i_{s+1}}, x_{i_{s+2}}, \cdots , x_{i_n} = 0 \tag{20}
\end{align}
```

とおき、$x$を

```math
x = (x_{i_1}, x_{i_2}, \cdots , x_{i_s}, \cdots , x_{i_n})^{T}  \tag{21}
```

のように決めます。もし、$x$が符号空間上の符号であるなら、式(19)に代入して、

```math
h_1 x_1 + h_2 x_2 + \cdots + h_s x_s = 0  \tag{22}
```

が成り立っていないといけません。前提より、${h_i}$のうち、どの$d-1$個をとっても線形独立なので、$s$は$d-1$よりも大きくなければいけません。なぜなら、$s$が$d-1$以下だとすると、線形独立の定義から、$x_1,x_2, \cdots , x_s$は全部$0$でないといけないからです。ということで、$s \geq d$でなければならないことがわかります。$s$は符号の中に含まれる$1$の個数（すなわち、Hamming重み）なので、それが$d$以上でなければならないということは、すなわち、この符号の距離$d(C)$は、$d(C) \geq d$を満たしていないといけません。

いま、$d(C)=d$を証明したいので、Hamming重み$wt(x) = d$となる符号$x$が１つでも存在することが示せれば良いです。前提より、$H$を構成する列ベクトルの中に、線形従属な$d$個の列ベクトルがあるということなので、それを、$h_{i_1}, h_{i_2}, \cdots , h_{i_d}$とおきます。そうすると、

```math
a_{i_1} h_{i_1} + a_{i_2} h_{i_2} + \cdots + a_{i_d} h_{i_d} = 0  \tag{23}
```

となる$a_{i_1}, a_{i_2}, \cdots , a_{i_d} \neq 0$が存在することになります（線形従属の定義より）。この$a_{i_1}, a_{i_2}, \cdots , a_{i_d}$を使って、

```math
x = (0, \cdots , 0, a_{i_1}, 0, \cdots , 0, a_{i_2}, 0, \cdots , 0, a_{i_d}, 0 , \cdots , 0) \tag{24}
```

を定義すると、この$x$は$Hx = 0$を満たすので、いま考えている符号空間の要素と言えます。そして、この$x$のHamming重み$wt(x)$は$d$です。これで、$d(C)=d$が証明できました。（証明終）

次に、大事なことの２つ目です。

[n,k,d]符号は$n-k \geq d-1$という関係を満たしていなければなりません。つまり、$n$と$k$が決まれば、可能な$d$の範囲が自動的に決まります。あるいは、$k$と$d$を決めれば、可能な$n$の範囲が自動的に決まります。この関係式のことを「Singletonの限界」と言います[^12]。では、証明してみます。

[^12]: [ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の「演習問題10.21」です。

【証明】

[n,k]符号のパリティ検査行列は$(n-k) \times n$行列なので、そのランクは、

```math
rank(H) \leq n-k  \tag{25}
```

でなければなりません。そして、ランクというのはその行列を構成する列ベクトル（または行ベクトル）の中で線形独立になっているものの最大個数です[^13]。ということは、$H$の中から$rank(H)+1$個の列ベクトルを選ぶと、必ず線形従属になります。先程証明したことから、符号の距離の最大値は$H$に含まれる線形独立な列ベクトルの最大個数に$1$を加えたものに等しいです[^14]。すなわち、

```math
d \leq rank(H) + 1  \tag{26}
```

です。式(25)と式(26)より、$rank(H)$を消去すると、

```math
n-k \geq d-1 \tag{27}
```

が成立ちます。（証明終）

[^13]: ん？と思われた方は、適当な線線形代数の本を参照してください。

[^14]: このあたり一読で理解するのは難しいかもしれません。先程証明した前提条件「任意の$d-1$個を選ぶと必ず線形独立になっていて、$d$個を選ぶと線形従属になる場合があるとします」という言葉の意味も含めて、じっくり吟味してみてください。「任意の$d-1$個を選ぶと必ず線形独立になっていて」というのは、$d$個を選ぶと必ず線形従属になるということとは違います。$d-1$の最大値が、$H$の中の線形独立なベクトルの最大個数（つまり、$rank(H)$）によって限界づけられているということです。


### 具体例

さて、基本的な説明ができたところで、符号方式の具体例を２つ紹介します。一つは「Hamming符号」、もう一つは「水平垂直パリティ符号」です。

#### Hamming符号

Hamming符号は、整数$r \geq 2$によって特徴付けられる$[2^{r}-1, 2^{r}-r-1]$符号です。$r=3$の場合のパリティ検査行列$H$を以下に示します。

```math
H =
\begin{pmatrix}
0 & 0 & 0 & 1 & 1 & 1 & 1 \\
0 & 1 & 1 & 0 & 0 & 1 & 1 \\
1 & 0 & 1 & 0 & 1 & 0 & 1
\end{pmatrix}  \tag{28}
```

全く規則性がなくランダムな$0,1$に見えるかもしれませんが、そうではありません。列ベクトルを左から順番に見てみると、$1$から$8$までの整数を3ビットの2進数に直した形に見えますよね。[n,k]符号のパリティ検査行列は$(n-k) \times n$行列だったので、Hamming符号のパリティ検査行列は$r \times (2^{r} - 1)$行列です。$r=3$とすると$3 \times 7$行列なので、確かに式(28)はそのようになっています。

これを標準形に直してみます。列ベクトルを入れ替えるだけで済みます[^15]。

[^15]: 符号の並び順の定義を変えることに相当します。

```math
H =
\begin{pmatrix}
0 & 1 & 1 & 1 & 1 & 0 & 0 \\
1 & 0 & 1 & 1 & 0 & 1 & 0 \\
1 & 1 & 0 & 1 & 0 & 0 & 1
\end{pmatrix}  \tag{29}
```

となります。この形になれば、対応する生成行列$G$は一瞬でわかります。

```math
G =
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
0 & 1 & 1 & 1 \\
1 & 0 & 1 & 1 \\
1 & 1 & 0 & 1
\end{pmatrix}  \tag{30}
```

となります。

先程説明したように、符号の距離はパリティ検査行列からわかるのでした。式(28)の形の方がわかりやすいので、そっちを見てください。任意の2つの列ベクトルを選ぶと必ず線形独立になっています。例えば、1列目と2列目を選ぶと線形独立です（し、他のどの2つを選んでもそうです）。が、その次の3列目を選ぶと、線形従属になってしまいます。ということで、この[7,3]-Hamming符号の距離は3ということになります[^16]。

[^16]: [ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の「演習問題10.22」です。

この符号で誤り訂正可能なビット数の最大値は、$3 = 2t + 1$の$t$なので、1ビットの誤りであれば訂正することが可能です、ということがわかります。

標準形になっている式(30)で符号化して、式(29)で誤りシンドロームを計算して、パリティ検査符号の列ベクトルたちと比較することで、何ビット目に誤りがあったかがわかりますので、それをビット反転します[^17]。

[^17]: 標準形ではなく、式(28)を使って誤りシンドロームを計算した方が本当は良いのかもしれません。誤りシンドロームの2進数列をそのまま10進数に直した数字が、誤りがあったビット番号に等しくなるので、実装が簡単そうです。標準形を使って誤りシンドロームを計算する場合、いちいちパリティ検査行列の何列目か判定しないといけないです。と思うかもしれませんが、多分、実用的には予めルックアップテーブルを作るのだと思います（実装がちょっとだけ面倒ですが）。

復号化は、標準形で考えていたので、最初の4ビットを取り出せば、元の情報を復元できます。

#### 水平垂直パリティ検査符号

水平垂直パリティ検査符号は、もとの情報をブロック化して、パリティを保持しておく符号化法です。$\\{ 0,1 \\}$系列において、$1$の数が偶数だったときパリティは$0$（または、偶パリティ）、奇数だったときパリティは$1$（または、奇パリティ）と定義します。このパリティをもとの情報データとは別に保持します。例えば、元の情報が6桁の$\\{ 0,1 \\}$系列とします。このときその情報を$2 \times 3$のブロックに並べます。例えば、

```
1 0 1
0 1 1
```

とします。ここで、縦横の数列を見て、各々のパリティを計算して、右または下に並べます。そうすると、

```
1 0 1 | 0
0 1 1 | 0
-------
1 1 0
```

となります。このトータル10個の$\\{ 0,1 \\}$を1次元系列にして伝送します。元の情報を左上から右下方向にスキャンして、そのあとに行のパリティ部分と列のパリティ部分をくっつけて送るとすると、全体の符号系列は$(1,0,1,0,1,1,0,0,1,1,0)^{T}$になります、この符号は[10,6]符号です。

もし、もとの情報の中の1ビットに誤りが発生したとすると、パリティ部分との不整合が生じます。例えば、一番右隅のビットが反転すると、

```
1 0 1 | 0
0 1 0 | 0
-------
1 1 0
```

こんな状態になります。そうすると、2行目のパリティ計算と3列目のパリティ計算が合わなくなります。合わなくなっている行と列から、どこに誤りが発生しているかがわかります。また、パリティ部分に誤りが発生して、1行目のパリティの値が反転してしまったとします。

```
1 0 1 | 1
0 0 0 | 0
-------
1 0 1
```

こうなりますね。そうすると、列のパリティ計算は全部合っているのに、1行目のパリティだけがおかしな値になっています。この場合、1行目のパリティ自身が反転したのだと判断できます。

こんな考え方で構成された符号なのですが、実はこれも線形符号の一種です。

まず、生成行列について考えてみます。

もとの情報が$(1,0,1,0,0,0)^{T}$で、それにパリティ部分$(1,0,1,0,1)$をマージしたような符号を生成すれば良いです。もとの情報部分はそのまま変化しないので、単位行列を置いておけば良いです。パリティ部分は、もとの情報に、

```math
A =
\begin{pmatrix}
1 & 1 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 1 & 1 \\
1 & 0 & 0 & 1 & 0 & 0 \\
0 & 1 & 0 & 0 & 1 & 0 \\
0 & 0 & 1 & 0 & 0 & 1
\end{pmatrix}  \tag{31}
```

こんな行列をかければ生成できます（もちろん、計算はmod 2ですよ。念のため）。$A$の1行目はもとの情報ブロックの1行目のパリティ計算、$A$の2行目はもとの情報ブロックの2行目のパリティ計算、以下順に、1列目のパリティ計算、2列目のパリティ計算、3列目のパリティ計算をするためのものです。したがって、生成行列$G$は、単位行列とこの$A$を合体させて、

```math
G =
\begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 \\
1 & 1 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 1 & 1 \\
1 & 0 & 0 & 1 & 0 & 0 \\
0 & 1 & 0 & 0 & 1 & 0 \\
0 & 0 & 1 & 0 & 0 & 1
\end{pmatrix}  \tag{32}
```

ということになります。標準形になっているので、パリティ検査行列は簡単ですね。

```math
H =
\begin{pmatrix}
1 & 1 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 1 & 1 & 0 & 1 & 0 & 0 & 0 \\
1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 \\
0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1
\end{pmatrix}  \tag{33}
```

になります。

この$H$を使って誤りシンドロームを計算し、その結果が$H$のどの列に相当するのかを見ることで、誤りビットが特定できるので、そのビットを反転すれば。誤り訂正ができます。後は最初の6ビットを取り出せば、もとの情報が復元できます。

### Gilbert-Varshamovの限界

さて、線形符号について具体例も含めて基本的な理解ができたところで、若干発展的な話題を取り上げてみます。「Gilbert-Varshamovの限界」と呼ばれている面白い定理があるので、証明してみます[^18]。

[^18]: [ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の「演習問題10.23」です（ただし、他の文献を参考に言い方を変えています）。「証明は簡単なので演習に残しておく」と書かれているのですが、全然簡単ではないです。数学の本を読んでいるとよく出会う、あるあるです。「容易にわかるように」とか「自明である」という言葉を信じてはいけません。

【定理：Gilbert-Varshamovの限界】

$0 \leq \delta \leq \frac{1}{2}$に対して、十分大きな符号長$n$をとれば、符号の距離$d = \delta n$である、符号化レート$R \geq 1 - H(\delta)$ の符号が存在します。ここで$H$は2値エントロピーで、$H(x) \equiv -x \log (x) - (1-x) \log (1-x)$で定義されます。

つまり、この定理は、十分大きな$n$をとれば、与えられた符号の距離$d = \delta n$を実現して、かつ、符号化レートが$1 - H(\delta)$以上となる符号が必ず存在する、ということを表しています。では、証明してみます。

【定理の証明】

まず、「Hamming球」というものを定義した上で、それに関する補題を２つ証明しておきます。

＜定義：Hammng球＞

$x \in \\{ 0,1 \\}^{n}$と実数$r$に対して、$x$を中心とした「Hamming球」$B(x,r)$とは、

```math
B(x,r) = \{y \in \{ 0,1 \}^{n}: d(x,y) \leq r \}  \tag{34}
```

で定義されるものとします。ここで$d(x,y)$は、Hamming距離なので、$B(x,r)というのは$要するに、$x$を中心に各要素の違いが$r$以下におさまっている点の集まりと思えば良いです。この点がいくつあるかを数えると、それはHamming球の体積というものになります。Hamming球の体積は実は中心$x$がどこであっても半径$r$さえ決まれば同じ値になります。ん？と思われるかもしれませんが、結局$x$がどんな値であったとしても、全体の要素数$n$個の中から、高々$r$個を取り出して反転させる場合の数がHamming球の体積に等しいです、と考えればわかると思います。その体積$Vol(n,r)$は、

```math
Vol(n,r) = \sum_{j=0}^{r}
\begin{pmatrix}
n \\
j
\end{pmatrix}  \tag{35}
```

で計算できます。

＜補題(1)＞

十分大きな$n$、$0 \leq p \leq \frac{1}{2}$に対して、

```math
Vol(n,pn) \leq 2^{H(p)n}  \tag{36}
```

が成立ちます。

＜補題(1)の証明＞

```math
\begin{align}
\frac{Vol(n,pn)}{2^{H(p)n}} &= \frac{\sum_{j=0}^{pn} \begin{pmatrix} n \\ j \end{pmatrix}}{p^{-pn} \cdot (1-p)^{-(1-p)n}} \\
&= \sum_{j=0}^{pn} \begin{pmatrix} n \\ j \end{pmatrix} p^{pn} (1-p)^{(1-p)n} \\
&= \sum_{j=0}^{pn} \begin{pmatrix} n \\ j \end{pmatrix} (1-p)^{n} (\frac{p}{1-p})^{pn} \tag{37}
\end{align}
```

ここで、$p \leq \frac{1}{2}$で、$\frac{p}{1-p} \leq 1$なので、

```math
\begin{align}
\frac{Vol(n,pn)}{2^{H(p)n}} &\leq \sum_{j=0}^{pn} \begin{pmatrix} n \\ j \end{pmatrix} (1-p)^{n} (\frac{p}{1-p})^{j} \\
&= \sum_{j=0}^{pn} \begin{pmatrix} n \\ j \end{pmatrix} (1-p)^{n-j} p^{j} \\
&= ((1-p)+p)^{n} = 1^{n} = 1 \tag{38}
\end{align}
```

これより、式(36)が導けます。（補題(1)の証明終）

＜補題(2)＞

```math
|C| \geq \frac{2^{n}}{Vol(n,d-1)} \tag{39}
```

が成立ちます。ここで、$|C|$は符号空間$C$に含まれる符号の総数です。

＜補題(2)の証明＞

$C$に含まれない、すべての$x \in \\{ 0,1 \\}^{n}$について、少なくとも1つの符号$c_x \in C$が存在して、

```math
d(x,c_x) \leq d-1 \tag{40}
```

が成り立ちます。なぜなら、成り立っていないとすると、$x$を符号空間$C$に含めても（距離が不変になるので）良いことになり、前提と矛盾するからです。

そうすると、集合$\\{ 0,1 \\}^{n}$の全体は、$c \in C$を中心とする半径$d-1$のすべての球の和集合に等しいと言えます。すなわち、

```math
\{ 0,1 \}^{n} = \bigcup_{c \in C} B(c,d-1) \tag{41}
```

したがって、

```math
\begin{align}
2^{n} &= | \{ 0,1 \}^{n} | = | \bigcup_{c \in C} B(c,d-1)| \\
&\leq \bigcup_{c \in C} | B(c,d-1) | \\
&= |C| \sum_{j=0}^{d-1} \begin{pmatrix} n \\ j \end{pmatrix} \\
&= |C| Vol(n,d-1)  \tag{42}
\end{align}
```

となり、

```math
|C| \geq \frac{2^{n}}{Vol(n,d-1)} \tag{39}
```

が成り立つことがわかります。（補題(2)の証明終）

それでは、本体の証明をします。ここまで来たら、一瞬です。符号化レートを$R$、$d=\delta n$、2値エントロピーを$H(\cdot)$として、補題(1)(2)を使うと、

```math
\begin{align}
R &= \frac{\log |C|}{n} \geq \frac{n - \log Vol(n,d-1)}{n} \\
&\geq \frac{n - H(\delta - 1/n)n}{n} \\
&= 1 - H(\delta - \frac{1}{n}) \\
&\geq 1 - H(\delta)  \tag{43}
\end{align}
```

が成立ちます。（定理の証明終）

この定理は、以下のように言い換えることもできます。

十分大きな$n$に対してある$k$を選ぶと、

```math
\frac{k}{n} \geq 1 - H(\frac{2t}{n})  \tag{44}
```

を満たすような$t$ビットの誤りから守る[n,k]誤り符号が存在します[^19]。符号の距離$d$と誤り可能ビット数$t$の間には$d=2t+1$という関係があることから、わかると思います。

[^19]: [ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)は、この言い方で説明がなされています。


### 双対符号

理論説明の締め括りとして「双対符号」というものを説明しておきます。CSS符号という量子符号を構築する上で重要な鍵となる概念とのことです。

$C$を生成行列$G$とパリティ検査行列$H$を持つ[n,k]符号であるとします。このとき、$C$に双対(dual)な[n,n-k]符号$C^{\perp}$を、生成行列$H^{T}$とパリティ検査行列$G^{T}$を持つものと定義します。つまり、$C$の双対は$C$に含まれる全符号と直交する符号から構成されます。もし、$C \subseteq C^{\perp}$ならば符号$C$は「弱い自己双対」と言い、$C = C^{\perp}$ならば「厳密に自己双対」と言います。

ここで、双対符号に関連して成り立つ２つの性質について説明します。

まず１つ目です。

生成行列$G$を持つ符号は$G^{T}G=0$のとき、そのときに限り弱い自己双対であると言えます[^20]。すなわち、

```math
G^{T}G = 0 \Leftrightarrow C \subseteq C^{\perp}  \tag{45}
```

が成立ちます。それでは、証明してみます。

[^20]: [ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の「演習問題10.24」です。

【証明】

$(\Rightarrow)$

すべての$x \in \\{ 0,1 \\}^{k}$に対して$y=Gx$を行うことで、すべての$y \in C$がもれなく生成されます。任意の$y \in C$を選んで、$C^{T}$のパリティ検査行列を適用すると、$G^{T}y = G^{T}Gx$となります。$G^{T}G = 0$が前提だったので、$G^{T}y=0$です。とすると、任意に選んだ$y$は$C^{\perp}$の要素でもあります。したがって、$C \subseteq C^{\perp}$が成立ちます。

$(\Leftarrow)$

$y \in C$ならば前提より$y \in C^{\perp}$でもあります。したがって、$y$に$C^{\perp}$のパリティ検査行列を作用させると$0$になります。すなわち、$G^{T}y = 0$です。$y=Gx$なので$G^{T}Gx=0$となり、これが任意の$x$に対して成り立っていなければならないので、$G^{T}G=0$です。

（証明終）

次に２つめです。

線形符号$C$に対して、

```math
\begin{align}
x \in C^{\perp} &\Rightarrow \sum_{y \in C} (-1)^{xy} = |C| \\
x \notin C^{\perp} &\Rightarrow \sum_{y \in C} (-1)^{xy} = 0  \tag{46}
\end{align}
```

が成立ちます[^21]。それでは、証明してみます。

[^21]: [ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の「演習問題10.25」です。$xy$と書いてあるところは$x$と$y$の内積です。なので、本当は$x^{T}y$と書くべきですが、$T$は省略しています。

【証明】

まず、前半です。

$C$の生成行列とパリティ検査行列を$G,H$、$C^{\perp}$の生成行列とパリティ検査行列を$H^{T},G^{T}$とします。$y \in C$にマッピングされる前の情報を$a \in \\{ 0,1 \\}^{k}$、$x \in C^{\perp}$にマッピングされる前の情報を$b \in \\{ 0,1 \\}^{n-k}$とします。このとき、$x=H^{T}b, y=Ga$で、$xy = b^{T} HG a = 0$なので、

```math
x \in C^{\perp} \Rightarrow \sum_{y \in C} (-1)^{xy} = \sum_{y \in C} 1 = |C|  \tag{47}
```

が成立ちます。

次に、後半です。

任意の$y \in C$は、$a \in \\{ 0,1 \\}^{k}$と$C$の生成行列$G$を使って、$y = Ga$と書くことができます。ここで、$a$を

```math
a = (a_1, a_2, \cdots , a_k)^{T}  \tag{48}
```

のようにその成分で表し、Gを、

```math
G = (g_1, g_2, \cdots , g_k)  \tag{49}
```

のように列ベクトル$g_1,g_2, \cdots , g_k$を使って表すことにします。そうすると、任意の$x \notin C^{\perp}$と$y$の内積は、

```math
xy = xGa = a_1 (x g_1) + a_2 (x g_2) + \cdots + a_k (x g_k)  \tag{50}
```

と書けます。ここで$(x g_{i})$の少なくともどれか１つは必ず$1$になります。なぜなら、すべて$0$だとすると$x$と$y$が直交することになり$x \notin C^{\perp}$という前提を否定することになるからです。いま仮に$(x g_1)$が$1$だったとします（どの$(x g_{i})$が1だったとしても、以下の議論は成り立ちます）。そうすると、

```math
xy = xGa = a_1 + a_2 (x g_2) + \cdots + a_k (x g_k)  \tag{51}
```

となり、

```math
\begin{align}
\sum_{y \in C} (-1)^{xy}
&= \sum_{a_1, a_2, \cdots , a_k = 0,1} (-1)^{a_1} (-1)^{a_2 (xg_2) + \cdots + a_k (xg_k)} \\
&= \sum_{a_2, \cdots , a_k = 0,1} (1+(-1)) (-1)^{a_2 (xg_2) + \cdots + a_k (xg_k)} \\
&= 0 \tag{52}
\end{align}
```

（証明終）


## シミュレーション

理論の説明が、随分長くなってしまいましたが、一通り完了したので、線形符号の動作をPythonプログラムで確認してみます。「繰り返し符号」「Hamming符号」「水平垂直符号」という３つの符号について理解できたはずので、全部試してみたいところではあるのですが、だいぶ疲れてきたので、一つだけやります。どれでも良かったのですが、なんとなく「Hamming符号」にしてみました。

### Hamming符号の実装

ランダムに情報データを生成して、それをHamming符号によって符号化して、ランダムに選んだ1ビットを反転します。それを誤り訂正して復号して誤りがちゃんと訂正されていることを確認します。全体のPythonコードは以下です。

```Python
import numpy as np

def decimal_to_binlist(decimal, digits): # ex) 6,3 --> [1,1,0]

    bin_str = "{:0{digits}b}".format(decimal, digits=digits)
    return  [int(s) for s in list(bin_str)]
    
def binlist_to_decimal(bin_list): # ex) [0,1,1] --> 3

    return int("".join([str(i) for i in bin_list]), 2)

def make_hamming_matrix(r):

    # parity check matrix (H)
    A = []
    for x in range(1,2**r):
        bin_list = decimal_to_binlist(x, r)
        if sum(bin_list) == 1: continue
        A.append(bin_list)
    A = np.array(A)
    I_H = np.eye(r, dtype=int)
    H = np.concatenate([A, I_H])
    
    # represent integer each row of H matrix (for error correction algorithm)
    H_int = [binlist_to_decimal(row) for row in H]
    
    # generator matrix (G)
    I_G = np.eye(2**r-r-1, dtype=int)
    G = np.concatenate([I_G, A], 1)
    
    return G, H, H_int

def generate_data(k, N):  # random k-bits data

    for _ in range(N):
        yield np.random.randint(2, size=k)

def add_noise(d_in):  # bit flip to one bit (select randomly)

    idx = np.random.randint(len(d_in))
    err = np.array([1 if i == idx else 0 for i in range(len(d_in))])
    d_out = (d_in + err) % 2
    return d_out
    
def correct_error(d_in, H_int):

    d_out = d_in.copy()
    p = (d_out @ H) % 2
    x = binlist_to_decimal(p)
    err_idx = H_int.index(x)
    d_out[err_idx] = (d_out[err_idx] + 1) % 2  # bit flip (recover)
    return d_out
    
if __name__ == '__main__':

    r = 3
    n = 2**r - 1
    k = 2**r - r - 1
    N = 10

    G, H, H_int  = make_hamming_matrix(r)

    print("* input(random) -> encode -> add noise(random 1-bit flip) -> correct -> decode:")
    err_count = 0
    for x in generate_data(k, N):
        y = (x @ G)%2
        y_error = add_noise(y)
        y_correct = correct_error(y_error, H_int)
        x_correct = y_correct[0:k]  # decode (= extract 1st to k-th elements)
        print("{0:} -> {1:} -> {2:} -> {3:} -> {4:}".format(x,y,y_error,y_correct,x_correct))
        
        if sum((x+x_correct)%2) == 1: # if x != x_correct --> 1
            err_count += 1

    err_rate = err_count / N
    print("* error rate = {0:} (count:{1:} / total:{2:})".format(err_rate, err_count, N))
```

何をやっているか、説明します。メイン処理部を見てください。

```Python
r = 3
n = 2**r - 1
k = 2**r - r - 1
N = 10
```

で、パラメータを設定します。rの値を設定すると（上の例では3ですが、任意に変えても良いです）、自動的に[n,k]符号のnとkが決まります。Nで何個のデータをランダムに生成するかを指定します。

```Python
G, H, H_int  = make_hamming_matrix(r)
```

で、Hamming符号の生成行列Gとパリティ検査行列Hを計算します。ここで要注意事項をひとつ。上の理論説明のときは[n,k]符号の生成行列は$n \times k$行列、パリティ検査行列は$(n-k) \times n$行列として、符号の生成を$Gx$のように書きました。が、プログラムの都合により、各々転置行列にします。したがって、符号の生成は$xG$のように逆から掛ける形になります。numpyを使うとベクトルは縦ベクトルではなく横ベクトルにした方が自然に書けるので、こうしました。ちょっと混乱させてしまうかもしれませんが、なんとか切り替えてください（すみません）。

H_intは、パリティ検査行列の各行（理論説明したときの記法では各列）のビット配列を十進数に直したリストです。誤り訂正の計算をするときに、これがあった方が楽なので、ついでに出力しています。

では、関数make_hamming_matrixの定義を見てください。

```Python
# parity check matrix (H)
A = []
for x in range(1,2**r):
    bin_list = decimal_to_binlist(x, r)
    if sum(bin_list) == 1: continue
    A.append(bin_list)
A = np.array(A)
I_H = np.eye(r, dtype=int)
H = np.concatenate([A, I_H])
```

標準形でパリティ検査行列を作ります。つまりr次元の単位行列I_Hとそれ以外の行列Aを作って、numpyのconcatenateで縦方向にマージします。Hamming符号のパリティ検査行列の各行（理論説明のときの記法では各列）は1から2**rまでの整数を2進数列にしたものでした。いま、標準形にしたいので単位行列に相当するもの以外の整数のみを取り出して2進数列にして、それを縦方向に並べて行列Aとします。上のforループでそれをやっています。binlist = decimal_to_binlist(x, r)で十進数xから2進数列を計算してbinlistに代入します。decimal_to_binlistでの処理は関数定義を参照してください。forループが終わったら、numpy行列にして、numpyのeyeで作成した単位行列とマージしてパリティ検査行列Hとします。

```Python
# represent integer each row of H matrix (for error correction algorithm)
H_int = [binlist_to_decimal(row) for row in H]
```

ここで、Hの各行に対応した十進数のリストを計算しています。関数binlist_to_decimalは2進数列を十進数列に変換する関数です。内部処理は関数定義を参照してください。

```Python
# generator matrix (G)
I_G = np.eye(2**r-r-1, dtype=int)
G = np.concatenate([I_G, A], 1)
```

で、生成行列Gを作ります。標準形なので簡単です。先程の行列Aと単位行列を横方向にマージするだけです。

それでは、メイン処理部に戻ります。

```Python
err_count = 0
for x in generate_data(k, N):
    ...
```

で、誤りカウンターerr_countを0に初期化して、forループに入ります。generate_data(k, N)で、kビットのデータをランダムにN個生成しています。その各データをxに代入してループ内の処理を実行します。関数generator_dataの処理内容は、関数定義を参照してください。

以下、ループの中身です。

```Python
y = (x @ G)%2
```

で、xに生成行列Gを作用させ、符号yを生成します。ただしmod 2の計算をしないといけないので'% 2'としています。

```Python
y_error = add_noise(y)
```

で、yにランダムにノイズを加え、y_errorとします。具体的にはランダムにビットを選んで反転させます。add_noiseの処理内容は関数定義を参照してください。

```Python
y_correct = correct_error(y_error, H_int)
```

誤り訂正を実行します。関数correct_errorの中身を見てみます。

```Python
d_out = d_in.copy()
p = (d_out @ H) % 2
x = binlist_to_decimal(p)
err_idx = H_int.index(x)
d_out[err_idx] = (d_out[err_idx] + 1) % 2  # bit flip (recover)
return d_out
```

1行目で入力データd_inをd_outにコピーします。2行目でパリティ検査行列を作用させます。結果は2進数列です。これがパリティ検査行列の何行目に相当するか判別します。そのため、3行目のbinlist_to_decimalで一旦十進数に直しておいて、4行目で、H_intのどの列に相当するかをindexメソッドで求めます。それをerr_idxに代入します。5行目で、d_outのerr_idx番目のビットを反転させます。最後にd_outをリターンします（indexメソッドでいちいち探索していますが、本当はLUTにした方が良いのだと思います。面倒なので実装省略しました）。

では、再びメイン処理部に戻ります。

```Python
x_correct = y_correct[0:k]  # decode (= extract 1st to k-th elements)
```

で、誤り訂正された符号y_correctから、もとの情報を復号します。標準形なので簡単です。最初からk番目までのビットを取り出すだけです。これで、一連の処理が完了したので、

```Python
print("{0:} -> {1:} -> {2:} -> {3:} -> {4:}".format(x,y,y_error,y_correct,x_correct))
```

で、各段階のデータを表示します。

```
if sum((x+x_correct)%2) == 1: # if x != x_correct --> 1
    err_count += 1
```

で、もし誤り訂正できていなかったら、誤りカウントに1を加えます。

最後に、

```Python
err_rate = err_count / N
print("* error rate = {0:} (count:{1:} / total:{2:})".format(err_rate, err_count, N))
```

で、誤り率を表示します。


### 動作の確認

それでは、上のコードを動作させた結果を見てみます。

```
* input(random) -> encode -> add noise(random 1-bit flip) -> correct -> decode:
[1 0 1 1] -> [1 0 1 1 0 1 0] -> [1 0 1 1 1 1 0] -> [1 0 1 1 0 1 0] -> [1 0 1 1]
[1 0 1 1] -> [1 0 1 1 0 1 0] -> [1 1 1 1 0 1 0] -> [1 0 1 1 0 1 0] -> [1 0 1 1]
[0 1 0 1] -> [0 1 0 1 0 1 0] -> [1 1 0 1 0 1 0] -> [0 1 0 1 0 1 0] -> [0 1 0 1]
[0 1 0 1] -> [0 1 0 1 0 1 0] -> [1 1 0 1 0 1 0] -> [0 1 0 1 0 1 0] -> [0 1 0 1]
[1 0 1 0] -> [1 0 1 0 1 0 1] -> [1 0 1 1 1 0 1] -> [1 0 1 0 1 0 1] -> [1 0 1 0]
[0 1 1 1] -> [0 1 1 1 1 0 0] -> [0 1 1 1 1 0 1] -> [0 1 1 1 1 0 0] -> [0 1 1 1]
[0 1 1 1] -> [0 1 1 1 1 0 0] -> [0 1 1 1 0 0 0] -> [0 1 1 1 1 0 0] -> [0 1 1 1]
[0 1 1 1] -> [0 1 1 1 1 0 0] -> [1 1 1 1 1 0 0] -> [0 1 1 1 1 0 0] -> [0 1 1 1]
[1 1 1 1] -> [1 1 1 1 1 1 1] -> [0 1 1 1 1 1 1] -> [1 1 1 1 1 1 1] -> [1 1 1 1]
[0 0 1 0] -> [0 0 1 0 1 1 0] -> [0 0 1 1 1 1 0] -> [0 0 1 0 1 1 0] -> [0 0 1 0]
* error rate = 0.0 (count:0 / total:10)
```

もとの情報はランダムに生成していて、ノイズもランダムなので、実行のたびに入出力データは変わるのですが、いつでも完璧に誤りは訂正できていました。めでたし、めでたし。

## おわりに

[ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)をベースに、他の文献を参考にしながら、行間を丁寧に埋めてみたので、かなりの長文になってしまいました。[ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の演習問題も、気づいたらたら、ほぼ全部制覇していました（正しく答えられているかは、、ですが）。おかげで、だいぶ基本事項の理解が進みました。とは言っても、これは「古典」の誤り訂正の話です。これをベースにした「量子」の話は、さらにもっと深くて広いです（多分）。引き続き、気まぐれ＆マイペースですが、頑張らせていただきます。

以上