This paper presents a computational method, with several variants, for fitting bicubic splines by least squares to data given at arbitrary points. Products of B-splines are used in the representation of the bicubic splines. The resulting observation equations are solved by means of Householder transformations. A stable method for imposing linear equality constraints is also described. The methods take account of rank-deficiency and are readily extended to more dimensions.