  We apply Gaussian processes and Hubble function data in $f(T)$ cosmology, to
reconstruct for the first time the $f(T)$ form in a model-independent way. In
particular, using $H(z)$ datasets coming from cosmic chronometers as well as
from the radial BAO method, alongside the latest released local value $H_{0} =
73.52 \pm 1.62$ km/s/Mpc, we reconstruct $H(z)$ and its derivatives, resulting
eventually in a reconstructed region for $f(T)$, without any assumption.
Although the cosmological constant lies in the central part of the
reconstructed region, the obtained mean curve follows a quadratic function.
Inspired by this we propose a new $f(T)$ parametrization, i.e. $f(T) =
-2\Lambda +\xi T^2$, with $\xi$ the sole free parameter that quantifies the
deviation from $\Lambda$CDM cosmology. Additionally, we confront three viable
one-parameter $f(T)$ models of the literature, which respectively are the
power-law, the square-root exponential, and the exponential one, with the
reconstructed $f(T)$ region, and then we extract significantly improved
constraints for their model parameters, comparing to the constraints that arise
from usual observational analysis. Finally, we argue that since we are using
the direct Hubble measurements and the local value for $H_0$ in our analysis,
with the above reconstruction of $f(T)$, the $H_0$ tension can be efficiently
alleviated.
