In this paper we present the numerical analysis of a three-field stabilized finite element formulation recently proposed to approximate viscoelastic flows. The three-field viscoelastic fluid flow problem may suffer from two types of numerical instabilities: on the one hand we have the two inf-sup conditions related to the mixed nature problem and, on the other, the convective nature of the momentum and constitutive equations may produce global and local oscillations in the numerical approximation. Both can be overcome by resorting from the standard Galerkin method to a stabilized formulation. The one presented here is based on the subgrid scale concept, in which unresolvable scales of the continuous solution are approximately accounted for. In particular, the approach developed herein is based on the decomposition into their finite element component and a subscale, which is approximated properly to yield a stable formulation. The analyzed problem corresponds to a linearized version of the Navier–Stokes/Oldroyd-B case where the advection velocity of the momentum equation and the non-linear terms in the constitutive equation are treated using a fixed point strategy for the velocity and the velocity gradient. The proposed method permits the resolution of the problem using arbitrary interpolations for all the unknowns. We describe some important ingredients related to the design of the formulation and present the results of its numerical analysis. It is shown that the formulation is stable and optimally convergent for small Weissenberg numbers, independently of the interpolation used.