# Is there a fast sparse Hadamard transform?

Suppose I give you an $$n$$-qubit state vector as a classical list of numbers (or as an oracle that can query the amplitudes). I tell you this state vector will contain exactly $$k$$ non-zero amplitudes, after you apply a Hadamard transform to it. You could determine what those non-zero amplitudes are in $$O(n2^n)$$ time by applying the fast Walsh-Hadamard transform. But that seems quite inefficient, given the promise that only $$k$$ of the $$2^n$$ outputs are actually needed.

Is there a way to compute a sparse representation of the the output, in time that scales polynomially in $$k$$ and $$n$$, like $$O(k^9 n^9)$$?

This would be the Hadamard analogue of the sparse fast Fourier transform. Or more specifically the high-dimension all-dimension-lengths-equal-to-2 sparse Fourier transform.

For example, in the $$k=1$$ case, you can just look at the amplitudes of $$|1\rangle$$, $$|2\rangle$$, $$|4\rangle$$, $$|8\rangle$$, ..., $$|2^{n-1}\rangle$$ to check if they are equal or opposite to the amplitude of $$|0\rangle$$. The equal-or-opposite-ness directly reads off the bits of the state that has the single non-zero amplitude in the output. Under the hood this is because $$HX = ZH$$: we're checking for the Z's negation to infer whether each qubit will be bit flipped away from 0 or not after the Hadamard. This takes time $$O(n)$$.

• I'm not sure I follow your argument for why the $k=1$ case has $O(\text{lg}\ n)$ time. Isn't it the case that you're looking at $n$ amplitudes, so its time is $O(n)$? Commented Feb 9 at 8:50
• @DaftWullie before Hadamarding you use the oracle to query $|0\rangle$,$|1\rangle$, $|2\rangle$, $|4\rangle$, etc. so as to get the bitstring of the basis state that, after Hadamarding, will be non-zero. The promise is on the Hadamarded state having $k$ nonzero entries, the queries are counted before Hadamarding. Commented Feb 9 at 12:14
• @MarkSpinelli Yes, but the bit string you're reading is length $n$ isn't it? Commented Feb 9 at 12:16
• @DaftWullie Whoops, you're right. I mixed up whether $n$ was the exponent or the number of amplitudes. Fixed. Commented Feb 9 at 16:35
• @MarkSpinelli For the $k=2$ equal-amplitude case, half of the input amplitudes will be zero. So no it won't just give you $x_0 \oplus x_1$ without at least adding what to do when you get a zero. Commented Feb 9 at 16:38

This partial answer provides a simple algorithm that generalizes the $$k=1$$ case and exploits sparsity promise to beat the Walsh-Hadamard transform when $$k$$ is constant. Specifically, the algorithm makes $$O(n^k)$$ calls to the oracle and runs in time $$O(n^k)$$.
The Hadamard transform of a $$k$$-sparse quantum state $$|\psi\rangle=\frac{1}{\sqrt k}\sum_{i=1}^k\alpha_i|b_i\rangle$$ where $$\alpha_i\in\mathbb{C}\setminus\{0\}$$ and $$b_i\in\{0,1\}^n$$ is $$H^{\otimes n}|\psi\rangle=\frac{1}{\sqrt{k2^n}}\sum_{y\in\{0,1\}^n}\left(\sum_{i=1}^k\alpha_i(-1)^{b_i\cdot y}\right)|y\rangle\tag1$$ where $$\cdot$$ denotes the dot product of vectors in $$\mathbb{Z}_2^n$$. Defining $$x_{i,j}=(-1)^{b_{i,j}}$$ we can rewrite $$(1)$$ as $$H^{\otimes n}|\psi\rangle=\frac{1}{\sqrt{k2^n}}\sum_{y\in\{0,1\}^n}\left(\sum_{i=1}^k\alpha_i\prod_{j=1}^nx_{i,j}^{y_j}\right)|y\rangle\tag2$$ where $$y_j$$ denotes the $$j$$th bit of the bit string $$y$$. Note that the amplitude in front of $$|y\rangle$$ is a polynomial in $$x_{i,j}$$ whose degree is the Hamming weight of $$y$$.
The algorithm proceeds in three steps. First, we query the oracle for the $$O(n^k)$$ amplitudes $$\beta_y:=\langle y|H^{\otimes n}|\psi\rangle$$ of all $$|y\rangle$$ whose Hamming weight is at most $$k$$. Next, we construct the system of polynomial equations \begin{align} \sum_{i=1}^k\alpha_i&=\beta_{0\dots 0}\tag3\\ \sum_{i=1}^k\alpha_ix_{i,u}&=\beta_{0\dots 010\dots 0}\tag4\\ \sum_{i=1}^k\alpha_ix_{i,u}x_{i,v}&=\beta_{0\dots 010\dots 010\dots 0}\tag5\\ \dots \end{align} where $$(3)$$ represents a single equation corresponding to the zero bit string, $$(4)$$ represents $$n$$ equations corresponding to $$y$$ with Hamming weight one, $$(5)$$ represents $$\frac{n(n-1)}{2}$$ equations corresponding to $$y$$ with Hamming weight two, and so on until $$y$$ with Hamming weight $$k$$. Note that equations $$(3)$$ and $$(4)$$ correspond to the algorithm for $$k=1$$ described in the question. Finally, in the third step of the algorithm, we solve the system of equations for $$\alpha_i$$ and $$x_{i,j}$$. This can be done for example by introducing new variables for products of $$x_{i,j}$$ in order to reduce the system of $$O(n^k)$$ polynomial equations in $$O(nk)$$ variables to a system of $$O(n^k)$$ linear equations in $$O(n^k)$$ variables. The resulting algorithm will run in time polynomial in $$n$$ and exponential in $$k$$.