Bayesian linear models are widely used as efficient approaches for nonparametric function estimation. In this paper, we present a Bayesian method for generating finite-dimensional linear models that can handle large datasets. This method is based on an efficient Markov chain Monte Carlo algorithm. The advantage of this approach is that sampling is performed before conditioning, rather than after. This enables the use of efficient samplers when the prior covariance matrix exhibits special properties, such as being Toeplitz, block-Toeplitz, or sparse. The performance of the proposed approach is evaluated in the context of nonparametric function estimation with large datasets. A numerical comparison with direct approaches based on Cholesky factorization is provided to illustrate its efficiency in terms of computational runtime.