This a very interesting problem, and I also wondered regarding these discrete conditions, a while ago. This is more of a comment/suggestion,
Firstly, keeping in mind the well-posedness of Elliptic PDE, I would start working with the case of $\lambda \geq 0$. And to make things easier have a Dirichlet boundary, Periodic boundary wouldn't be very different, we will the have to consider the point on the boundary of the given space to also be active nodes, but there is a subtle problem here, how do we define connectivities of these nodes with the other boundary nodes? By the taking the nodes from the internal of the other boundary being place at the same spacing on the empty side of the other periodic boundary, like ghost points and then making the connectivities.
Second, taking in the conditions to be satisfied internally also as boundary conditions, we have a domain with "holes" (even then it is a connected domain) and it is a Lipschitz domain, so Lax-Milgram theorem can be applied, resulting in a well-posed problem. I am pretty sure the solution in our case would be continuous and has weak derivatives from the well-posedness of the variational form, I am not sure of elliptic regularity possibly because of the "holes" but I feel it should be there as well!?
Third, for the numerical solution, I would start working with conformal FEM or the Ritz-Galerkin method, but making sure our discretization or Triangulation of the space is such that, these given point are nodes and not internal to triangles (lets call them special nodes), so that it is conformal and the related conditions can be easily applied. Now for simplicity, if we take up only linear basis, we write the numerical solution in terms of the basis of hat functions, see that the special nodes also have hat functions associated whose coefficients are known from the given conditions and be put on the RHS. The number of unknowns would be the rest of the nodes, and a system can be formed of that required size by testing with the hat functions (the remaining ones). And, if the triangulation is conformally (i.e. special nodes are nodes and not internal to triangles) refined, I would expect convergence similar to the general advection-diffusion problem.
Although the abstraction is the same, it would get complicated to implement for higher dimensional case, like it is the case usually.
Regarding references: I searched for these "hole" kind of dirichlet boundaries but didn't find any!