Fast Walsh Transform

Introduction

The Walsh transform (Walsh Transform)¹ is a method widely used in signal processing as an alternative to the discrete Fourier transform in spectral analysis. FFT uses double-precision values, whereas Walsh decomposes a signal into square waves of different oscillation frequencies, so all coefficients are integers with the same absolute value. This avoids floating-point multiplication and improves the computation speed.

Thus, FWT and FFT have the same core idea: transforming an array. Let the result of applying the fast Walsh transform to array $A$ be $FWT[A]$.

The core idea of FWT is as follows:

We need a new sequence $C$, obtained from sequences $A$ and $B$ by some operation, namely $C = A \cdot B$.

First compute $FWT[A], FWT[B]$ by the forward transform. Then, from $FWT[C]=FWT[A] \cdot FWT[B]$, compute $O(n)$ in $FWT[C]$ time, where $\cdot$ denotes multiplication at corresponding positions of the sequences.

Finally, apply the inverse operation to obtain the original sequence $C$. The time complexity is $O(n \log{n})$.

In algorithm competitions, FWT is used to solve convolutions whose indices are combined by bitwise operations.

The formula is $C_{i} = \sum_{i=j \oplus k}A_{j} B_{k}$.

Here $\oplus$ is one of the binary bitwise operations.

We use $\cup$ (bitwise OR), $\cap$ (bitwise AND), and $\oplus$ (bitwise XOR) as examples.

FWT Operations

OR Convolution

If $k=i\cup j$, then the positions where the binary digits of $i$ are $1$ and the positions where the binary digits of $j$ are $1$ must form a subset of the positions where the binary digits of $k$ are $1$.

To obtain $FWT[C] = FWT[A] \cdot FWT[B]$, we must construct the rule for this FWT.

By definition, we can construct $FWT[A]_i = A'_i = \sum_{i=i\cup j}A_{j}$ to represent the indices $j$ whose positions with bit value $1$ form a subset of the positions with bit value one in $i$.

Therefore:

Now let us see how to compute $FWT[A]$.

We clearly cannot enumerate all pairs, which would take $O(n^2)$. Since global enumeration is too slow, we use divide and conquer.

Split the whole interval into two halves. After this split, the binary representation of an index follows a regular pattern.

Let $A_0$ denote the first half of $A$, and $A_1$ the second half. Then $A_0$ has the highest bit of its largest index equal to $0$, so its subsets are just subsets within itself because the highest bit is $0$. In contrast, $A_1$ has highest bit $1$; its valid subsets include not only itself but also subsets whose highest bit is $0$. Therefore,

Here merge concatenates the two arrays like strings, $+$ is ordinary addition, and it adds corresponding binary positions.

The divide-and-conquer construction completes a merge in $O(\log{n})$ time. Since each merge performs one operation, it obtains $O(n\log{n})$ in $FWT[A]$ time.

The inverse is also simple. The subsets of $A_0$ are itself ($A_0 = FWT[A_0]$), while the subsets of $A_1$ are $FWT[A_0] + FWT[A_1]$; hence the inverse recurrence is

Here is the implementation. The forward and inverse transforms can share one function: use $\text{type}=1$ for the forward transform and $\text{type}=-1$ for the inverse transform.

void Or(ll *a, ll type) {  // Iterative implementation with a smaller constant factor
  for (ll x = 2; x <= n; x <<= 1) {
    ll k = x >> 1;
    for (ll i = 0; i < n; i += x) {
      for (ll j = 0; j < k; j++) {
        (a[i + j + k] += a[i + j] * type) %= P;
      }
    }
  }
}

AND Convolution

The same reasoning as for OR gives analogous conclusions.

Here is the implementation. Use $\text{type}=1$ for the forward transform and $\text{type}=-1$ for the inverse transform.

void And(ll *a, ll type) {
  for (ll x = 2; x <= n; x <<= 1) {
    ll k = x >> 1;
    for (ll i = 0; i < n; i += x) {
      for (ll j = 0; j < k; j++) {
        (a[i + j] += a[i + j + k] * type) %= P;
      }
    }
  }
}

XOR Convolution

XOR convolution is based on the following principle.

Let $x\circ y$ denote the parity of the number of 1s in $x\cap y$, namely $1$. The definition is $x\circ y=\text{popcnt}(x\cap y)\bmod 2$, and it is easy to see that $(x\circ y)\oplus (x\circ z)=x\circ(y\oplus z)$.

The operation for $FWT[A]$ follows in the same way.

Set $FWT[A]_i=\sum_{i\circ j=0}A_j-\sum_{i\circ j=1}A_j$. We prove that $FWT[C] = FWT[A] \cdot FWT[B]$.

Let us see how to compute the values of $A,B$ quickly; again, use divide and conquer.

For the subsequence $FWT[A_0]$ whose current bit of $i$ is $0$, the $\circ$ operation gives the same result when combined with $0$ or $1$, because $0\cap 0=0,0\cap1=0$. Thus, in $FWT[A]=\sum_{i\circ j=0}A_j-\sum_{i\circ j=1}A_j$, the term $\sum_{i\circ j=1}A_j=0$ is zero.

For the subsequence $A_1$ whose current bit of $i$ is $1$, the $\circ$ operation gives $0$ when combined with $0$, and $1$ when combined with $1$, because $1\cap 0=0,1\cap1=1$.

Combining these observations gives:

That is:

The inverse transform follows immediately:

The code is given below. Use $\text{type}=1$ for the forward transform and $\text{type}=\frac{1}{2}$ for the inverse transform.

void Xor(ll *a, ll type) {
  for (ll x = 2; x <= n; x <<= 1) {
    ll k = x >> 1;
    for (ll i = 0; i < n; i += x) {
      for (ll j = 0; j < k; j++) {
        (a[i + j] += a[i + j + k]) %= P;
        (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;
        (a[i + j] *= type) %= P;
        (a[i + j + k] *= type) %= P;
      }
    }
  }
}

XNOR Convolution

By analogy with XOR, we obtain the formula:

$FWT[A]_{i} = \sum_{C_1}A_{j} - \sum_{C_2}A_{j}$ ($C_1$ denotes the indices for which $\text{popcnt}(x\cup y)\bmod 2$ is $0$, while $C_2$ denotes the indices for which $\text{popcnt}(x\cup y)\bmod 2$ is $1$).

Another View of FWT

Let $c(i,j)$ be the contribution coefficient of $A_j$ to $FWT[A]_i$. We can describe the FWT process as follows:

Since:

A direct proof gives $c(i,j)c(i,k)=c(i,j\odot k)$, where $\odot$ is any bitwise operation.

The function $c$ also has an important property: it can be processed bit by bit.

For example, during the transform we have:

This direct expression is inefficient, so split it as follows:

Comparing the preceding and following expressions, the only difference in $i,j$ is the highest bit.

Remove the highest bit of $i,j$ and call the result $i',j'$; let $i_0$ be the highest bit of $i$. Then:

If $i_0=0$, then:

If $i_0=1$, then:

In other words, we only need:

These four numbers suffice to complete the transform. We call this matrix the bit matrix.

For the inverse transform, we need the inverse of the bit matrix above.

If the inverse matrix is $c^{-1}$, a similar argument gives the original values:

The inverse matrix does not always exist. For example, if there is a row $0$ or a column $0$, the matrix has no inverse, so the construction requires care.

Bitwise OR

We can construct the matrix

It satisfies $c(i,j)c(i,k)=c(i,j\cup k)$. This is exactly the same as the previously derived $FWT[A]=\text{merge}(FWT[A_0], FWT[A_0]+FWT[A_1])$. Similarly, the following matrix also satisfies the condition, but we generally use the first one:

Although the following matrix also satisfies $c(i,j)c(i,k)=c(i,j\cup k)$, it contains a row $0$, has no inverse, and is therefore invalid:

For the inverse transform, invert the matrix. For the top matrix, the result is:

Then substitute this inverse matrix into the forward-transform procedure.

Bitwise AND

We can construct the matrix

It satisfies $c(i,j)c(i,k)=c(i,j\cap k)$.

The inverse matrix is:

Bitwise XOR

We can construct the matrix

It satisfies $c(i,j)c(i,k)=c(i,j\oplus k)$.

The inverse matrix is:

FWT Is a Linear Transform

FWT is a linear transform. In other words, it satisfies:

and:

K-Dimensional FWT

At its core, a bitwise operation acts on an $n$-dimensional $\{0,1\}$ vector. OR takes the $\max$ in each dimension, AND takes the $\min$, and XOR adds corresponding coordinates modulo $\bmod 2$.

A key feature of bitwise operations is that the coordinates of the vector are independent.

Extend $\{0,1\}$ to $[0,K)\cap \mathbf{Z}$, that is, to base-$K$, and see what happens.

Max Operation

Extend the $\cup$ operation to base-$K$. Define $i\cup j$ to take the $\max$ in each digit; then

If $j=k$, the equation becomes:

Thus, in every row, the $1$s must occur before the $0$s; placing them afterward is invalid. One valid construction is:

Its inverse is:

Min Operation

Extend the $\cap$ operation to base-$K$. Define $i\cap j$ to take the $\min$ in each digit; then

If $j=k$, the equation again becomes:

Thus, in every row, the $1$s must occur after the $0$s; placing them beforehand is invalid. One valid construction is:

Its inverse is:

The first two operations are used less often. The more common one is:

Carryless Addition

Extend the $\oplus$ operation to base-$K$. Define $i\oplus j$ as digitwise addition modulo $\bmod K$; then

Choosing $c(i,j)=\omega_{K}^j$ satisfies the requirement:

However, identical rows give a singular matrix, so instead choose $c(i,j)=\omega_{K}^{(i-1)j}$.

This gives the following matrix:

This is a Vandermonde matrix; its inverse is:

If the modulus in the problem contains a suitable root of unity, the transform is straightforward to implement.

However, a root of unity may not exist modulo the given modulus. Consider an extension field: define $x$ with $x^K=1$, and substitute $x$ into the computation. Every value then becomes a polynomial in $x$ of degree $k-1$. We only need to compute modulo $\bmod {x^K-1}$. The matrix can therefore be written as:

This may introduce zero divisors, meaning one value has multiple representations, so its actual value is ambiguous.

Instead of reducing $\bmod {x^K-1}$, take $\bmod$ modulo the cyclotomic polynomial $\Phi_{K}(x)$. It makes $x$ have order $k$ and is irreducible over $\mathbb{Q}$, so the above computation can be performed modulo $\bmod {\Phi_{K}(x)}$.

One issue is that the constants in $\bmod \Phi_{K}(x)$ are large because $\Phi$ itself is a polynomial. Since $\Phi_{K}(x)\mid x^k-1$, we can reduce modulo $\bmod x^k -1$ during the computation and reduce modulo $\bmod \Phi_{K}(x)$ only at the end.

Examples

References

• Taojiang's algorithm notes (Taojiang's algorithm notes)
• ZnPdCo's blog (ZnPdCo's blog)