Whereas there are many references on univariate boundary kernels, the construction of boundary kernels for multivariate density and curve estimation has not been investigated in detail. The use of multivariate boundary kernels ensures global consistency of multivariate kernel estimates as measured by the integrated mean-squared error or sup-norm deviation for functions with compact support. We develop a class of boundary kernels which work for any support, regardless of the complexity of its boundary. Our construction yields a boundary kernel for each point in the boundary region where the function is to be estimated. These boundary kernels provide a natural continuation of non-negative kernels used in the interior onto the boundary. They are obtained as solutions of the same kernel-generating variational problem which also produces the kernel function used in the interior as its solution. We discuss the numerical implementation of the proposed boundary kernels and their relationship to locally weighted least squares. Along the way we establish a continuous least squares principle and a continuous analogue of the Gauss–Markov theorem.