In this paper we formulate and analyze a least squares boundary element method for the weakly singular boundary integral equation which is related to the solution of a Dirichlet boundary value problem for a second order partial differential equation, with the Laplacian as model problem. In particular we may assume less regular boundary data \(g \not \in H^{1/2}(\Gamma )\) but \(g \in L^2(\Gamma )\) . For this we consider the single layer boundary integral operator \(V: H^{-1}(\Gamma ) \rightarrow L^2(\Gamma )\) , i.e., we will solve the boundary integral equation \(Vw=f\) by minimizing \(\frac{1}{2} \, \Vert V w - f \Vert _{L^2(\Gamma )}^2\) . This results in a mixed variational formulation where we use piecewise constant approximations to discretize both the primal unknown \(w \in H^{-1}(\Gamma )\) and the adjoint \(p:=f-Vw \in L^2(\Gamma )\) . Using nested boundary element spaces \(S_H^0(\Gamma ) \subseteq S_h^0(\Gamma )\) we can prove stability and related error estimates for both the primal and adjoint approximations, \(w_H\) and \(p_h\) , respectively. When considering the approximate adjoint \(p_h\) on a finer mesh than the primal \(w_H\) , we can use \(\Vert p_h \Vert _{L^2(\Gamma )}\) as a posteriori error indicator for the error \(\Vert w-w_H \Vert _{H^{-1}(\Gamma )}\) to drive an adaptive mesh refinement. Note that this defines an adaptive boundary element method also for regular boundary data \(g \in H^{1/2}(\Gamma )\) . Numerical examples confirm the theoretical results.