A new method of multivariate smooth fitting of scattered, noisy data using cubic B-splines was developed. An optimum smoothing function was defined to minimize the l 2 norm composed of the data residuals and the first and the second derivatives, which represent the total misfit, fluctuation, and roughness of the function, respectively. The function is approximated by a cubic B-spline expansion with equispaced knots. The solution can be interpreted in three ways. From the stochastic viewpoint, it is the maximum-likelihood estimate among the admissible functions under the a priori information that the first and second derivatives are zero everywhere due to random errors, i.e., white noise. From the physical viewpoint, it is the finite-element approximation for a lateral displacement of a bar or a plate under tension which is pulled to the data points by springs. From a technical viewpoint, it is an improved spline-fitting algorithm. The additional condition of minimizing the derivative norms stabilizes the linear equation system for the expansion coefficients.