  A finite element methodology for large classes of variational boundary value
problems is defined which involves discretizing two linear operators: (1) the
differential operator defining the spatial boundary value problem; and (2) a
Riesz map on the test space. The resulting linear system is overdetermined. Two
different approaches for solving the system are suggested (although others are
discussed): (1) solving the associated normal equation with linear solvers for
symmetric positive-definite systems (e.g. Cholesky factorization); and (2)
solving the overdetermined system with orthogonalization algorithms (e.g. QR
factorization). The finite element assembly algorithm for each of these
approaches is described in detail. The normal equation approach is usually
faster for direct solvers and requires less storage. The second approach
reduces the condition number of the system by a power of two and is less
sensitive to round-off error. The rectangular stiffness matrix of second
approach is demonstrated to have condition number $\mathcal{O}(h^{-1})$ for a
variety of formulations of Poisson's equation. The stiffness matrix from the
normal equation approach is demonstrated to be related to the monolithic
stiffness matrices of least-squares finite element methods and it is proved
that the two are identical in some cases. An example with Poisson's equation
indicates that the solutions of these two different linear systems can be
nearly indistinguishable (if round-off error is not an issue) and rapidly
converge to each other. The orthogonalization approach is suggested to be
beneficial for problems which induce poorly conditioned linear systems.
Experiments with Poisson's equation in single-precision arithmetic as well as
the linear acoustics problem near resonance in double-precision arithmetic
verify this conclusion.
