Abstract
The definition of the propagator $\Delta(x)$ in the path integral formalism (PI) is different from the definition in the operator formalism (OF). In general the definitions agree, but it is easy to write down theories where they do not. In those cases, are the PI and OF actually inequivalent, or is it reasonable to expect that the $S$ matrices of both theories agree?
For definiteness I'll consider a real scalar field $\phi$, with action $$ S_0=\int\mathrm dx\ \frac{1}{2}(\partial \phi)^2-\frac{1}{2}m^2\phi^2 \tag{1} $$
Path integral formalism
In PI we insert a source term into the action, $$ S_J=S_0+\int\mathrm dx\ \phi(x)\ J(x) \tag{2} $$
The mixed term $\phi J$ can be simplified with the usual trick: by a suitable change of variables, the action can be written as two independent terms
$$ S_J=S_0+\int \mathrm dx\;\mathrm dy\ J(x)\Delta_\mathrm{PI}(x-y)J(y) \tag{3} $$ and this relation defines $\Delta_\mathrm{PI}(x)$: to get the action into this form we have to solve $(\partial^2+m^2)\Delta_{\mathrm{PI}}=\delta(x)$, i.e., in PI the propagator is defined as the Green function of the Euler-Lagrange equations of $S_0$. This definition is motivated by the the fact that when $S_J$ is written as $(3)$ the partition function can be factored as $$ Z[J]=Z[0]\exp\left[i\int \mathrm dx\;\mathrm dy\ J(x)\Delta_\mathrm{PI}(x-y)J(y) \right] \tag{4} $$ which makes functional derivatives trivial to compute. For example, if we differente $Z[J]$ two times we get $$ \langle0|\mathrm T\ \phi(x)\phi(y)|0\rangle=\Delta_\mathrm{PI}(x-y) \tag{5} $$
In this formalism, the propagator is always the Green function of the differential operator of the theory.
Operator formalism
In OF the propagator is defined as the contraction of two fields: $$ \Delta_\mathrm{OF}(x-y)\equiv \overline{\phi(x)\phi}(y)\equiv \begin{cases}[\phi^+(x),\phi^-(y)]&x^0>y^0\\ [\phi^+(y),\phi^-(x)]&x^0
In general, $\Delta_\mathrm{OF}$ is an operator, but if it commutes with everything (or, more precisely, if the propagator is in the centre of the operator algebra) we can prove Wick's theorem, which in turns means that $$ \langle 0|\mathrm T\ \phi(x)\phi(y)|0\rangle=\Delta_\mathrm{OF}(x-y) \tag{7} $$ i.e., the propagator coincides with the two-point function. This makes it very easy to see, for example, that $$ \Delta_\mathrm{OF}=\Delta_\mathrm{PI} \tag{8} $$
In this theory, the fact that the propagator is a Green function is a corollary and not a definition. The theorem may fail if the assumptions are not satisfied.
The discrepancy
The positive/negative frequency parts of $\phi$ are the creation and annihilation operators, which in OF usually satisfy $$ [\phi^+(x),\phi^-(y)]\propto\delta(x-y)\cdot1_\mathcal H \tag{9} $$ and therefore $\overline{\phi(x)\phi}(y)$ is a c-number. This means that the assumptions of Wick's theorem are satisfied and $(8)$ holds.
The relation $(9)$ can be derived from one of the basic assumptions of OF: the canonical commutation relations: $$[\phi(x),\pi(y)]=\delta(x-y)\cdot1_\mathcal H\tag{10} $$
But if we use any non-trivial operator in the r.h.s. of $(10)$ instead of a $1_\mathcal H$, Wick's theorem is violated and in general $\Delta_\mathrm{PI}\neq \Delta_\mathrm{OF}$. One could argue that the r.h.s. of $(10)$ is fixed by the Dirac prescription $\{\cdot,\cdot\}_\mathrm{D}\to\frac{1}{i\hbar}[\cdot,\cdot]$, where $\{\cdot,\cdot\}_\mathrm{D}$ is the Dirac bracket. In the Standard Model, it's easy to prove that $\{\cdot,\cdot\}_\mathrm{D}$ is always proportional to the identity, but in a more general theory we may have complex constraints which would make the Dirac bracket non-trivial - read, not proportional to the identity - and therefore $\Delta_\mathrm{OF}\neq\Delta_\mathrm{PI}$.
Scalar, spinor and vector QFT's always satisfy a relation similar to $(10)$, and therefore OF and PI formalisms agree. But in principle it is possible to study more general QFT's where we use a commutation relation more complex than $(10)$. I don't know of any practical use of this, but to me it seems to me that we can have a perfectly consistent theory where PI and OF formalisms predict different results. Is this correct? I hope someone can shed any light on this.
EDIT
I think it can be useful to add some details to what I said about Dirac brackets, defined as $$ \{a,b\}_\mathrm{DB}=\{a,b\}_\mathrm{PB}-\{a,\ell^i\}_\mathrm{PB} M_{ij}\{\ell^j,b\}_\mathrm{PB}\tag{11} $$ where $\ell^i$ are the constraints and $M^{ij}=\{\ell^i,\ell^j\}_\mathrm{PB}$. As $\{q,p\}_\mathrm{PB}=\delta(x-y)$, the only way to get non-trivial Dirac brackets is through the second term. This may happen if we have non-linear constrains such that the second term is a function of $p,q$. If in some case we have non-linear constrains the matrix $M$ will depend on $p,q$ and $\{p,q\}_\mathrm{DB}$ will be a function of $p,q$. If we translate this into operators, we shall find $$ [\pi,\phi]=\delta(x-y)\cdot1_\mathcal H+f(\pi,\phi)\tag{12} $$ and so $[\pi,\phi]$ won't commute with neither $\pi$ nor $\phi$ as required by Wick's theorem. (This term $f$ is perhaps related to the $\hbar^2$ terms in QMechanic's answer, and the higher order terms in, e.g., the Moyal bracket).
General comments to the question (v1):
Any textbook derivation of the correspondence between $$\tag{1} \text{Operator formalism}\qquad \longleftrightarrow \qquad \text{Path integral formalism}$$ is just a formal derivation, which discards contributions in the process, cf. e.g. this Phys.SE post.
Rather than claiming complete understanding and existence of the correspondence (1), it is probably more fair to say that we have a long list of theories (such as e.g. Yang-mills, Cherns-Simons, etc.), where both sides of the correspondence (1) have been worked out.
The correspondence (1) is mired with subtleties. Example: Consider a non-relativistic point particle on a curved target manifold $(M,g)$ with classical Hamiltonian $$\tag{2} H_{\rm cl} ~=~\frac{1}{2} p_i p_jg^{ij}(x), $$ which we use in the Hamiltonian action of the phase space path integral. Then one may show that the corresponding Hamiltonian operator is $$\tag{3}\hat{H}~=~ \frac{1}{2\sqrt[4]{g}} \hat{p}_i\sqrt{g}~ g^{ij} ~\hat{p}_j\frac{1}{\sqrt[4]{g}}+ \frac{\hbar^2R}{8} +{\cal O}(\hbar^3),$$ cf. Refs. 1 & 2. The first term in eq. (3) is the naive guess, cf. my Phys.SE answer here. The two-loop correction proportional to the scalar curvature $R$ is a surprise, which foretells that a full understanding of the correspondence (1) is going to be complicated.
References:
F. Bastianelli and P. van Nieuwenhuizen, Path Integrals and Anomalies in Curved Space, 2006.
B. DeWitt, Supermanifolds, Cambridge Univ. Press, 1992.
