Optimizing $\epsilon$ in $\epsilon$-kernel
Jeff Philips has a related result showing the stability of kernels: http://www.cs.utah.edu/~jeffp/papers/stable-kernel.pdf.

As for the problem itself, it can be formulated as a hitting set problem. Indeed, you need to pick a point into the kernel for each "complement of a slab" that contains $\geq \epsilon_k$ fraction of the input points. Here $\epsilon_k$ is the minimum $\epsilon$ that has a kernel of size $k$. Since this is just a hitting set problem, it can be approximated up to $O( \log k)$ factor in geometric settings. You need to guess the value of $\epsilon_k$, but naively just do a binary search on epsilon. As such, one can get two results:


A set of size $O( k \log k)$ that is an $\epsilon_k$ kernel.
A set of size k that is  an $\epsilon_{O(k/\log k)}$ kernel.


The algorithm for computing the set cover is the reweighting algorithm. In fact, the original paper by Clarkson solves a related problem about polytope approximation:

@inproceedings{c-apca-93,
  author       = "Kenneth L. Clarkson",
  title        = "Algorithms for Polytope Covering and Approximation",
  booktitle    = "Proc. 3rd Workshop Algorithms Data Struct.",
  series       = LNCS,
  volume       = 709,
  publisher    = "Springer-Verlag",
  year         = 1993,
  pages        = "246--252"
}

As for what to implement, i would use an incremental algorithm that performs quite well in practice, and add points an incremental fashion. See for example, the paper "Robust Shape Fitting via Peeling and Grating Coresets" and references therein.

it seems believable BTW that one should be able to do constant factor approximation instead of logarithmic.

I dont know of a formal proof that this problem is NP-hard, but it would be very surprising if it is not, IMHO.

The problem is probably solvable exactly in the plane using dynamic programming, I think....