Abstract: We propose a computationally efficient method for data-driven adaptive spline smoothing. Finite differences are adopted to characterize the roughness pattern of the underlying function. These differences are then clustered into several groups, each of which is modeled to have a separate penalty parameter. State space representation is developed for implementation. Simulation shows that the proposed method is fast in computation, with median computational time as 5%~30% of existing methods. It also shows that the proposed method works well in several typical functions for adaptive smoothing with respect to mean square errors. Application to a shock wave lithotripsy data example shows that the proposed method generates function estimates that agree with the corresponding physical properties.