Why is it so coincident that Palatini variation of Einstein-Hilbert action will obtain an equation that connection is Levi-Civita connection? There are two ways to do the variation of Einstein-Hilbert action.
First one is Einstein formalism which takes only metric independent. After variation of action, we get the Einstein field equation.
Second one is Palatini formalism which takes metric and connection are independent. After variation, we get two equations, first is field equation and second is that connection is Levi-Civita connection.
So my question is why it is so coincident that Palatini variation of Einstein-Hilbert action will obtain an equation that connection is Levi-Civita connection and Palatini formalism coincides with Einstein formalism? While for $f(R)$ action they are generally different. Are there some deeper mathematical or physical structures of Einstein-Hilbert action which can account for it.
 A: I) In Palatini $f(R)$ gravity, the Lagrangian density is$^1$ 
$$ {\cal L}(g,\Gamma)~=~ \frac{1}{2\kappa}\sqrt{-g} f(R) + {\cal L}_{\rm m}; \tag{1}$$
with matter Lagrangian density ${\cal L}_{\rm m}$;
with scalar curvature
$$R~:=~ g^{\mu\nu} R_{\mu\nu}(\Gamma);\tag{2}$$ 
with Ricci curvature $R_{\mu\nu}(\Gamma)$; and where 
$$\Gamma^{\lambda}_{\mu\nu}~=~\Gamma^{\lambda}_{\nu\mu}\tag{3}$$ 
is an arbitrary torsionfree$^2$ connection. 
II) As OP mentions, the word Palatini refers to that the metric $g_{\mu\nu}$ and the connection $\Gamma^{\lambda}_{\mu\nu}$ are independent variables$^3$. We therefore get two types of EL equations:


*

*The EL equations 
$$ f^{\prime}(R)R_{\mu\nu} -\frac{1}{2}f(R)g_{\mu\nu}~\stackrel{(1)+(5)}{\approx}~\kappa T_{\mu\nu} \tag{4}$$
for the metric $g_{\mu\nu}$ are the generalization of EFE, where 
$$T^{\mu\nu}~:=~\frac{2}{\sqrt{-g}}\frac{\delta S_{\rm m}}{\delta g_{\mu\nu}} \tag{5} $$
is the matter Hilbert stress-energy-momentum (SEM) tensor. 
[In eq. (4) the $\approx$ symbol means equality modulo equations of motion. In this answer, we use $(-,+,\ldots,+)$ Minkowski sign convention in $d$ spacetime dimensions.]

*If the matter action $S_{\rm m}$ doesn't depend on the connection $\Gamma^{\lambda}_{\mu\nu}$, then the EL equations 
$$  \nabla_{\lambda}\hat{\mathfrak{g}}^{\mu\nu}
~\stackrel{(1)}{\approx}~0,\qquad 
\hat{\mathfrak{g}}^{\mu\nu}
~:=~\sqrt{-g} f^{\prime}(R) g^{\mu\nu}
~\stackrel{(8)}{=}~\sqrt{-\hat{g}}  \hat{g}^{\mu\nu}, \tag{6}$$
for the connection $\Gamma^{\lambda}_{\mu\nu}$ turn out to be the metric compatibility condition 
$$\nabla_{\lambda}\hat{g}_{\mu\nu} ~\stackrel{(6)+(8)}{\approx}~0\tag{7}$$
for a conformally equivalent metric
$$ \hat{g}_{\mu\nu}~:=~f^{\prime}(R)^{\frac{2}{d-2}} g_{\mu\nu}, \tag{8}$$
known as the Einstein frame metric.
In other words, the classical solution for $\Gamma^{\lambda}_{\mu\nu}$ is the Levi-Civita connection for the Einstein frame metric $\hat{g}_{\mu\nu}$.
III) So Einstein gravity (GR) with a possible cosmological constant 
$$ f(R)~=~R-2\Lambda,\tag{9} $$
or equivalently
$$ f^{\prime}(R)~=~1,\tag{10}$$ 
corresponds to the special case where the two metrics $g_{\mu\nu}$ and $\hat{g}_{\mu\nu}$ coincide, and hence $\Gamma^{\lambda}_{\mu\nu}$ becomes the Levi-Civita connection for $g_{\mu\nu}$.
--
$^1$ It is natural to replace the Lagrangian density (1) with the extended Lagrangian density
$$ \tilde{\cal L}(g,\Gamma,\Phi)~=~ \frac{1}{2\kappa}\sqrt{-g}\{\Phi R-V(\Phi)\} + {\cal L}_{\rm m}; \tag{11}$$
with auxiliary scalar dilaton field $\Phi$; and where the potential
$$V(\Phi)~:=~\sup_r(\Phi r -f(r))\tag{12}$$
is the Legendre transform of the function $f$. If we integrate out the auxiliary scalar field $\Phi$, we then return to the $f(R)$ Lagrangian density (1) which we started from! The EL equations 
$$  \nabla_{\lambda}\hat{\mathfrak{g}}^{\mu\nu}
~\stackrel{(1)}{\approx}~0,\qquad 
\hat{\mathfrak{g}}^{\mu\nu}
~:=~\sqrt{-g} \Phi g^{\mu\nu}
~\stackrel{(14)}{=}~\sqrt{-\hat{g}}  \hat{g}^{\mu\nu}, \tag{13}$$
for the connection $\Gamma^{\lambda}_{\mu\nu}$ become the metric compatibility condition (7) for the Einstein frame metric
$$ \hat{g}_{\mu\nu}~:=~\Phi^{\frac{2}{d-2}} g_{\mu\nu}. \tag{14}$$
After the connection $\Gamma^{\lambda}_{\mu\nu}$ has been integrated out, the Lagrangian density (11) becomes
$$ \tilde{\cal L}(g,\Phi)~=~ \frac{1}{2\kappa}\sqrt{-\hat{g}}\{R(\hat{g})-\Phi^{\frac{d}{2-d}}V(\Phi)\} + {\cal L}_{\rm m} , \tag{15}$$
where eq. (14) is implicitly assumed. 
However, unfortunately the Legendre transform $V$ does not exist for Einstein gravity (9), so we will not consider the extended Lagrangian density (11) further in this answer.
$^2$ One could allow a non-dynamical torsion piece, but we will not pursue this here for simplicity. For more on torsion, see e.g. also this Phys.SE post.
$^3$ Normally in non-Palatini formulations, we integrate out the connection $\Gamma^{\lambda}_{\mu\nu}$ and keep the metric $g_{\mu\nu}$. Eddington & Schrödinger proposed the opposite! Let us analyze this possibility here. Define for later convenience a double-index notation $M=\mu\mu^{\prime}$ and the following shorthand notation
$$\frac{f(R)}{2f^{\prime}(R)}~=:~\hat{f}(R) ~\equiv~ \hat{f}_0+\hat{f}_1R +\hat{f}_2(R).\tag{16}$$ 
Let us only consider the vacuum 
$$T_{\mu\nu}~=~0.\tag{17}$$ 
from now on. We then have 
$$g^M~\stackrel{(4)+(16)+(17)}{\approx}~\hat{f}(R)R^M,\tag{18}$$ 
where $R^M$ is the inverse Ricci tensor. Equivalently, we have
$$\left(\delta^M_N - \hat{f}_1 R^M R_N\right) g^N~\stackrel{(16)+(18)}{\approx}~\left(\hat{f}_0+f_2(R) \right)R^M.\tag{19}$$
So we get a fixed-point equation for the inverse metric
$$g^N~\stackrel{(19)}{\approx}~\left(\delta^N_M +\frac{\hat{f}_1}{1-d\hat{f}_1} R^N R_M\right)\left(\hat{f}_0+\hat{f}_2(R) \right)R^M$$
$$~=~ \frac{1}{1-d\hat{f}_1}\left(\hat{f}_0+\hat{f}_2\left(g^MR_M\right) \right)R^N.\tag{20}$$
Let us specialize to Einstein gravity (9). Then  
$$\hat{f}_0~=~-\Lambda;\qquad\hat{f}_1~=~\frac{1}{2};\qquad\hat{f}_2(R)~=~0.\tag{21}$$ 
The inverse metric becomes 
$$g^N~\stackrel{(20)+(21)}{\approx}~\frac{2\Lambda}{d-2}R^N.\tag{22}$$ 
And hence 
$$R~\stackrel{(22)}{\approx}~\frac{2d}{d-2}\Lambda,\tag{23}$$ 
and 
$$g_{\mu\nu}~\stackrel{(2)+(22)}{\approx}~\frac{d-2}{2\Lambda}R_{\mu\nu}(\Gamma).\tag{24}$$ 
So the EH Lagrangian density becomes Born-Infeld-like: 
$${\cal L}(\Gamma)~\stackrel{(1)+(17)+(23)+(24)}{\approx}~\frac{1}{\kappa}\left(\frac{d-2}{2\Lambda}\right)^{\frac{d}{2}-1}\sqrt{-\det(R_{\mu\nu}(\Gamma))}.\tag{25}$$ 
Note that the Eddington-Schrödinger action (25) only works for a non-zero cosmological constant $\Lambda\neq 0$.
