They are indeed equivalent. It is quite illuminating and satisfying to see the proof of this. Consider a general action of the form
$$ Z = \int \exp \left( -S_0[\varphi]+\frac{\lambda}{2} \;\mathcal O[\varphi] \mathcal O[\varphi] \right) \; \mathrm D[\varphi] $$
where $S_0$ is the action which we perturb around (we will not presume this to be a free theory).
1. Mean field theory. The essential assumption is that
$$\mathcal O([\varphi]) = M + \underbrace{(\mathcal O([\varphi]) -M)}_{= \textrm{ small}}.$$
Squaring this gives us $\mathcal O[\varphi] \mathcal O[\varphi] \approx 2M \mathcal O[\varphi] - M^2$. Plugging this into the partition function gives us the mean field approximation:
$$ \boxed{ Z_\textrm{mf}[M] = e^{- \lambda M^2 / 2} \int \exp \left( -S_0[\varphi]+\lambda M\mathcal O[\varphi] \right) \; \mathrm D[\varphi] }$$
with the self-consistency relation
$$ M = \langle \mathcal O[\varphi] \rangle_\textrm{mf} = \frac{1}{\lambda} \partial_M \ln \left( e^{\lambda M^2/2} Z_\textrm{mf}[M] \right). $$
Note that the latter is equivalent to $\boxed{ \partial_M \ln Z_\textrm{mf}[M] = 0 }$ (which one can interpret as saying that the mean field extremizes the free energy).
2. Hubbard-Stratonovich transformation. The essential identity is
$$ \exp\left( \frac{\lambda}{2} \; \mathcal O[\varphi] \mathcal O[\varphi] \right) = \int \exp \left( - \frac{\lambda \mathcal M^2}{2} + \lambda \mathcal M \mathcal O[\varphi] \right) \; \mathrm D[\mathcal M]. $$
Plugging this into the original partition function, we get
$$ \boxed{ Z = \int \mathrm D[\mathcal M] \; e^{-\lambda \mathcal M^2/2} \int \mathrm D[\varphi] e^{-S_0[\varphi] + \lambda \mathcal M \mathcal O[\varphi]} = \int e^{-S_\textrm{HS}[\mathcal M]} \mathrm D [\mathcal M] }, $$
with the Hubbard-Stratonovich action
$$ \boxed{ S_\textrm{HS}[\mathcal M] = \frac{\lambda \mathcal M^2 }{2} - \ln \left( \int e^{-S_0[\varphi] + \lambda \mathcal M \mathcal O[\varphi]} \; \mathrm D[\varphi] \right) }. $$
3. Connection between mean field and Hubbard-Stratonovich. From the above, we can directly see that $S_\textrm{HS}[\mathcal M ] = - \ln Z_\textrm{mf}[\mathcal M] $. Equivalently, we can write
$$ \boxed{ Z = \int Z_\textrm{mf}[\mathcal M] \; \mathrm D[\mathcal M] } \; . $$
The saddle-point approximation is thus
$$ \boxed{ Z \approx e^{-S_\textrm{HS}[\mathcal M_0]} = Z_\textrm{mf}[\mathcal M_0] \qquad \textrm{with } S'_\textrm{HS}[\mathcal M_0] = 0 },$$
but note that the Hubbard-Stratonovich action being extremal is exactly equivalent to the extremality of the mean field free energy, i.e., $\mathcal M_0 = M$. QED.
4. Beyond mean field theory. What do we gain from this? We can directly incorporate the quadratic corrections around this saddle-point approximation. I.e., a better approximation is
$$ Z \approx \frac{Z_\textrm{mf}[M]}{\sqrt{S''_\textrm{HS}[M]}}. $$
In fact, this can be used to test how reliable/stable the mean field appoximation is.