  We consider the numerical solution of an abstract operator equation $Bu=f$ by
using a least-squares approach. We assume that $B: X \to Y^*$ is an
isomorphism, and that $A : Y \to Y^*$ implies a norm in $Y$, where $X$ and $Y$
are Hilbert spaces. The minimizer of the least-squares functional $\frac{1}{2}
\, \| Bu-f \|_{A^{-1}}^2$, i.e., the solution of the operator equation, is then
characterized by the gradient equation $Su=B^* A^{-1}f$ with an elliptic and
self-adjoint operator $S:=B^* A^{-1} B : X \to X^*$. When introducing the
adjoint $p = A^{-1}(f-Bu)$ we end up with a saddle point formulation to be
solved numerically by using a mixed finite element method. Based on a discrete
inf-sup stability condition we derive related a priori error estimates. While
the adjoint $p$ is zero by construction, its approximation $p_h$ serves as a
posteriori error indicator to drive an adaptive scheme when discretized
appropriately. While this approach can be applied to rather general equations,
here we consider second order linear partial differential equations, including
the Poisson equation, the heat equation, and the wave equation, in order to
demonstrate its potential, which allows to use almost arbitrary space-time
finite element methods for the adaptive solution of time-dependent partial
differential equations.
