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 H0 = 73.52 ± 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Λ + ξT 2 , with ξ the sole free parameter that quantifies the deviation from Λ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 H0 in our analysis, with the above reconstruction of f (T ), the H0 tension can be efficiently alleviated.