Here is a sketch of an approach that can can be generalized to much more extensive settings. A good reference is Michael E. Taylor's book Partial Differential Equations, Part 1, chapter 5, but it is quite standard and can be found probably in most books about linear elliptic PDEs.
Let $D:C^{\infty}(\Omega)\rightarrow\mathfrak{X}(\Omega)$ be a first order differential operator that assigns to each smooth function a vector field. Denote by $D^{*}:\mathfrak{X}(\Omega)\rightarrow C^{\infty}(\Omega)$ its $L^{2}$-dual, fully charactrized by operations on compactly supported functions. This gives rise to a boundary operator $\sigma:C^{\infty}(\Omega)\rightarrow \mathfrak{X}(\Omega)|_{\partial\Omega}$ such that,
$$
(Du,v)=(u,D^{*}v)+(\sigma u,v)_{\partial\Omega}
$$
e.g. when $D=\nabla$ then $D^{*}=\nabla^{*}$ is the divergence operator and $\sigma u=u|_{\partial\Omega} N$. The result is a second-order differential operator $D^{*}D:C^{\infty}(\Omega)\rightarrow C^{\infty}(\Omega)$. In the case where $D=\nabla$, we have Poincare's inequality which reads for $u$ with $\sigma u=0$:
$$
(D^{*}Du,u)=(Du,Du)\geq C\|u\|_{L^{2}}^{2}
$$
This is the sense in which $D^{*}D$ is a positive definite operator on $L^2(\Omega)$. The problem here is that $D^{*}D$ is not actually defined on the whole of $L^{2}(\Omega)$, only on the dense subset $H^{2}(\Omega)$. There are several ways to contain this, also taking into account the demand on the boundary $\sigma u=0$. The most general one, I think, is the method of unbounded operators, which you can find either Rudin's Functional Analysis chapter 13 or in the Appendix of Taylor's book. In any case, once you establish the manner in which $D^{*}D:L^{2}(\Omega)\rightarrow L^{2}(\Omega)$ is positive definite, which goes through the fact that it is "self adjoint", then you can produce a left inverse for it (up to finite dimensional effects, maybe), with nonnegative eigenvalues. This inverse, which is the so called Green Function, must be compact by Relich's theorem since it takes functions from $L^{2}(\Omega)$ to the domain of $D^{*}D$, which is $H^{2}(\Omega)$. This in turn implies that $I+\lambda G$ is always invertible.