2020-03-21

陽性適中率を上げたければ検査陽性に対して繰り返して検査すればよい

以下の記事のような「（有病率が小さければ）むやみに検査を拡大しても陽性適中率が低いので却って問題を引き起こす」という言説をしばしば見かけるが，これはちょっとミスリードで，（検査戦略は別の話として）とりあえず陽性適中率を上げたければ，検査陽性に対して繰り返して検査すればよい． jp.quora.com note.com bunshun.jp

各検査の感度を ${p}$ ，特異度を ${q}$ とし，各検査が独立ならば，1回目の検査陽性に対して追加でもう1回検査した場合，クロス集計表は次のようになる．

	罹患あり	罹患なし
検査陽性	${p^{2}\equiv s_{e}}$	${(1-q)^{2}}$
検査陰性	${1-p^{2}}$	${1-(1-q)^{2}\equiv s_{p}}$

有病率を ${r}$ とすれば陽性適中率 ${\mathrm{PPV}}$ は，

$\begin{align} \mathrm{PPV}&=\frac{rs_{e}}{rs_{e}+(1-r)(1-s_{p})}\\ &=\frac{rp^{2}}{rp^{2}+(1-r)(1-q)^{2}} \end{align}$

なので，例として文春の記事の数値 ${p=0.7, q=0.99, r=0.001}$ を使うと， ${\mathrm{PPV}\approx 0.831}$ (83.1%) となる．2回とも検査陽性に対して，さらにもう1回検査を追加すれば，

	罹患あり	罹患なし
検査陽性	${p^{3}}$	${(1-q)^{3}}$
検査陰性	${1-p^{3}}$	${1-(1-q)^{3}}$

となるので，先ほどの数値を使うと ${\mathrm{PPV}\approx 0.997}$ (99.7%) である．

ちなみに検査陽性・検査陰性に関係なく3回検査して多数決で決める場合は，

	罹患あり	罹患なし
検査陽性	${3p^{2}-2p^{3}}$	${1-3q^{2}+2q^{3}}$
検査陰性	${1-3p^{2}+2p^{3}}$	${3q^{2}-2q^{3}}$

となり，先ほどの数値例では ${\mathrm{PPV}\approx 0.725}$ (72.5%) となる．

2020-03-04

物理量の空間分布を VTK ファイルに出力する [Julia]

数値計算で得られる結果は数値の列でしかないので，何かしらの方法で可視化しないといけないことが多い．2次元の簡単なものなら Excel でもいいが，3次元のデータとなると結構厄介である．幸いなことに ParaView という非常に高機能な可視化ソフトウェアがフリーで公開されており，たいていの可視化は ParaView ひとつで事足りる． www.paraview.org

ParaView が扱うファイルは VTK というフォーマットだが，これにはテキストベースの単純なもの（レガシーフォーマット）と XML ベースのフォーマットがある．自分でレガシーフォーマットを出力するプログラムを書くのも手だが，Python や C++ のライブラリが整備されているのでこれを使うとよい．Julia でも WriteVTK というパッケージがあるので，これを使うと比較的簡単に XML ベースの VTK ファイルを出力できる． github.com

WriteVTK の使い方は上のサイトに書いてあるとおりだが，ひとつ引っかかったのは何らかの（スカラー・ベクトル）量を座標を指定して出力する場合である．VTK ファイルではこのような場合，データセットのフォーマットとして Unstructured Grid（非構造格子）を指定する必要があるが，単に物理量の空間分布を可視化したい場合は構造要素（VTK ではセルという）の設定は不要である．ところが，WriteVTK では必ずセルの設定をしないといけないことになっている（Python, C++ のライブラリではどうだろう）．調べてみると VTK_EMPTY_CELL という型があったので，結局以下のようにすれば通った．

vtkfile = vtk_grid(filename, x, y, z, [MeshCell(VTKCellTypes.VTK_EMPTY_CELL, Vector{Integer}(undef, 0))]);

あとは vtkfile["Velocity", VTKPointData()] = v などとして出力するデータを追加していけばよい．他に気をつけるところとしては，PointData としてベクトルを設定する場合，各ベクトルは各列に入れておく必要がある点だろうか．

2019-12-23

プロセス並列化＋リダクション演算による配列の和の計算 [Julia]

Julia プログラミング memo tips

配列要素の総和を求めるような演算は，数値計算においてしばしば表れる．たとえば Julia で書かれた次のようなコードがあるとする．

# 配列 a, b の確保
sum_a = zero(eltype(a));
sum_b = zero(eltype(b));
for j in 1:N
    for i in 1:M
        sum_a += a(i,j);
        sum_b += b(i,j);
    end
end

このようなコードを並列化する場合，リソース（メモリ）の競合を避けるための対策が必要となる．たとえば，C++/OpenMP でスレッド並列化する場合は，reduction指示子があるので比較的簡単に並列化できる．

// 配列 a, b の確保
double sum_a = 0.0;
double sum_b = 0.0;
#pragma omp parallel for reduction(+:sum_a, sum_b)
for ( int i = 0; i < M; i++) {
    for ( int j = 0; j < N; j++ ) {
        sum_a += a[i][j];
        sum_b += b[i][j];
    }
}

残念ながら，現時点で Julia のスレッド並列にはリダクションが実装されていない．一時変数やアトミック操作を使う手もあるが，なるべく簡単に並列化したい．幸いプロセス並列にはリダクションが実装されているので，上のコードは次のようにすることでプロセス並列化できる．

using Distributed
addprocs(); # 複数プロセス立ち上げ
# 配列 a, b の確保
sum_a, sum_a = @distributed (+) for j in 1:N
    sum_a = zero(eltype(a));
    sum_b = zero(eltype(b));
    for i in 1:M
        sum_a += a(i,j);
        sum_b += b(i,j);
    end
    [sum_a, sum_b];
end

2018-12-28

ジョーンズ計算法によるセナルモン法の説明

物理 memo

複屈折媒質（リターダ）に光が入射すると，直交する偏光成分間に位相差が生じる．この位相差のことをリタデーションといい，実験的に測定する手法としてセナルモン法 (Senarmont method) がある．意外とセナルモン法について書かれた書籍が少ないので，備忘録を兼ねてジョーンズ計算法 (Jones calculus) によりセナルモン法を説明してみる．

セナルモン法の手順

セナルモン法の手順は以下の通りである．割と面倒．

レーザーの光軸上に偏光子を置く
光が透過しないように偏光子の回転角を調整する（偏光面と透過軸を直交させる）
偏光子の手前に ${\lambda/4}$ 板を置き，光が透過しないように回転角を調整して（高速軸を偏光面と一致させ）一旦取り出す
リターダを偏光板の手前に置き，光が透過しないように回転角を調整する（高速軸を偏光面と一致させる）
リターダを ${45^{\circ}}$ だけ回転させる
リターダと偏光子の間に一旦取り出した ${\lambda/4}$ 板を再び入れる
光が透過しないように偏光子の回転角を調整する．その回転角の2倍がリタデーションである

ジョーンズ行列による偏光素子の表記

偏光が絡む問題を数学的に扱う方法としてはジョーンズ計算法とミュラー計算法 (Mueller calculus) があるが，完全偏光だけを扱うのであればジョーンズ計算法が簡単でよい．

リターダのジョーンズ行列

リターダの高速軸と低速軸をそれぞれ ${x, y}$ 軸として，リタデーション ${\phi}$ を生じるリターダのジョーンズ行列 ${\mathbf{J}(\phi)}$ は次の式で表せる．

$\begin{equation} \mathbf{J}(\phi) = \begin{bmatrix} \mathrm{e}^{\mathrm{i}\frac{\phi}{2}} & 0 \\ 0 & \mathrm{e}^{-\mathrm{i}\frac{\phi}{2}} \end{bmatrix} \end{equation}$

高速軸が ${x}$ 軸から ${\theta}$ 傾いたリターダのジョーンズ行列 ${\mathbf{J}(\phi, \theta)}$ は，2次元の回転行列を ${\mathbf{R}(\theta)}$ として ${\mathbf{J}(\phi, \theta) = \mathbf{R}(-\theta)\mathbf{J}(\phi)\mathbf{R}(\theta)}$ で表せるので，

$\begin{equation} \mathbf{J}(\phi, \theta) = \begin{bmatrix} \cos\dfrac{\phi}{2} + \mathrm{i}\sin\dfrac{\phi}{2}\cos 2\theta & \mathrm{i}\sin\dfrac{\phi}{2}\sin 2\theta \\ \mathrm{i}\sin\dfrac{\phi}{2}\sin 2\theta & \cos\dfrac{\phi}{2} - \mathrm{i}\sin\dfrac{\phi}{2}\cos 2\theta \end{bmatrix} \end{equation}$

となる．

${\lambda/4}$ 波長板と偏光子のジョーンズ行列

セナルモン法で用いる ${\lambda/4}$ 板と偏光子のジョーンズ行列を以下に示しておく．

$\begin{align} \mathbf{J}_{\mathrm{QWP}}(\theta) = \mathbf{J}\left(\dfrac{\pi}{2}, \theta\right) &= \begin{bmatrix} \dfrac{\sqrt{2}}{2} + \mathrm{i}\dfrac{\sqrt{2}}{2}\cos 2\theta & \mathrm{i}\dfrac{\sqrt{2}}{2}\sin 2\theta \\ \mathrm{i}\dfrac{\sqrt{2}}{2}\sin 2\theta & \dfrac{\sqrt{2}}{2} - \mathrm{i}\dfrac{\sqrt{2}}{2}\cos 2\theta \end{bmatrix} \\ \mathbf{J}_{\mathrm{POL}}(\theta) &= \begin{bmatrix} \cos^{2}\theta & \sin\theta\cos\theta \\ \sin\theta\cos\theta & \sin^{2}\theta \end{bmatrix} \end{align}$

偏光子の ${\theta}$ は透過軸の ${x}$ 軸からの傾きである．

ジョーンズ計算法によるセナルモン法の説明

入射光は ${[1~0]^{T}}$ 偏光であるとする．上に書いたセナルモン法の手順は結局のところ， ${[1~0]^{T}}$ 偏光の光をリターダ ${\mathbf{J}(\phi, 45^{\circ})}$ ， ${\lambda/4}$ 板 ${\mathbf{J}_{\mathrm{QWP}}(0^{\circ})}$ ，偏光子 ${\mathbf{J}_{\mathrm{POL}}(\theta+90^{\circ})}$ の順に入射させることに相当するので，透過光 ${\mathbf{u}}$ はジョーンズ計算法により，

$\begin{align} \mathbf{u}&=\mathbf{J}_{\mathrm{POL}}(\theta+90^{\circ})\mathbf{J}_{\mathrm{QWP}}(0^{\circ})\mathbf{J}(\phi, 45^{\circ})\begin{bmatrix}1 \\ 0\end{bmatrix}\\ &= \begin{bmatrix} \sin^{2}\theta & -\sin\theta\cos\theta \\ -\sin\theta\cos\theta & \cos^{2}\theta \end{bmatrix} \begin{bmatrix} \mathrm{e}^{\mathrm{i}\frac{\pi}{4}} & 0\\ 0 & \mathrm{e}^{-\mathrm{i}\frac{\pi}{4}} \end{bmatrix} \begin{bmatrix} \cos\dfrac{\phi}{2} & \mathrm{i}\sin\dfrac{\phi}{2} \\ \mathrm{i}\sin\dfrac{\phi}{2} & \cos\dfrac{\phi}{2} \end{bmatrix} \begin{bmatrix}1 \\ 0\end{bmatrix}\\ &= \begin{bmatrix} \mathrm{e}^{\mathrm{i}\frac{\pi}{4}}\sin\theta\left(\cos\dfrac{\phi}{2}\sin\theta-\sin\dfrac{\phi}{2}\cos\theta\right) \\ \mathrm{e}^{\mathrm{i}\frac{\pi}{4}}\cos\theta\left(\sin\dfrac{\phi}{2}\cos\theta-\cos\dfrac{\phi}{2}\sin\theta\right) \end{bmatrix}\\ &= \begin{bmatrix} \mathrm{e}^{\mathrm{i}\frac{\pi}{4}}\sin\theta\sin\left(\theta-\dfrac{\phi}{2}\right) \\ -\mathrm{e}^{\mathrm{i}\frac{\pi}{4}}\cos\theta\sin\left(\theta-\dfrac{\phi}{2}\right) \end{bmatrix} \end{align}$

と表され，光強度は

$\begin{align} |\mathbf{u}|^{2}=\sin^{2}\left(\theta-\dfrac{\phi}{2}\right) \end{align}$

となる．上の手順 7. は ${\sin^{2}(\theta-\phi/2)=0}$ となる ${\theta}$ を見つけることに他ならない．これは当然 ${\theta=\phi/2}$ なので，偏光子の回転角 ${\theta}$ の2倍がリタデーション ${\phi}$ となる．

2018-09-23

シングルピクセルイメージング

物理 memo

イメージングというと普通は CCD や CMOS などのアレイセンサを使うが，フォトダイオードのように空間分化能を持たない点型検出器を使ったイメージング手法を圧縮イメージングとかシングルピクセルイメージング (SPI) という．この方法の何がうれしいかというと，高感度なので SNR が低くても解像可能らしい．他にも点型検出器を使うので普通のイメージセンサが感度を持たないテラヘルツ波などでイメージングできる，結像光学系が必要ない，といった利点もある．シングルピクセルイメージングの基本的なアイデアは空間的に変調された照明光を使う点にある．

問題の定式化

照明光の空間パターンを1次元に並べたベクトルを ${\mathbf{a}_{i}}$ ，イメージング対象物の透過率/反射率の空間分布を1次元に並べたベクトルを ${\mathbf{x}}$ ，透過光/反射光の強度を ${b_{i}}$ とすると，

$\begin{equation} \langle\mathbf{a}_{i}, \mathbf{x}\rangle = \mathbf{a}_{i}^{T}\mathbf{x} = b_{i} \end{equation}$

となるのはいいだろう．あるいは連続系で表した方が分かりやすいかもしれない．

$\begin{equation} \iint a_{i}(\mathbf{r})x(\mathbf{r})\mathrm{d}^{2}r =b_{i} \end{equation}$

添字の ${i=1, 2, \dots}$ は測定回を表す．空間パターンの生成には空間光変調器 (SLM) やデジタルミラーデバイス (DMD) などが使われるようである． ${\mathbf{a}_{i}^{T}}$ を列方向に並べた行列を ${\mathbf{A}}$ ， ${b_{i}}$ を並べたベクトルを $\mathbf{b}$ とすると，まとめて

$\begin{equation} \mathbf{A}\mathbf{x} = \mathbf{b} \end{equation}$

と書ける．何のことはない，線形方程式系である．ただしこの線形方程式系は普通は劣決定系である．というのも，対象物の解像度を ${N^{2}}$ ，測定回数を ${M}$ とすると， ${\mathbf{A}}$ は ${M}$ 行 $N^{2}$ 列となるが，以下に説明するように ${N^{2}>M}$ でもイメージング可能なのである．それに（優）決定系にしようと思えば，解像度を少し高くしただけで多数の測定が必要になることはすぐ分かる．

解析的アプローチ

劣決定系の線形方程式系の解法としてすぐ思い浮かぶのは，正規方程式から ${\|\mathbf{A} \mathbf{x}-\mathbf{b}\|_{2}}$ を最小化する ${\mathbf{x}}$ 求めるという方法である．この場合，普通 ${\mathbf{A}^{T}\mathbf{A}}$ が正則行列ではないので，最急降下法や共役勾配法等の反復法が使われる．反復法で最小二乗解を求める方法はそこそこ上手くいくが，測定結果にノイズが乗っている場合いわゆる「過剰適合」を生じる．そこでもう少し凝った方法として正則化を使う手がある．スパースモデリングに基づく方法だと， ${\mathbf{x}}$ がある種の基底変換でスパース（ほとんどがゼロ成分）になるという性質を利用する．基底変換の行列を ${\mathbf{D}}$ ，変換後のベクトル（基底に関する成分）を ${\mathbf{c}}$ とすると， ${\mathbf{D}\mathbf{x}=\mathbf{c}}$ となるが， ${\mathbf{c}}$ の ${\ell_{1}}$ -norm を最小化してやるのである（本当は ${\ell_{0}}$ -norm を最小化すべきだが，この場合は組合せ最適化の問題となってしまう．なお， ${\ell_{0}}$ -norm は通常の意味でのノルムではない）．問題としては次の様に定式化できる．

$\begin{equation} \mathrm{min}~\|\mathbf{c}\|_{1}\quad\mathrm{s. t.}\quad\mathbf{A}\mathbf{x}=\mathbf{b},~\mathbf{D}\mathbf{x}=\mathbf{c} \end{equation}$

この手の問題は Alternating Direction Method of Multipliers (ADMM) で割と効率的に解ける（「割と」と書いたのは逆行列を求める必要があるため）．ADMM では次の拡張ラグランジュ関数の最小化問題を解く．

$\begin{equation} \mathrm{min}~L=\|\mathbf{c}\|_{1}+\mathbf{h}^{T}_{1}(\mathbf{D}\mathbf{x}-\mathbf{c})+\frac{\mu_{1}}{2}\|\mathbf{D}\mathbf{x}-\mathbf{c}\|_{2}+\mathbf{h}^{T}_{2}(\mathbf{A}\mathbf{x}-\mathbf{b})+\frac{\mu_{2}}{2}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_{2} \end{equation}$

ここで ${\mathbf{h}}$ はラグランジュ乗数， ${\mu}$ は調整パラメータである． ${\mathbf{D}}$ には ${\mathbf{x}}$ の勾配を取る行列を使うとよい結果が得られるようである（この場合， ${\|\mathbf{c}\|_{1}}$ が Total Variation となる）．

統計光学的アプローチ

上記とは違ったアプローチとして ${\mathbf{A}}$ と $\mathbf{b}$ の相関を利用する方法がある．詳しい原理はこちらを見てもらうとして，手順としては次の計算を行うだけである．

$\begin{equation} \mathbf{g} = \frac{1}{M}\mathbf{A}^{T}(\mathbf{b}-\langle\mathbf{b}\rangle) \end{equation}$

ここで ${\langle\rangle}$ はアンサンブル平均である．この計算は ${\mathbf{a}_{i}}$ のゆらぎ ${\Delta\mathbf{a}_{i}=\mathbf{a}_{i}-\langle\mathbf{a}_{i}\rangle}$ と， ${b_{i}}$ のゆらぎ ${\Delta b_{i}=b_{i}-\langle b_{i}\rangle}$ の2次コヒーレンス ${\langle\Delta\mathbf{a}_{i}\Delta b_{i}\rangle}$ を求めることに相当する．この方法は特にゴーストイメージングと呼ばれているようである．反復計算を必要としないのはよいが， ${M\to\infty}$ での誤差の漸近挙動があまりよろしくない．ただしノイズには強い．

直交基底を使う方法

これまで係数行列（照明光の空間パターン） ${\mathbf{A}}$ について特に言及しなかったが，こちらを工夫する手もある．何も考えずランダムなパターンを用いるという手もあるが，直交基底を用いると空間周波数ごとに測定できて効率的である．実際，アダマール行列に基づいた空間パターンを使うことでゴーストイメージングより鮮明度の高い像が得られるようである．あるいは事後的に ${\mathbf{A}}$ を直交化するという方法も考えられる．グラム・シュミットの方法で ${\mathbf{A}}$ を正規直交化すると（ ${\mathbf{b}}$ にも同じ操作を施す），かなりよい像が得られるようである．ただしこの方法ではノイズに弱くなる模様．

回折の影響

変調器で生成された照明光は回折するので，必要に応じてフレネル回折やら角スペクトル法で回折パターンを求める必要がある．ただしリレー光学系を組み込めばこの処理は必要ない．

参考文系

あとで書く．

2018-09-19

ベクトル化した2次元データに対する離散フーリエ変換行列

数学 memo

機械学習では画像のような2次元のデータをベクトル化して1次元のデータとして扱うことが多い．このようにベクトル化した2次元データに対する離散フーリエ変換 (DFT) 行列を定義できれば便利である．

長さ ${N}$ のベクトル ${\mathbf{x}}$ に対する1次元の DFT は， ${\mathbf{W}_{N}}$ を ${N{\times}N}$ の DFT 行列として次の様に表せる．

$\begin{equation} \mathrm{DFT}_{1}(\mathbf{x}) = \mathbf{W}_{N}\mathbf{x} \end{equation}$

同様に ${M}$ 行 ${N}$ 列の2次元データ ${\mathbf{X}}$ に対する2次元の DFT は次の様に表せる．

$\begin{equation} \mathrm{DFT}_{2}(\mathbf{X}) = \mathbf{W}_{M}\mathbf{X}\mathbf{W}_{N} \end{equation}$

さて，2次元のデータ（行列）をベクトル化する操作を ${\mathrm{vec}(\cdot)}$ と表すと，次の関係が成り立つ．

$\begin{equation} \mathrm{vec}(\mathbf{AXB}) = \left(\mathbf{B}^{T}\otimes\mathbf{A}\right)\mathrm{vec}(\mathbf{X}) \end{equation}$

ここで ${\otimes}$ はクロネッカー積である．2番目の式とあわせて考えると，次の関係を得る*1．

$\begin{align} \mathrm{vec}(\mathrm{DFT}_{2}(\mathbf{X})) &= \mathrm{vec}(\mathbf{W}_{M}\mathbf{X}\mathbf{W}_{N})\\ &= \left(\mathbf{W}_{N}^{T}\otimes\mathbf{W}_{M}\right)\mathrm{vec}(\mathbf{X})\\ &= \left(\mathbf{W}_{N}\otimes\mathbf{W}_{M}\right)\mathrm{vec}(\mathbf{X})\\ \end{align}$

というわけで，ベクトル化した2次元データに対する DFT 行列は ${\mathbf{W}_{N}\otimes\mathbf{W}_{M}}$ で定義できることが分かる．ちなみに ${\mathbf{X}}$ を ${128{\times}128}$ ，DFT 行列の要素を std::complex< double > とすると，この DFT 行列は 4 GB のメモリを必要とする．

*1:DFT 行列の性質 ${\mathbf{W}^T=\mathbf{W}}$ に注意．

2018-08-28

シルベスターの方法で生成したアダマール行列からウォルシュ型アダマール行列を得る方法 [C++]

数学プログラミング memo

信号処理や統計学で使われるアダマール行列は，以下のようにシルベスターの方法で再帰的に生成するのが（おそらく）一般的．

$\begin{align} \mathbf{H}_{2} &=\left[ \begin{array}{rr} 1 & 1 \\\ 1 & -1 \end{array}\right]\\\ \mathbf{H}_{2^{k}} &= \mathbf{H}_{2}\otimes\mathbf{H}_{2^{k-1}} \end{align}$

ここで $\otimes$ はクロネッカー積を表す．例として $\mathbf{H}_{8}$ を示す．

$\begin{equation} \mathbf{H}_{8} =\left[ \begin{array}{rrrrrrrr} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\\ 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1\\\ 1 & 1 & -1 & -1 & 1 & 1 & -1 & -1\\\ 1 & -1 & -1 & 1 & 1 & -1 & -1 & 1\\\ 1 & 1 & 1 & 1 & -1 & -1 & -1 & -1\\\ 1 & -1 & 1 & -1 & -1 & 1 & -1 & 1\\\ 1 & 1 & -1 & -1 & -1 & -1 & 1 & 1\\\ 1 & -1 & -1 & 1 & -1 & 1 & 1 & -1 \end{array}\right] \end{equation}$

さて，この方法で生成したアダマール行列について，各行の符号反転の回数（周波数）を小さい順に並べ替えた方が都合がいい場合がある．そのような行列をウォルシュ型アダマール行列と呼んだりするが，次のような方法で並べ替えることができる．

/**
 * @brief ビット逆転（32 bit 用）
 * @param[in] a 32 bit 符号なし整数
 * @return ビット逆転した 32 bit 符号なし整数
 */
unsigned int bitreverse( const unsigned int a ) {
    unsigned int x = a;
    x = ( ( ( x & 0xaaaaaaaa ) >> 1 ) | ( ( x & 0x55555555 ) << 1 ) );
    x = ( ( ( x & 0xcccccccc ) >> 2 ) | ( ( x & 0x33333333 ) << 2 ) );
    x = ( ( ( x & 0xf0f0f0f0 ) >> 4 ) | ( ( x & 0x0f0f0f0f ) << 4 ) );
    x = ( ( ( x & 0xff00ff00 ) >> 8 ) | ( ( x & 0x00ff00ff ) << 8 ) );
    return( ( x >> 16 ) | ( x << 16 ) );
}

/**
 * @brief 並び替えのための配列生成
 * @param[in] k 行列の次数 (2^k)
 * @param[out] v 並び替えのための配列
 */
void permutation( const unsigned int k,
                  std::vector< unsigned int > &v ) {
    const unsigned int N = 1 << k;
    v.resize(N);
    for ( unsigned int i = 0; i < N; i++ ) {
        unsigned int j = bitreverse(i) >> ( 31 - k );
        j = j ^ ( j >> 1 );
        j = ( j << ( 32 - k ) ) >> ( 32 - k );
        v[i] = j;
    }
}

たとえばpermutation( 3, v )とした場合はv = [0 4 6 2 3 7 5 1]となるので，この順に行を並べ替えれば8次のウォルシュ型アダマール行列を得られる．上のプログラムは当然ながらsizeof( unsigned int ) == 4という前提．