We define and analyse a least-squares finite element method for a first-order reformulation of the obstacle problem. Moreover, we derive variational inequalities that are based on similar but non-symmetric bilinear forms. A priori error estimates including the case of non-conforming convex sets are given and optimal convergence rates are shown for the lowest-order case. We provide also a posteriori bounds that can be be used as error indicators in an adaptive algorithm. Numerical studies are presented.Date: October 14, 2018. 2010 Mathematics Subject Classification. 65N30, 65N12, 49J40. Key words and phrases. First-order system, least-squares method, variational inequality, obstacle problem, a priori analysis, a posteriori analysis.Acknowledgment. This work was supported by CONICYT through FONDECYT project "Least-squares methods for obstacle problems" under grant 11170050. 1 arXiv:1801.09622v1 [math.NA] 29 Jan 2018on. For simplicity assume a zero obstacle (the remainder of the paper deals with general non-zero obstacles). Then, the problem reads −∆u ≥ f, u ≥ 0, (−∆u − f )u = 0 in some domain Ω and u| ∂Ω = 0. Introducing the Lagrange multiplier (or reaction force) λ = −∆u − f and σ = ∇u, we rewrite the problem as a first-order system, see also [2,3,9,16],Note that f ∈ L 2 (Ω) does not imply more regularity for u so that λ ∈ H −1 (Ω) is only in the dual space in general. However, observe that div σ + λ = −f ∈ L 2 (Ω) and therefore the functionalwhere · , · denotes a duality pairing, is well-defined for div σ + λ ∈ L 2 (Ω). We will show that minimizing J over a convex set with the additional linear constraints u ≥ 0, λ ≥ 0 is equivalent to solving the obstacle problem. We will consider the variational inequality associated to this problem with corresponding bilinear form a(·, ·). An issue that arises is that a(·, ·) is not necessarily coercive. However, as it turns out, a simple scaling of the first term in the functional ensures coercivity on the whole space. In view of the aforementioned properties, this means that our method is unconstrained stable. The recent work [16] based on a Lagrange formulation (without reformulation to a first-order system) considers augmenting the trial spaces with bubble functions (mixed method) resp. adding residual terms (stabilized method) to obtain stability.Furthermore, we will see that the functional J evaluated at some discrete approximation (u h , σ h , λ h ) with u h , λ h ≥ 0 is an upper bound for the error. Note that for λ h ∈ L 2 (Ω) the duality λ h , u h reduces to the L 2 inner product. Thus, all the terms in the functional can be localized and used as error indicators.Additionally, we will derive and analyse other variational inequalities that are also based on the first-order reformulation. The resulting methods are quite similar to the least-squares scheme since they share the same residual terms. The only difference is that the compatibility condition λu = 0 is incorporated in a different, non-symmetric, way. We will present a uniform analysis that covers the least-squa...