We propose a new mathematical model for $n-k$-dimensional non-linear correlations with intrinsic scatter in $n$-dimensional data. The model is based on Riemannian geometry, and is naturally symmetric with respect to the measured variables and invariant under coordinate transformations. We combine the model with a Bayesian approach for estimating the parameters of the correlation relation and the intrinsic scatter. A side benefit of the approach is that censored and truncated datasets and independent, arbitrary measurement errors can be incorporated. We also derive analytic likelihoods for the typical astrophysical use case of linear relations in $n$-dimensional Euclidean space. We pay particular attention to the case of linear regression in two dimensions, and compare our results to existing methods. Finally, we apply our methodology to the well-known $M_\text{BH}$-$σ$ correlation between the mass of a supermassive black hole in the centre of a galactic bulge and the corresponding bulge velocity dispersion. The main result of our analysis is that the most likely slope of this correlation is $\sim 6$ for the datasets used, rather than the values in the range $\sim 4$-$5$ typically quoted in the literature for these data.