Why not implement quantum circuits on classical computers? 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?
 A: Matrix multiplication is polynomial in the number of matrix elements. In quantum computing the number of matrix elements is basically the number of elements in the quantum state vector, which is exponential in the size of the input (the number of qubits).
A: The false point in the premise is that a composition of single-qubit and two-qubit gates will be represented by some sort of composition of 2×2 and 4×4 matrices. In fact, if a single-qubit gate (for two-qubit this is similar) acts on the $i$-th qubit out of $n$, the matrix must be cast to
$$(Id_2)^{\otimes (i-1)} \otimes A \otimes (Id_2)^{\otimes (n-i)},$$
where $A$ is the original $2×2$ matrix, in order to describe its action on the $n$-qubit state. This is a matrix of dimensions $2^n × 2^n$ regardless of the type of the gate.
While the tensor product with identities looks trivial, let's just recall what it means for an application on a state vector using the most straightforward algorithm:
for each index u, 0 ≤ u < 2^n:
  define u_i = bit value of i-th bit of u
  define v = u XOR 2^i
  if u_i = 0:
    new (psi[u], psi[v]) = a*psi[u] + b*psi[v]
  else:
    new (psi[u], psi[v]) = d*psi[u] + c*psi[v]

Yes, in the core of this pseudo-algorithm there's just a 2-dimensional multiplication. But it is done $2^n$ times.
Even before one gets to perform the actual classical simulation, one can run out of available memory in even storing the state at any particular instant. Take $n = 100$ qubits, that's $2^{100}$ complex double-precision amplitudes (128 bits each). I'm getting an estimate exceeding the data storage capacity of all computers on Earth by some 14 orders of magnitude so we're not going to see that anytime soon. Meanwhile, a quantum computer with 100 bits would just about start to be interesting for applications, in theory it's not uncommon to see many more than that envisaged.
A: 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.
A: 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.
