  The increase of discrepancy in the standard procedure to choose the arbitrary
functional form of the Lagrangian $f(Q)$ motivates us to solve this issue in
modified theories of gravity. In this regard, we investigate the Gaussian
process (GP), which allows us to eliminate this issue in a $f(Q)$
model-independent way. In particular, we use the 57 Hubble measurements coming
from cosmic chronometers and the radial Baryon acoustic oscillations (BAO) to
reconstruct $H(z)$ and its derivatives $H'(z)$, $H''(z)$, which resulting lead
us to reconstruct region of $f(Q)$, without any assumptions. The obtained mean
curve along $\Lambda$CDM constant in the reconstructed region follows a
quadratic behavior. This motivates us to propose a new $f(Q)$ parametrization,
i.e., $f(Q)= -2\Lambda+ \epsilon Q^2$, with the single parameter $\epsilon$,
which signifies the deviations from $\Lambda$CDM cosmology. Further, we probe
the widely studied power-law and exponential $f(Q)$ models against the
reconstructed region and can improve the parameter spaces significantly
compared with observational analysis. In addition, the direct Hubble
measurements, along with the reconstructed $f(Q)$ function, allow the $H_0$
tension to be alleviated.
