This article proposes a scalable, matrix-free approach to kernel-based regularization for finite impulse response estimation. Our methodology is based on Bayesian optimization, a gradient-free global optimization strategy that can handle noisy objective function evaluations and is suitable for indirect or stochastic hyperparameter estimation. A key assumption is that the kernel function only has a small number of parameters and that matrix-vector products with the kernel matrix are inexpensive. Popular kernels, such as the tuned-correlated (TC), diagonal-correlated (DC), and stable spline (SS) kernels, meet these criteria. To achieve scalability, we simultaneously exploit the structure of both the kernel matrix and the regressor matrix, which reduces memory requirements and computational cost. Our empirical results, which are based on randomly generated systems, demonstrate that the computational cost scales approximately linearly with the model order for the TC, DC, and SS kernels.