Is there a fast sparse Hadamard transform?

Question

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)$.

Answer 1

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$.