Apart from the formal result about #P-hardness, there's something worth touching on, about the nature of strong simulation itself. I'll comment first on strong simulation, and then specifically on the quantum case.
1. Strong simulation even of classical randomised computation is hard
Strong simulation is a very powerful concept — not only in the fact that it is a useful concept to consider, but in the more practical sense that it would allow you to do very powerful computations, even if you could do it in a purely classical setting.
At issue here is that a process with random outcomes — whether or not it is quantum — does not automatically come equipped with a way to compute the probabilities of its events; not even for the probabilities of the events which you actually see realised. The power of randomised computation, so to speak, is that you don't have to worry much about what the precise probability of an event is: it suffices to sample and attempt to realise that event sufficiently often to be confident of roughly what that probability is (or, for instance, just that the probability is non-zero).
Asking for actual probabilities is a very, well, strong requirement to ask of a simulation, in the following sense:
Proposition. Strongly simulating a polynomial-time randomized computation with zero error is #P-hard.
Proof: For any non-deterministic Turing machine N, consider the question of how many branches it accepts on for an input $x \in \{0,1\}^m$. It is enough for us to consider the case where N makes a non-deterministic choice at every transition, and runs for time $m \in \mathrm{poly}(n)$ in every branch, so that the branches of the computation can be indexed by the strings $z \in \{0,1\}^m$. The computation performed by N can be represented similarly as a deterministic computation (a function $f$) depending on the input $x$ and a given branch string $z$. If we represent the status of 'accept' vs. 'reject' by a bit $a \in \{0,1\}$ which is computed as $a = f(x,z)$.
Determining the number of branches $z \in \{0,1\}^m$, for which $f(x,z) = a$ for a given $x \in \{0,1\}^n$, is essentially the canonical #P-complete problem, using this connection to non-deterministic Turing machines.
Consider a polynomial-time randomised computation. We can describe this computation as a deterministic classical computation, performed with the help of a uniformly random bit-string of length $m$.
Suppose that we are interested in whether a particular bit yields the outcome '1'.
For an input $x \in \{0,1\}^n$ and random bit-string $z \in \{0,1\}^m$, let $f(x,z) \in \{0,1\}$ be the value that this bit takes: then $f(x,z)$ can be computed in polynomial time.
The probability $P(a)$ that this computation gives the result $a \in \{0,1\}$ is then
$$ P(a) = \frac{\# \bigl\{ z \in \{0,1\}^m \;\big|\; f(x,z) = a \bigr\}}{2^m}\,. $$
Because we could choose $f$ to be the function determining the acceptance condition of a non-deterministic Turing machine on input $x$ in branch $z$, then it is #P-complete to compute $P(a)$ exactly.
Proposition. Strongly simulating a polynomial-time randomized computation, with any relative error, is NP-hard.
- Proof. An important corner-case of the problem of approximating a probability with relative error, is the case where the exact probability is equal to zero. In this case, any process which gives the correct probability up to any multiplicative factor must correctly produce the exact probability, if that probability happens to be 0. Similarly, any process which approximates an event with positive probability up to a positive scalar factor, must yield a probability estimate which is greater than zero. That process can then be used to determine whether the number of accepting branches of a non-deterministic Turing machine is zero or non-zero.
From these two observations, you should take away the idea that strong simulation is a strong requirement — in many cases, unfairly strong — to make of a simulation method: it allows you to do much more powerful things than the computational model itself might be capable of.
2. Strong simulation of quantum computation is very hard
One difference between classical and quantum computation is on how difficult we think it is (in a complexity-theoretic sense) to strongly simulate them.
We know that it is NP-hard to strongly simulate a polynomial-time classical randomised process to relative error less than 1. Simulating a quantum process isn't going to be any easier. However, there is reason to believe that even with unreasonably powerful computational resources which would allow us to strongly simulate classical processes with bounded relative error, it will still be difficult to strongly simulate quantum processes.
Theorem (Stockmeyer [1]). For any counting problem in #P, and any constant $d \geqslant 0$ constant, the problem of computing the counting problem within to a $(1 + O(n^{-d}))$ factor is in $\mathbf{FP^{\:\!NP^{\:\!NP}}}$.
The difficulty in the q quantum case is that while classical probabilities correspond to counting the number of accepting branches of a non-deterministic Turing machine, quantum computation corresponds more closely to counting the difference between the number of accepting and rejecting branches of a non-deterministic Turing machine. That is, where a classical probability corresponds to a #P function, a quantum amplitude corresponds to a GapP function, which is the set of functions which may be expressed as a difference between two #P functions. The connection between quantum computation and GapP may be made formal, but on a high level, it is essentially because quantum computation can involve destructive interference between amplitudes associated with different events. More formally:
Proposition (Liberal paraphrase of Theorem 3.2 of [2]). For any $f \in \mathbf{GapP}$, there is a polynomial-time quantum algorithm $Q$, with an accepting configuration $c$, and a polynomial $p$ such that, for all $x \in \{0,1\}^n$, $\langle c | Q | x \rangle = -f(x) \cdot 2^{-p(n)/2}$.
Approximating GapP functions even to constant factors is difficult, because computing a GapP function to within even a constant factor determines whether it is positive or negative. If you can do this, you can immediately solve the PP-complete problem of determining whether or not the number of accepting paths is greater than the number of rejecting paths; and if you do it repeatedly, with a number of related quantum computations in which you artificially inflate the number of accepting or rejecting paths by doing your quantum computation conditionally, you can compute the exact GapP function essentially by binary search (the same way you can find a satisfying solution to a boolean formula, if one exists, given an oracle which simply tells you whether a solution exists).
References
In addition to the references mentioned above, Maarten Van den Nest's article [3] mentioned by Martin Schwarz is noteworthy for presenting the first definition of 'strong simulation' of quantum systems (to distinguish it from the more reasonable standard of weak simulation), and also for presenting a number of ideas of the links between classical and quantum computation in the context of simulation.
Stockmeyer.
The complexity of approximate counting.
Proceedings of STOC '83 (pp.&thinsp.118-126), 1983.
[PDF available at acm.org]
Fenner, Green, Homer, and Pruim.
Determining Acceptance Possibility for a Quantum Computation is Hard
for the Polynomial Hierarchy.
Proceedings of the Royal Society London A vol. 455 (pp. 3953–3966), 1999.
[arXiv:quant-ph/9812056]
Van den Nest.
Classical simulation of quantum computation, the Gottesman-Knill theorem, and slightly beyond.
Quantum Information and Computation, vol. 10 (pp. 258-271), 2010.
[arXiv:0811.0898].