Summary: The Hybrid Spline method \((H\)-spline) is a method of density estimation which involves regression splines and smoothing splines methods. Using basis functions \((B\)-splines), this method is much faster than Smoothing Spline Density Estimation approach. Simulations suggest that with more structured data (e.g., several modes) \(H\)-spline method estimates the modes as well as Logspline. The \(H\)-spline algorithm is designed to compute a solution to the penalized likelihood problem. The smoothing parameter is updated jointly with the estimate via a cross-validation performance estimate, where the performance is measured by a proxy of the symmetrized Kullback-Leibler. The initial number of knots is determined automatically based on an estimate of the number of modes and the symmetry of the underlying density. The algorithm increases the number of knots by 1 until the symmetrized Kullback-Leibler distance, based on two consecutives estimates, satisfies a condition which was determined empirically.