For a final project in my class, I decided to try to simulate a quantum computer and implement Grover's algorithm. I followed this excellently written blog post by Craig Gidney, and was successful in getting Grover's algorithm to work using matrix multiplication.
However, I wanted to figure out how to generate the quantum circuit that implements the oracle. The blog post by Craig Gidney intentionally skips over how this is done, and instead says that an oracle can always be represented as a matrix such as
$$ \pmatrix{ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & -1 & 0\\ 0 & 0 & 0 & 1 } $$
In this case, Grover's algorithm would be searching over a "database" of 4 items, and since the -1 is in the 3rd row, the entry we are searching for would be 3.
Using this matrix as the oracle, my implementation of Grover's algorithm works, with the correct "entry" coming up as the result with a >98% probability for all oracles I've tried.
The blog post links to this mathoverflow answer by Robin Kothari on how to actually go about constructing the oracle. Using it as a guide, I was also able to write code that finds the appropriate quantum circuit to simulate a classical circuit (which represents the search function $f$), and generates the matrix which performs
$$ \left | x \right \rangle \left | 0^{k} \right \rangle ↦ \left | f(x) \right \rangle \left | \mathrm{junk}(x) \right \rangle $$
From here, using the "uncomputation trick," it's easy to get
$$ \left | x \right \rangle \left | 0^{k} \right \rangle \left | b \right \rangle ↦ \left | x \right \rangle \left | 0^{k} \right \rangle \left | b \oplus f(x) \right \rangle $$
The resulting matrix should correspond to the oracle function.
According to the Robin Kothari's answer, my next step would be to set $\left | b \right \rangle$ to $\frac{1}{\sqrt{2}} \left ( \left | 0 \right \rangle - \left | 1 \right \rangle \right )$. As far as I understand, this is equivalent to setting $\left | b \right \rangle$ to $\left | 1 \right \rangle$, and then applying a Hadamard gate. This is because
$$ H (\left | 1 \right \rangle) = \pmatrix{0 & 1} \frac{1}{\sqrt{2}} \pmatrix{1 & 1 \\ 1 & -1}=\frac{1}{\sqrt{2}} \pmatrix{1 & -1} $$
However, when I use the matrix I generated and apply the necessary Hadamard gate, Grover's function does not work as before.
For example, let's say $f$ can be represented with the boolean expression $(x_{1} \land x_{2})$. This is a valid circuit because there is only one solution. My code simulates this with a single Toffoli gate which acts on a third qubit. Then a CNOT gate copies the state of the third qubit onto a fourth. Then the Toffoli gate is reversed. This is represented by the matrix
$$ \pmatrix{ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\ } $$
As far as I can tell, this matrix does the trick. It maps $\left | 0000 \right \rangle ↦ \left | 0000 \right \rangle$, $\left | 0100 \right \rangle ↦ \left | 0100 \right \rangle$, and $\left | 1100 \right \rangle ↦ \left | 1101 \right \rangle$.
However, when I use this matrix as the oracle, the vector representing the state of the qubits comes out to be
$$ \frac{1}{2^{\frac{3}{2}}}\pmatrix{-1 & 0 & 0 & -1 & -1 & 0 & 0 & -1 & -1 & 0 & 0 & -1 & 0 & -1 & -1 & 0} $$
Of course, when measured, this will not result in a high likelihood of the desired answer of $\left | 3 \right \rangle$.
What's the reason that my oracle matrix is in a different form than the one described in Craig Gidney's blog post? Did I miss a step in the process for generating it, or should they be different? If they should be different, why should they both work, and what did I miss?
