Fast Walsh-Hadamard Transform

Why Fast Walsh-Hadamard Transform Matters

Two datasets record frequencies of feature masks. You want the distribution of the result when choosing one mask from each dataset and combining them with XOR, OR, or AND. If each dataset has $N$ possible masks, checking every ordered pair takes $\Theta(N^2)$ time.

A suitable bitwise transform turns this convolution into pointwise multiplication: transform both arrays, multiply matching entries, and apply the inverse transform. With $N=2^d$ masks, the work becomes $\Theta(N\log N)$ for $N\ge2$. This is useful when a pair operation acts independently on each bit, as in parity-state counting, feature unions, feature intersections, and Boolean-function analysis.

The name needs one distinction. The Walsh–Hadamard transform diagonalizes XOR convolution. The corresponding OR and AND algorithms use subset and superset zeta transforms with Möbius inverses. They share a butterfly-shaped implementation, but their formulas and inverse normalization differ. This lesson covers all three and traces a complete XOR example before comparing the OR and AND results.

Core Concepts and Terminology

Let $a$ and $b$ be arrays of the same positive power-of-two length $N=2^d$. Index $i$ is a $d$-bit mask, and $a_i$ is its weight or frequency. For a selected operation $\star\in\{\oplus,\lor,\land\}$, define

$$ c_k=\sum_{\substack{0\le i,j\lt N\\i\star j=k}}a_i b_j. $$

This counts ordered cross-dataset pairs. It does not impose that the two selected masks differ, and it does not identify swapped choices. All three operations keep their result within the same mask universe.

For XOR, the natural-order Walsh–Hadamard transform is

$$ \widehat a_s=\sum_{x=0}^{N-1}(-1)^{\operatorname{popcount}(s\land x)}a_x. $$

Each output is a signed sum; popcount counts set bits. On one bit, the transform takes $(u,v)$ to $(u+v,u-v)$. Repeating these two-point transforms across every bit dimension gives the full transform. Applying it twice gives $Na$, so the inverse repeats the butterflies and divides every entry by $N$.

For OR, use subset sums: $\widehat a_S=\sum_{T\subseteq S}a_T$. For AND, use superset sums: $\widehat a_S=\sum_{T\supseteq S}a_T$. Their inverses use subtraction, with no global division. These are the zeta and Möbius transforms from subset aggregation.

A butterfly processes two entries differing in one chosen bit. In the following formulas, $u$ belongs to the lower index, where that bit is clear, and $v$ belongs to the upper index, where it is set. Both symbols mean the values before either write.

How Fast Walsh-Hadamard Transform Works, Step by Step

The implementation uses the prime modulus $M=1{,}000{,}000{,}007$. Every butterfly and product is computed with exact modular arithmetic. Since $N$ is a power of two and $M$ is odd, $\gcd(N,M)=1$, so the XOR inverse exists. Fermat's little theorem computes it as $N^{M-2}\bmod M$.

The transform routine works as follows:

1. Validate the operation name and require a positive power-of-two length. Copy the input while normalizing every value modulo $M$.
2. Start with half-width $h=1$. Partition the array into consecutive blocks of length $2h$.
3. In each block, pair each lower-half entry with the corresponding upper-half entry. Save both old values in temporaries $u$ and $v$ before writing.
4. Apply the pair rule for the selected operation and direction. Reduce additions and subtractions modulo $M$.
5. Double $h$ and repeat until $h=N$.
6. For inverse XOR only, multiply every entry by $N^{-1}\bmod M$. Return the new array.

The convolution wrapper checks equal lengths, transforms both inputs forward, multiplies corresponding entries modulo $M$, and transforms that product inversely. The transform–multiply–inverse structure is also visible in Nyaan's XOR-convolution implementation; the programs below include their own complete implementations.

For XOR input [1, 2, 3, 4], the first pass pairs indices $0$ with $1$ and $2$ with $3$. Its completed result is [3, -1, 7, -1] in signed arithmetic. The second pass pairs indices $0$ with $2$ and $1$ with $3$, producing [10, -2, -4, 0]. Negative values are shown in signed form to make the trace readable; the code stores their modular residues, such as $M-2$ for $-2$.

flowchart TD
    A["Read a, b, and operation"] --> B{"Equal positive power-of-two lengths and valid operation?"}
    B -- No --> C["Reject input"]
    B -- Yes --> D["Forward-transform both arrays with selected pair rule"]
    D --> E["Multiply corresponding transformed entries modulo M"]
    E --> F["Start inverse with half-width h = 1"]
    F --> G{"h less than N?"}
    G -- Yes --> H["Process every block and pair using saved u and v"]
    H --> I["Double h"]
    I --> G
    G -- No --> J{"Operation is XOR?"}
    J -- Yes --> K["Scale every entry by inverse of N modulo M"]
    J -- No --> L["Return convolution"]
    K --> L

Worked Example with a Full Trace

Take $a=[1,2,3,4]$ and $b=[5,6,7,8]$, with mask indices 00, 01, 10, and 11. We first compute XOR convolution. The following intermediate arrays use ordinary signed representatives; reducing each operation modulo $M$ gives the same final residues.

Transform the first input. At half-width $h=1$, the first butterfly uses $(1,2)$ and writes $(3,-1)$, giving [3, -1, 3, 4]. The second uses $(3,4)$ and writes $(7,-1)$, giving [3, -1, 7, -1]. At $h=2$, the pair at indices $0,2$ uses $(3,7)$ and writes $(10,-4)$, giving [10, -1, -4, -1]. The pair at indices $1,3$ uses $(-1,-1)$ and writes $(-2,0)$, giving $\widehat a=[10,-2,-4,0]$.

Transform the second input. Its successive butterfly states are [11, -1, 7, 8], [11, -1, 15, -1], [26, -1, -4, -1], and $\widehat b=[26,-2,-4,0]$. The first two states complete $h=1$; the last two complete $h=2$.

Multiply pointwise. The four products are $10\times26=260$, $(-2)(-2)=4$, $(-4)(-4)=16$, and $0\times0=0$. Thus $p=[260,4,16,0]$.

Apply the XOR butterflies again. From $p$, the successive states are [264, 256, 16, 0], [264, 256, 16, 16], [280, 256, 248, 16], and [280, 272, 248, 240]. Each state follows one butterfly, in the same block order as before.

Normalize. Since $N=4$, divide by $4$ to obtain

$$ a*_{\oplus}b=[70,68,62,60]. $$

The modular implementation multiplies by $4^{-1}\bmod M=250{,}000{,}002$. A direct check of output mask $0$ uses $i\oplus j=0$ exactly when $i=j$, giving $1\times5+2\times6+3\times7+4\times8=70$.

For comparison, OR forward passes give $a$ states [1, 3, 3, 7] and [1, 3, 4, 10], and $b$ states [5, 11, 7, 15] and [5, 11, 12, 26]. Their pointwise product is [5, 33, 48, 260]. The two inverse passes give [5, 28, 48, 212] and $[5,28,43,184]$.

AND forward passes give $a$ states [3, 2, 7, 4] and [10, 6, 7, 4], and $b$ states [11, 6, 15, 8] and [26, 14, 15, 8]. Their product is [260, 84, 105, 32]. The two inverse passes give [176, 84, 73, 32] and $[103,52,73,32]$. Neither OR nor AND performs global scaling.

All three outputs sum to $260=(\sum_i a_i)(\sum_j b_j)$, because each ordered input pair contributes to exactly one output mask.

Implementation in Python, C++, and Java

Each panel provides a transform routine and a convolution wrapper for operation names xor, or, and and. The wrapper requires equal, nonempty power-of-two lengths; a single entry is valid. The routines copy their inputs and normalize signed values, so the caller's arrays are preserved. The common portable input values fit signed $64$-bit integers, although Python and the Java arithmetic can represent larger values. Arrays must fit available memory.

All three programs print the same seven lines: modulus = 1000000007; a = [1, 2, 3, 4]; b = [5, 6, 7, 8]; xor = [70, 68, 62, 60]; or = [5, 28, 43, 184]; and = [103, 52, 73, 32]; xor round trip = [1, 2, 3, 4].

Python uses exact integers, C++ uses long long, and Java uses BigInteger. With this modulus, a product of reduced operands is at most $(M-1)^2=1{,}000{,}000{,}012{,}000{,}000{,}036$, safely within signed $64$-bit range. Additions and normalized subtractions are at most $2M-1$. Each multiplication is reduced immediately. There is no floating-point arithmetic. Save the Java panel as Main.java.

MOD = 1_000_000_007


def transform(values: list[int], kind: str, inverse: bool = False) -> list[int]:
    size = len(values)
    if kind not in ("xor", "or", "and"):
        raise ValueError("kind must be one of: xor, or, and")
    if size == 0 or size & (size - 1):
        raise ValueError("length must be a positive power of two")
    out = [value % MOD for value in values]
    half = 1
    while half < size:
        for base in range(0, size, 2 * half):
            for offset in range(half):
                left = base + offset
                right = left + half
                u, v = out[left], out[right]
                if kind == "xor":
                    out[left], out[right] = (u + v) % MOD, (u - v) % MOD
                elif kind == "or":
                    out[right] = (v - u if inverse else u + v) % MOD
                else:
                    out[left] = (u - v if inverse else u + v) % MOD
        half *= 2
    if inverse and kind == "xor":
        inverse_size = pow(size, MOD - 2, MOD)
        out = [value * inverse_size % MOD for value in out]
    return out


def convolution(a: list[int], b: list[int], kind: str) -> list[int]:
    if len(a) != len(b):
        raise ValueError("input lengths must match")
    left = transform(a, kind)
    right = transform(b, kind)
    product = [u * v % MOD for u, v in zip(left, right)]
    return transform(product, kind, True)


def main() -> None:
    a, b = [1, 2, 3, 4], [5, 6, 7, 8]
    print(f"modulus = {MOD}")
    print(f"a = {a}")
    print(f"b = {b}")
    for kind in ("xor", "or", "and"):
        print(f"{kind} = {convolution(a, b, kind)}")
    print(f"xor round trip = {transform(transform(a, 'xor'), 'xor', True)}")


if __name__ == "__main__":
    main()

#include <cstddef>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>

constexpr long long MOD = 1000000007LL;

long long mod_pow(long long base, long long exponent) {
    long long result = 1;
    base %= MOD;
    while (exponent > 0) {
        if (exponent % 2 == 1) result = result * base % MOD;
        base = base * base % MOD;
        exponent /= 2;
    }
    return result;
}

std::vector<long long> transform(const std::vector<long long>& values,
                                 const std::string& kind, bool inverse = false) {
    const std::size_t size = values.size();
    if (kind != "xor" && kind != "or" && kind != "and") {
        throw std::invalid_argument("kind must be one of: xor, or, and");
    }
    if (size == 0 || (size & (size - 1)) != 0) {
        throw std::invalid_argument("length must be a positive power of two");
    }
    std::vector<long long> out(size);
    for (std::size_t i = 0; i < size; ++i) {
        out[i] = values[i] % MOD;
        if (out[i] < 0) out[i] += MOD;
    }
    for (std::size_t half = 1; half < size; half *= 2) {
        for (std::size_t base = 0; base < size; base += 2 * half) {
            for (std::size_t offset = 0; offset < half; ++offset) {
                const std::size_t left = base + offset, right = left + half;
                const long long u = out[left], v = out[right];
                if (kind == "xor") {
                    out[left] = (u + v) % MOD;
                    out[right] = (u + MOD - v) % MOD;
                } else if (kind == "or") {
                    out[right] = inverse ? (v + MOD - u) % MOD : (u + v) % MOD;
                } else {
                    out[left] = inverse ? (u + MOD - v) % MOD : (u + v) % MOD;
                }
            }
        }
    }
    if (inverse && kind == "xor") {
        const long long inverse_size = mod_pow(static_cast<long long>(size % MOD), MOD - 2);
        for (long long& value : out) value = value * inverse_size % MOD;
    }
    return out;
}

std::vector<long long> convolution(const std::vector<long long>& a,
                                   const std::vector<long long>& b,
                                   const std::string& kind) {
    if (a.size() != b.size()) throw std::invalid_argument("input lengths must match");
    auto left = transform(a, kind);
    const auto right = transform(b, kind);
    for (std::size_t i = 0; i < left.size(); ++i) left[i] = left[i] * right[i] % MOD;
    return transform(left, kind, true);
}

void print_values(const std::string& label, const std::vector<long long>& values) {
    std::cout << label << " = [";
    for (std::size_t i = 0; i < values.size(); ++i) {
        if (i != 0) std::cout << ", ";
        std::cout << values[i];
    }
    std::cout << "]\n";
}

int main() {
    const std::vector<long long> a = {1, 2, 3, 4}, b = {5, 6, 7, 8};
    std::cout << "modulus = " << MOD << '\n';
    print_values("a", a);
    print_values("b", b);
    for (const std::string kind : {"xor", "or", "and"}) {
        print_values(kind, convolution(a, b, kind));
    }
    print_values("xor round trip", transform(transform(a, "xor"), "xor", true));
    return 0;
}

import java.math.BigInteger;
import java.util.Arrays;

public class Main {
    private static final BigInteger MOD = BigInteger.valueOf(1_000_000_007L);

    public static BigInteger[] transform(BigInteger[] values, String kind, boolean inverse) {
        int size = values.length;
        if (!kind.equals("xor") && !kind.equals("or") && !kind.equals("and")) {
            throw new IllegalArgumentException("kind must be one of: xor, or, and");
        }
        if (size == 0 || (size & (size - 1)) != 0) {
            throw new IllegalArgumentException("length must be a positive power of two");
        }
        BigInteger[] out = new BigInteger[size];
        for (int i = 0; i < size; ++i) out[i] = values[i].mod(MOD);
        for (int half = 1; half < size; half *= 2) {
            for (int base = 0; base < size; base += 2 * half) {
                for (int offset = 0; offset < half; ++offset) {
                    int left = base + offset, right = left + half;
                    BigInteger u = out[left], v = out[right];
                    if (kind.equals("xor")) {
                        out[left] = u.add(v).mod(MOD);
                        out[right] = u.subtract(v).mod(MOD);
                    } else if (kind.equals("or")) {
                        out[right] = inverse ? v.subtract(u).mod(MOD) : u.add(v).mod(MOD);
                    } else {
                        out[left] = inverse ? u.subtract(v).mod(MOD) : u.add(v).mod(MOD);
                    }
                }
            }
        }
        if (inverse && kind.equals("xor")) {
            BigInteger inverseSize = BigInteger.valueOf(size).modPow(
                    MOD.subtract(BigInteger.valueOf(2)), MOD);
            for (int i = 0; i < size; ++i) out[i] = out[i].multiply(inverseSize).mod(MOD);
        }
        return out;
    }

    public static BigInteger[] convolution(BigInteger[] a, BigInteger[] b, String kind) {
        if (a.length != b.length) {
            throw new IllegalArgumentException("input lengths must match");
        }
        BigInteger[] left = transform(a, kind, false);
        BigInteger[] right = transform(b, kind, false);
        for (int i = 0; i < left.length; ++i) left[i] = left[i].multiply(right[i]).mod(MOD);
        return transform(left, kind, true);
    }

    public static void main(String[] args) {
        BigInteger[] a = {BigInteger.ONE, BigInteger.valueOf(2),
                BigInteger.valueOf(3), BigInteger.valueOf(4)};
        BigInteger[] b = {BigInteger.valueOf(5), BigInteger.valueOf(6),
                BigInteger.valueOf(7), BigInteger.valueOf(8)};
        System.out.println("modulus = " + MOD);
        System.out.println("a = " + Arrays.toString(a));
        System.out.println("b = " + Arrays.toString(b));
        for (String kind : new String[]{"xor", "or", "and"}) {
            System.out.println(kind + " = " + Arrays.toString(convolution(a, b, kind)));
        }
        System.out.println("xor round trip = "
                + Arrays.toString(transform(transform(a, "xor", false), "xor", true)));
    }
}

Complexity, Correctness, and Invariants

For $N=2^d$, every transform has $d$ stages and $N/2$ butterflies per stage. Copying and normalization cost $\Theta(N)$; inverse XOR additionally computes one modular inverse in $O(\log M)$ arithmetic operations and scales $N$ entries. Thus a convolution takes $O(N(1+\log N)+\log M)$ time; with the fixed modulus this is $\Theta(N\log N)$ for $N\ge2$ and $\Theta(1)$ for $N=1$.

For a fixed length and selected operation, input values do not change the loop structure. Best, average, and worst running times therefore have the same asymptotic bounds. There is no sparsity shortcut in these implementations. The copying APIs and convolution wrapper use $\Theta(N)$ auxiliary storage; an in-place transform kernel alone can use $O(1)$ extra space beyond its array. The arithmetic bounds treat the fixed modulus as constant-sized; variable-size moduli require accounting for multiplication costs.

For XOR correctness, define $\chi_s(x)=(-1)^{\operatorname{popcount}(s\land x)}$. It satisfies $\chi_s(i\oplus j)=\chi_s(i)\chi_s(j)$, so expanding the convolution sum proves $\widehat c_s=\widehat a_s\widehat b_s$. Orthogonality of these sign functions gives $H^2=NI$, proving the inverse scaling. The butterfly factors this transform into one two-point transform per bit dimension. Serre and Püschel's study of Walsh–Hadamard algorithms discusses this staged butterfly structure.

For OR, both input masks lie inside $S$ exactly when their union lies inside $S$. Therefore the subset sum of the OR convolution equals the product of the two subset sums. Möbius inversion recovers the exact union counts. For AND, both masks contain $S$ exactly when their intersection contains $S$; superset sums give the analogous identity. This explains why the three operations need different pair rules.

Variants, Pitfalls, and Optimizations

Preserve both old values. In an XOR butterfly, assigning the sum to the lower entry and then using that updated entry to form the difference corrupts the transform. For $(u,v)=(1,2)$, the correct output is $(3,-1)$; using the new lower value gives $(3,1)$. The diagram below isolates this write dependency.

Use the correct inverse. XOR needs the same sum/difference butterflies and one global factor $N^{-1}$. OR and AND need subtraction butterflies and no global factor. Dividing each XOR pair by two at every inverse stage is an equivalent convention when two is invertible, but do not combine that convention with an additional final division by $N$.

Choose a supported modulus. OR and AND inverses use no division, so they work over any modular ring with normalized subtraction. XOR inversion requires $\gcd(N,M)=1$. For power-of-two lengths, an odd modulus suffices, but Fermat's exponent $M-2$ specifically requires a prime modulus. With an odd composite modulus, use an appropriate inverse algorithm; with an even modulus and $N\gt1$, this XOR inverse is unavailable. The supplied prime avoids that problem.

Pad the mask universe correctly. For masks up to a maximum value $V$, choose a common power-of-two length $N\gt V$ and fill missing frequencies with zero. Unlike polynomial convolution, bitwise convolution does not need an output length equal to the sum of input lengths minus one. The code requires arrays already padded to the same universe; it does not infer padding.

Distinguish pair conventions. Convolving a frequency array with itself counts ordered pairs, including selecting the same item. To count unordered distinct-item XOR pairs from one dataset, first subtract the number of original items from output mask zero, then divide every count by two in a setting where that division is justified. Equal-valued but distinct items still form valid pairs. Such corrections depend on the application's selection model and are not part of the generic convolution routine.

Ordering, sparsity, and repeated work. This code uses natural binary index order. Sequency-ordered Walsh conventions permute transform coordinates and must be used consistently. For many queries against one fixed distribution, cache its forward transform. If only a few masks have nonzero weights, compare sparse pair enumeration against the dense transform. Independent butterflies within a stage also allow parallel processing, but every stage must finish before the next begins.

Practice Problems and Real-World Use

Specify whether mask combination means parity, union, or intersection before choosing a transform. State whether pairs are ordered, whether they come from different datasets, and whether self-selection is permitted. For small universes, direct pair enumeration provides an independent check of the entire output.

In addition to matching direct convolutions, verify that forward followed by inverse recovers every input modulo $M$, each operation's identity distribution behaves correctly, and output sums equal $(\sum a)(\sum b)\bmod M$. Include a one-entry universe, negative weights, modular wraparound, and invalid lengths. These checks test different properties and help catch mistakes that a round trip alone might miss.