Question: I have somehow lost sight of the big picture in my study of quantum computing. I understand that we still don't know if quantum computers are more powerful than classical ones, in the sense of computational complexity. But I'm wondering: 

Since every quantum gate can be seen as a (unitary) matrix multiplication, and matrix multiplication can be done classically in polynomial time, where is the supposed advantage coming from? What prevents me from taking a polynomial size quantum circuit and implementing it as a polynomial number of matrix multiplications (that require polynomial time)? what is the operation that is classically unavailable?

Answer: 
To put slightly differently the Hilbert Space (H) for your quantum circuit is growing exponentially ($ \text{dim H} = 2^n$) where n is the number of qubits. There are a class of quantum circuits that can be implemented in polynomial time but they require that the gates in the quantum circuit are restricted to the Pauli group and/or the normalizer of the Pauli group(Gottesman Knill Theorem). This also restricts what kind of error models one employs in the quantum circuit. The basic reason is that there is an isomorphism one can define between the Pauli group and vector space over the finite field with elements. Multiplication of the Pauli group elements becomes a symplectic product of the vector space elements.
Because the matrices are bigger than you think.

Consider some simple circuit on 8 qubits; perhaps a Fourier transform. It uses ~40 gates, and each gate multiplies the 8-qubit system state vector by some matrix. To simulate this system classically. you might expect us to have to invest something like $40\cdot 8=320$ arbitrary units of computational work (AUCW).

But the system state vector has an amplitude for every possible classical state of the system. Eight bits can be assigned $2^8 = 256$ possible values, so our state vector must have 256 entries not eight or sixteen. Correspondingly, a matrix that transforms this system must have size $256 \times 256$, though just $256$ is really a much better rough estimate of the "work size" since a single gate's matrix will be sparse. Our guess of $40 \cdot 8 = 320$ AUCW was way off! A much better guess is $40 \cdot 2^8 = 40 \cdot 256 = 10240$. So $10$ kiloAUCW.

This problem gets way *way* worse as the number of qubits gets larger. With 50 qubits the state vector has $2^{50}$ entries, the matrices are $2^{50} \times 2^{50}$, the Fourier transform uses ~1250 gates, and we have to spend a million *million* **million** AUCWs to do the classical simulation. Whereas the quantum computer does a mere 1250 AUQW.