Modeling and identification of high dimensional systems, involving signals with many components, poses severe challenges to off-the-shelf techniques for system identification. This is particularly so when relatively small data sets, as compared to the number signal components, have to be used. It is often the case that each component of the measured signal can be described in terms of a few other measured variables and these dependences can be encoded in a graphical way via so called "Dynamic Bayesian Networks". The problem of finding the interconnection structure as well as estimating the dynamic models can be posed as a system identification problem which involves variable selection. While this variable selection could be performed via standard selection techniques, computational complexity may however be a critical issue, being combinatorial in the number of inputs and outputs. In this paper we introduce two new nonparametric techniques which borrow ideas from a recently introduced kernel estimator called "stable-spline" as well as from sparsity inducing priors which use l(1)-type penalties. Numerical experiments regarding estimation of large scale sparse (ARMAX) models show that this technique provides a definite advantage over a group LAR algorithm and state-of-the-art parametric identification techniques based on prediction error minimization. (C) 2012 Elsevier Ltd. All rights reserved.