Hadamard 変換とbitwise XOR convolution

最近の ABC は 8 問制になって知らない高度典型がたくさん出るようになりました。(ありがたいけど勉強が追いつかず困ったものです) そんな高度典型はあまり知られていないためでしょうか、解説している記事がとても少ないことが多いです。 XOR による畳み込みもその1つです。(出題は 6 問制のときですが)

↓この記事における XOR による畳み込み

judge.yosupo.jp

Hadamard 変換で畳み込みが可能という記述は見かけますが日本語でその証明を詳しく記してある記事はないようでしたのでここにできるだけ詳しく記しておきたいと思います。

ついでに高速 Hadamard 変換や他の bit 演算による畳み込みとの関係性についても紹介したいと思います。

Kronecker 積

まずは、Herdamard変換を説明する上で重要となる Kronecker 積を導入していきたいと思います。

Kronecker 積とは 2 つの行列$A = \lbrace a_{i, j}\rbrace_{n\times m}, B = \lbrace b_{i, j}\rbrace_{p\times q}$に対して 2 項演算$\otimes$を用いて次のように定義されます。

$$ \begin{align} \begin{split} A\otimes B&= \begin{pmatrix} a_{1,1} B & \cdots & a_{1, m} B \\ \vdots & \ddots & \vdots \\ a_{n,1} B & \cdots & a_{n, m} B \\ \end{pmatrix} \\ &= \begin{pmatrix} a_{1, 1}b_{1, 1} & \cdots & a_{1, 1}b_{1, q} & \cdots & a_{1, m}b_{1, 1} & \cdots & a_{1, m}b_{1, q} \\ \vdots & \ddots & \vdots & \cdots & \vdots & \ddots & \vdots \\ a_{1, 1}b_{p, 1} & \cdots & a_{1, 1}b_{p, q} & \cdots & a_{1, m}b_{p, 1} & \cdots & a_{1, m}b_{p, q} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ a_{n, 1}b_{1, 1} & \cdots & a_{n, 1}b_{1, q} & \cdots & a_{n, m}b_{1, 1} & \cdots & a_{n, m}b_{1, q} \\ \vdots & \ddots & \vdots & \cdots & \vdots & \ddots & \vdots \\ a_{n, 1}b_{p, 1} & \cdots & a_{n, 1}b_{p, q} & \cdots & a_{n, m}b_{p, 1} & \cdots & a_{n, m}b_{p, q} \\ \end{pmatrix} \end{split} \end{align} $$

複数の Kronecker 積を簡単に表現するために

$$A_1\otimes A_2 \otimes \cdots \otimes A_n = \bigotimes_{i = 1}^n A_n$$

という記法も導入しておきます。

この Kronecker 積には$A, B, C$を適当なサイズの行列、$k$をスカラーとして以下のような性質があります。

$$ \begin{align} \begin{split} & k(A\otimes B) = (kA)\otimes B = A \otimes (kB) \\ & (A\otimes B) \otimes C = A \otimes (B \otimes C) \\ & (A\otimes B)(C\otimes D) = (AC) \otimes (BD) \\ \end{split} \end{align} $$

これらは要素を書き下すなどして簡単に証明することができます。

2 つ目の式は結合則が成り立つことを示しています。

また、1 つ目と 3 つ目の式をそれぞれ利用して次の式が成り立つことも示すことができます。

$$ \begin{align} \begin{split} \bigotimes_{i = 1}^n{(k_i A_i)} &= \left(\prod_{i = 1}^n{k_i}\right) \left(\bigotimes_{i = 1}^n{A_i}\right) \\ \left( \bigotimes_{i = 1}^n{A_i} \right)\left( \bigotimes_{i = 1}^n{B_i} \right) &= \left( \bigotimes_{i = 1}^n{A_i B_i} \right) \end{split} \end{align} $$

もう一つだけ性質を示しておきます。

$J_{n, m}$を要素がすべて 1 のである${n\times m}$行列、行列に対する 2 項演算子$\circ$を行列の各要素ごとの積($(A\circ B)_{i,j} = a_{i, j} b_{i, j}$)とします。 このとき$A = \lbrace a_{i, j}\rbrace_{n\times m}, B = \lbrace b_{i, j}\rbrace_{p\times q}$に対して、

$A\otimes B = (A\otimes J_{p, q}) \circ (J_{n, m}\otimes B)$

が成り立ちます。これも要素を書き下せば明らかだと思います。

Hadamard 変換

Hadamard 変換は$2^n$個の実数から$2^n$個の実数への線形写像であり Hadamard 行列$H_n$によって記述されます。 信号解析で結構使われるみたいです。

Hadamard 行列

Hadamard 行列は以下のように定義されます。

$$H_n = \begin{cases} J_{1, 1} & n = 0 \\ \frac{1}{\sqrt{2}}\begin{pmatrix}1 & 1 \\ 1 & -1 \end{pmatrix} \otimes H_{n-1} = H_1 \otimes H_{n-1} & otherwise \end{cases}$$

すなわち

$$H_n = \overbrace{H_1 \otimes H_1 \otimes \cdots \otimes H_1}^{n個} = \bigotimes_{i = 1}^n{H_1}$$

です。

具体的に書き下すと

$$H_1 = \frac{1}{\sqrt{2}}\begin{pmatrix}1 & 1 \\ 1 & -1 \end{pmatrix}, H_2 = \frac{1}{2} \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -1 & 1 & -1 \\ 1 & 1 & -1 & -1 \\ 1 & -1 & -1 & 1 \\ \end{pmatrix}, H_3 = \frac{1}{2^{3/2}} \begin{pmatrix} 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{pmatrix} $$

です。

Hadamard 行列の性質

まず、$H_1$が対称行列であることから 任意の Hadamard 行列も対称行列であることがわかります。

さらに、

$$ \begin{align} \begin{split} H_nH_n &= \left(\bigotimes_{i = 1}^n{H_1}\right)\left(\bigotimes_{i = 1}^n{H_1}\right) \\ &= \bigotimes_{i = 1}^n{(H_1H_1)} \\ &=\bigotimes_{i = 1}^n{I_2} = I_{2^n} \end{split} \end{align} $$

であることから、Hadamard 行列は逆行列が自分自身であり、ユニタリ行列でもあることがわかります。

xor による畳み込みとの関係を示す上で重要なものが次の性質です。

$$(H_n)_{i, j} = \frac{1}{2^{n/2}} (-1)^{\mathrm{popcount(i\&j)}}$$

$i\&j$は$i, j$の bitwise-AND、$\mathrm{popcount}$は 2 進数で表わしたときの 1 の個数を表わしています。この性質を示していきます。

まず、Hadamard 行列は次のように分解できます。

$$H_{n} = \left({J_{2^{n-1}} \otimes H_1 \otimes J_{2^{0}}}\right) \circ \left({J_{2^{n-2}} \otimes H_1 \otimes J_{2^{1}}}\right) \circ \cdots \circ \left({J_{2^{n-i}} \otimes H_1 \otimes J_{2^{i-1}}}\right) \cdots \circ \left({J_{2^{0}} \otimes H_1 \otimes J_{2^{n-1}}}\right)$$

これだけでは少しわかりにくいので、$H_3$の場合を図にしました。

これを見ると $J_{2^{3-k}} \otimes H_1 \otimes J_{2^{k - 1}}$ の$i$行$j$ 列目の要素は、$i$と$j$の bitwise-AND の下から$k$桁目の値によって決まっていることが分かります。具体的には値が 0 なら要素の値は $1/\sqrt{2}$ 、 1 なら $-1/\sqrt{2}$ となっています。

一般の$H_n$についても$\left({J_{2^{n-k}} \otimes H_1 \otimes J_{2^{k - 1}}}\right)_{i, j} $は$i\&j$の下から$k$桁目によって決まります。 したがって、$(H_{n})_{i, j}$は$i\&j$の 1 の数だけ符号が反転することになります。

XOR による畳み込み

いよいよ Hadamard 変換による畳み込みについて解説していきます。

長さ$2^n$のベクトル$\boldsymbol{a}, \boldsymbol{b}, \boldsymbol{c}$が

$$c_{k} = \sum_{i, j, i\oplus j = k}{a_i b_j}$$

の関係にあるとします。($\oplus$は bitwise-XOR)このとき、

$$H_n \boldsymbol{c} = 2^{n/2}H_n \boldsymbol{a} \circ H_n \boldsymbol{b}$$

が成立します。

証明

変換後の$k$番目の要素を考えます。

$S = \{ j | (H_n)_{k, j} = 1/2^{n/2}, 0 \le j < 2^{n} \}, T = \{ j | (H_n)_{k, j} = -1/2^{n/2}, 0 \le j < 2^{n} \}$とします。

すると、

$$(H_n \boldsymbol{a})_k = \frac{1}{2^{n/2}}\sum_{i \in S}{a_i} - \frac{1}{2^{n/2}}\sum_{i \in T}{a_i}$$

となります。

$(H_n \boldsymbol{b})_k$についても同様なので、

$$2^{n/2}(H_n \boldsymbol{a})_k \circ (H_n \boldsymbol{b})_k = \frac{1}{2^{n/2}}\sum_{(i, j)\in (S\times S)\cup(T\times T)}{a_i b_j} - \frac{1}{2^{n/2}}\sum_{(i, j)\in (S\times T)\cup(T\times S)}{a_ib_j}$$

ここで、非負整数$x, y$について

$$\mathrm{popcount}( (x\oplus y)\&k) = \mathrm{popcount}(x\& k) + \mathrm{popcount}(y \& k) - 2\mathrm{popcount}( (x\& k)\& (y \& k))$$

が成立する*1ことから、パリティに注目すると

$$(-1)^{\mathrm{popcount}(i\& k)}(-1)^{\mathrm{popcount}(j\& k)} = (-1)^{\mathrm{popcount}( (i\oplus j)\& k)}$$

が得られます。

したがって、

$$(S\times S)\cup(T\times T) = \{(i, j) | (i\oplus j)\& k \in S \}$$ $$(S\times T)\cup(T\times S) = \{(i, j) | (i\oplus j)\& k \in T \}$$

です。

すなわち、$2^{n/2}(H_n \boldsymbol{a})_k(H_n \boldsymbol{b})_k$は$H_n \boldsymbol{c}$そのものであることがわかりました。

(追記)tokusakurai 氏から$H_n$が$n$個の$H_1$の Kronecker 積で表わすことができる事実を使えば$n = 1$のときのみ考えれば十分というという指摘をもらいました。

高速 Hadamard 変換

単純に Hadamard 行列を掛けるだけでは Hadamard 変換の計算量は $\Theta(4^n)$ であり、XOR の畳み込みを行うには愚直に計算した方が速いです。

これを$\Theta(n2^n)$で行うのが高速 Hadamard 変換です。

アルゴリズム

Hadamard 行列は以下のように分解できます。

$$H_n = \prod_{i = 1}^{n}{I_{2^{n-i}} \otimes H_1 \otimes I_{2^{i-1}}}$$

$H_{3}$の場合はわかりやすさのため$1/2^{n/2}$を分けて書くと、

$$H_3 = \frac{1}{2^{n/2}} \begin{pmatrix} 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & -1 \\ \end{pmatrix} \begin{pmatrix} 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & -1 \\ \end{pmatrix} \begin{pmatrix} 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 \\ \end{pmatrix} $$

です。

$$\boldsymbol{a'} = \left(I_{2^{n-i-1}} \otimes H_1 \otimes I_{2^{i}}\right)\boldsymbol{a}$$

は、$i$ bit 目のみが異なる$x, y$に対し、

$$ \begin{pmatrix} a'_{x} \\ a'_{y} \end{pmatrix} = H_1 \begin{pmatrix} a_{x} \\ a_{y} \end{pmatrix} $$

とする操作に当たります。

したがって、分解によって得られる$n$個の行列を掛ける操作は、行列 1 つあたり$\Theta(2^n)$で可能なので全体で$\Theta(n 2^n)$で Hadamard 変換を行うことができました。

この高速 Hadamard 変換と、前に示した Hadamard 変換と XOR の畳み込みの関係より、長さ$2^n$のベクトルに対するXORの畳み込みも$\Theta(n 2^n)$でできることが分かります。

実装例

剰余環では$2^{-n/2}$を掛ける操作はできない場合があるので、高速 Hadamard 変換では$2^{n/2}H_n$を掛ける操作、逆変換では$2^{-n/2}H_n$を掛ける操作を行うことでプログラム中で足し算、引き算、掛け算と整数の割り算しかしなくともよくしています。

2021/11/20追記

上のプログラムで行っている変換を$\tilde{H}_n = 2^{n/2}H_n$とすると、$\boldsymbol{a}, \boldsymbol{b}$をたたみ込んだあとの配列を$\boldsymbol{c}$として、

$$\tilde{H}_n \boldsymbol{c} =\tilde{H}_n \boldsymbol{a} \circ \tilde{H}_n\boldsymbol{b}$$

が成立します。 したがって、何度も畳み込みを行いたいときでもいちいち逆変換を行う必要がありません。 さらに線形性にも注意すると様々な問題に応用できます。(特に Nim が絡む問題 *2

追記終

// fast Walsh–Hadamard transform
template<class T>
std::vector<T> fast_hadamard_transform(std::vector<T> vec) {
    using hadamard_size_type = typename std::vector<T>::size_type;

    auto vec_size = vec.size();

    // check vec_size is power of 2
    assert(((vec_size - 1)&vec_size) == 0);
    
    for(hadamard_size_type i = 1; i < vec_size; i = i << 1) {
        auto mask = ~i;
        for(auto j = i; j < vec_size; j = (j+1)|i) {
            T a = vec[j&mask];
            T &b = vec[j];
            vec[j&mask] += b;
            b = a - b;
        }
    }

    return vec;
}

// inverse fast Walsh–Hadamard transform
template<class VecType>
auto inv_fast_hadamard_transform (VecType &&vec) {
    auto vec_size = vec.size();
    auto &&ret = fast_hadamard_transform(std::forward<VecType>(vec));
    for(auto &i : ret) i /= vec_size;
    return ret;
}

// bitwise xor convolution
// using fast-Walsh–Hadamard-transform
template<class T>
std::vector<T> xor_convolution
    (const std::vector<T> &a, const std::vector<T> &b) {
    using xorconv_size_type = typename std::vector<T>::size_type;

    assert(a.size() == b.size());

    auto vec_size = a.size();
    std::vector<T> &&transa = fast_hadamard_transform(a),
                   &&transb = fast_hadamard_transform(b);
    
    for(xorconv_size_type i = 0; i < vec_size; i++) {
        transa[i] *= transb[i];
    }
    return inv_fast_hadamard_transform(transa);
}

Library Checkerへの提出

AND、OR による畳み込みへの応用

アルゴリズム

まず、AND による畳み込みを考えます。

いきなりですが$G_n$を以下のように定義します。*3 *4

$$G_n = \bigotimes_{i = 1}^n \begin{pmatrix} 1 & 1 \\ 0 & 1 \\ \end{pmatrix} $$

長さ$2^n$のベクトル$\boldsymbol{a}, \boldsymbol{b}, \boldsymbol{c}$が

$$c_{k} = \sum_{i, j, i\& j = k}{a_i b_j}$$

の関係にあるとき、

$$G_n \boldsymbol{c} = G_n \boldsymbol{a} \circ G_n \boldsymbol{b}$$

が成立します。

詳しく証明は記しませんが、

$$(G_n)_{i, j} = \begin{cases} 1 & (i\&j = i) \\ 0 & otherwise \end{cases} $$

であることを利用すると Hadamard 変換のときと同様に証明できます。

また、逆行列

$$G_n^{-1} = \bigotimes_{i = 1}^n \begin{pmatrix} 1 & 1 \\ 0 & 1 \\ \end{pmatrix}^{-1} = \bigotimes_{i = 1}^n \begin{pmatrix} 1 & -1 \\ 0 & 1 \\ \end{pmatrix}$$

となります。

AND による畳み込みが理解できれば OR による畳み込みも bit を反転させたようなものなので、

$$\bigotimes_{i = 1}^n \begin{pmatrix} 1 & 0 \\ 1 & 1 \\ \end{pmatrix}$$

でできることが分かると思います。

実装例

今までのことを踏まえて AND, OR, XOR の畳み込みが全部できるようなライブラリを書きました。

もっと簡潔になると思っていたのですが 100 行を超えてしまいました。*5

template<class T>
class Bitwise_convolution {
    using Vec = std::vector<T>;
    using conv_size_type
    = typename Vec::size_type;

    template<T (*a)(), T (*b)(),
             T (*c)(), T (*d)()>
    static constexpr std::pair<T, T>
    liner_map(const T &x, const T &y) {
        return {a()*x + b()*y, c()*x + d()*y};
    }

    template<std::pair<T, T> 
            (*liner_map)(const T&, const T&)>
    static constexpr Vec 
    kronecker_product(Vec vec) {
        using kp_size_type =
        typename Vec::size_type;

        auto vec_size = vec.size();

        // check vec_size is power of 2
        assert(((vec_size - 1)&vec_size) == 0);
        
        for(kp_size_type i = 1; i < vec_size; i = i << 1) {
            auto mask = ~i;
            for(auto j = i; j < vec_size; j = (j+1)|i) {
                T &a = vec[j&mask];
                T &b = vec[j];
                std::pair<T, T> &&tmp = liner_map(a, b);
                a = tmp.first;
                b = tmp.second;
            }
        }
        return vec;
    }

    static constexpr T conv_zero() { return T(0);  }
    static constexpr T conv_one()  { return T(1);  }
    static constexpr T conv_neg()  { return T(-1); }

    template<class VecType>
    static auto and_map (VecType& vec) { 
        return kronecker_product
               <liner_map
                <conv_one , conv_one ,
                 conv_zero, conv_one >>(vec);
    }
    template<class VecType>
    static auto and_map_inv (VecType& vec) { 
        return kronecker_product
               <liner_map
                <conv_one , conv_neg ,
                 conv_zero, conv_one >>(vec);
    }

    template<class VecType>
    static auto or_map (VecType& vec) { 
        return kronecker_product
               <liner_map
                <conv_one , conv_zero,
                 conv_one , conv_one >>(vec);
    }
    template<class VecType>
    static auto or_map_inv (VecType& vec) { 
        return kronecker_product
               <liner_map
                <conv_one , conv_zero,
                 conv_neg , conv_one >>(vec);
    }

    template<class VecType>
    static auto xor_map (VecType& vec) { 
        return kronecker_product
               <liner_map
                <conv_one , conv_one ,
                 conv_one , conv_neg >>(vec);
    }
    template<class VecType>
    static auto xor_map_inv (VecType& vec) { 
        auto ret = xor_map(vec);
        for(auto &&i : ret) {
            i /= vec.size();
        }
        return ret;
    }

    template<Vec (*Transform)(const Vec&),
             Vec (*Inverse)  (const Vec&)>
    static auto convolution (const Vec &a,
                             const Vec &b) {
        auto vec_size = a.size();
        auto &&transa = Transform(a),
             &&transb = Transform(b);
        for(conv_size_type i = 0; i < vec_size; i++) {
            transa[i] *= transb[i];
        }
        return Inverse(transa);
    }

    public:

    static auto and_convolution(const Vec&a, const Vec&b) {
        return convolution<and_map, and_map_inv>(a, b);
    }
    static auto xor_convolution(const Vec&a, const Vec&b) {
        return convolution<xor_map, xor_map_inv>(a, b);
    }
    static auto or_convolution(const Vec&a, const Vec&b) {
        return convolution<or_map, or_map_inv>(a, b);
    }
};

Library Checkerへの提出 - Bitwise And Convolution

Library Checkerへの提出 - Bitwise Xor Convolution

最後に

最後の方は少し雑になってしましましたが、ここまで読んでいただきありがとうございました。

今後も解説記事の少ない高度典型を理解した場合は何かしら書くかもしれません。*6

この記事に限らず記事中にミスを見つけた場合や、質問等があればご連絡ください。

参考資料

アダマール変換 - Wikipedia

クロネッカー積 - Wikipedia

アダマール変換, 高速 Walsh-Hadamard 変換, xor-畳み込み

更新履歴

2021/11/20更新 実装例に問題例のリンクなどを追記

2022/02/26更新 tokusakurai 氏からの指摘により一部追記、修正

2022/09/06更新 誤植を修正

*1:各bitごとに考えればいいです

*2:問題例 :https://yukicoder.me/problems/no/1753

*3:気持ちとしては (G_n)_i,j において i の k bit目が 1 なら j の k bit目も1です。

*4:これはメビウス変換だかゼータ変換だかと類似の変換となります。

*5:テンプレートを多用したためです

*6:今回で筆が遅いことを自覚したのであまりたくさんはできない気がしますが