Not so fast, there! You're jumping to conclusions. Who said you couldn't write the action as quadratic in a gauge field strength?
More precisely, though you may not necessarily be able to do that for the Einstein-Hilbert action without a cosmological coefficient, that's not believed to be physically relevant anymore: there's a non-zero cosmological coefficient present. You can do it, if the cosmological coefficient is not zero, which is the physically relevant case.
In that case, the requirement that the action be quadratic in the field strengths carries, with it, the prediction (after the fact, so "retrodiction") that the cosmological coefficient be non-zero.
An Einstein-Hilbert action with a cosmological coefficient has the form:
$$S = \int ε_{abcd} \left(K Ω^{ab} + L θ^aθ^b\right)θ^cθ^d,$$
for some coefficients $K$ and $L$, where $θ^a = h^a_μ dx^μ$ are the frame one-forms, $Ω^{ab} = ½ {R^{ab}}_{μν} dx^μ dx^ν$ the curvature two-forms, indexes are raised and lowered with the frame metric $η_{ab}$ and $ε_{abcd}$ is completely anti-symmetric with $ε_{0123} = 1$.
For convenience I'm writing Grassmann products as juxtaposition without the wedge here and below; e.g. $θ^aθ^b$, instead of $θ^a ∧ θ^b$; i.e. as anti-symmetric products $θ^aθ^b = -θ^bθ^a$, in the case of odd-degree forms.
The Cartan structure equations, with upstairs indexes, are:
$$dθ^a + ω^{ab}η_{bc}θ^c = Θ^a, \hspace 1em dω^{ab} + ω^{ac}η_{cd}ω^{db} = Ω^{ab},$$
where
$$ω^{ab} = Γ^α_{μc}η^{cb}dx^μ, \hspace 1em Θ^a = ½ T^a_{μν} dx^μdx^ν$$
are, respectively, the connection one-forms and and torsion two-forms.
In Riemann-Cartan geometries, the frame metric $η_{ab}$(which is constant) is non-degenerate, with inverse $η^{ab}$ and is constrained to have zero covariant derivatives, which entails the following equations:
$$Η^{ab} ≡ dη^{ab} + ω^{ab} + ω^{ba} = 0 \hspace 1em⇒\hspace 1em 0 = dH^{ab} + ω^{ac}η_{cd}H^{db} + ω^{bc}η_{cd}H^{ad} = Ω^{ab} + Ω^{ba},$$
and makes the connection one-form and curvature two-form anti-symmetric in their frame indices: $ω^{ab} = -ω^{ba}$ and $Ω^{ab} = -Ω^{ba}$.
This also entails zero covariant derivative on the space-time metric $g_{μν} = η_{ab}h^a_μh^b_ν$.
You can either constrain the torsion to be zero, in which case you have the Palatini action with a cosmological coefficient - which is just Einstein-Hilbert, itself, (with cosmological coefficient) transcribed from Riemannian geometry to Riemann-Cartan geometry, or else you can leave it unconstrained, in which case you have the Einstein-Cartan action with cosmological coefficient ... where zero torsion will follow dynamically, at least outside of matter sources.
The gauge potential, here, is
$$A = θ^aP_a + ½ω^{ab}S_{ab},$$
where we'll adopt the convention of treating products of the underlying Lie basis $P_a$ and $S_{ab}$ as freely interspersed with $dx^μ$'s, for convenience, and will work within an algebra where the Lie bracket can be expressed as $[u,v] = uv - vu$. Then, we may write
$$A^2 = ½θ^aθ^c\left[P_a,P_c\right] + ⅓θ^aω^{cd}\left[P_a,S_{cd}\right] + ¼ω^{ab}θ^c\left[S_{ab},P_c\right] + ⅛ω^{ab}ω^{cd}\left[S_{ab},S_{cd}\right],$$
with the anti-symmetry of the form products, e.g.
$$θ^aP_aθ^cP_c = θ^aθ^cP_aP_c = ½\left(θ^aθ^c - θ^cθ^a\right)P_aP_c = ½θ^aθ^c\left(P_aP_c - P_cP_a\right) = ½θ^aθ^c\left[P_a,P_c\right].$$
Adopting the following Lie brackets
$$\begin{align}
\left[P_a,P_c\right] &= λS_{ac}, & \left[P_a,S_{cd}\right] &= η_{ac}P_d - η_{ad}P_c, \\
\left[S_{ab},P_c\right] &= η_{bc}P_a - η_{ac}P_b, & \left[S_{ab},S_{cd}\right] &= η_{bc}S_{ad} - η_{bd}S_{ac} - η_{ac}S_{bd} + η_{ad}S_{bc},
\end{align}$$
for some constant $λ$, then for the field strength, we have:
$$F = dA + A^2 = Θ^aP_a + ½\left(Ω^{ab} + λθ^aθ^b\right)S_{ab}.$$
Therefore, an action of a form quadratic in the field strengths
$$\begin{align}
S &= \int M ε_{abcd}\left(Ω^{ab} + λθ^aθ^b\right)\left(Ω^{cd} + λθ^cθ^d\right) \\
&= \int ε_{abcd}\left(M Ω^{ab}Ω^{cd} + 2Mλ Ω^{ab}θ^cθ^d + Mλ^2 θ^aθ^bθ^cθ^d\right)
\end{align}$$
reduces equivalently to
$$
S = \int ε_{abcd}\left(2Mλ Ω^{ab}θ^cθ^d + Mλ^2 θ^aθ^bθ^cθ^d\right)
= \int ε_{abcd}\left(K Ω^{ab}θ^cθ^d + L θ^aθ^bθ^cθ^d\right)
$$
provided that
$$M = \frac{K^2}{4L}, \hspace 1em λ = \frac{2L}{K}.$$
because $ε_{abcd}Ω^{ab}Ω^{cd}$ is a boundary term and drops out from the analysis.
I'm not the only one to notice this. Just name the Lie algebra, A, whose Lie brackets I just listed above and do a web search on "A-gravity", whatever "A" might be. There's going to be something out there on it.
This also provides room for further expansion. You can just as well consider more fully quadratic actions of the form:
$$S = \int \left(M ε_{abcd} + N η_{ac}η_{bd}\right)\left(Ω^{ab} + λθ^aθ^b\right)\left(Ω^{cd} + λθ^cθ^d\right) + O η_{ab}Θ^aΘ^b$$
and see what results from that. Now you have an "Immirzi" coefficient somewhere in there, too.