  Aims. To develop a fully Bayesian least squares deconvolution (LSD) that can
be applied to the reliable detection of magnetic signals in noise-limited
stellar spectropolarimetric observations using multiline techniques. Methods.
We consider LSD under the Bayesian framework and we introduce a flexible
Gaussian Process (GP) prior for the LSD profile. This prior allows the result
to automatically adapt to the presence of signal. We exploit several linear
algebra identities to accelerate the calculations. The final algorithm can deal
with thousands of spectral lines in a few seconds. Results. We demonstrate the
reliability of the method with synthetic experiments and we apply it to real
spectropolarimetric observations of magnetic stars. We are able to recover the
magnetic signals using a small number of spectral lines, together with the
uncertainty at each velocity bin. This allows the user to consider if the
detected signal is reliable. The code to compute the Bayesian LSD profile is
freely available.
