This paper presents a stabilized finite element method (FEM) and a spherical harmonics method to discretize the space and angle of the Boltzmann transport equation. The FEM is based on the subgrid-scale (SGS) model, which decomposes the unknowns into resolvable scale and SGS with an approximation for the SGS and then embeds it into a resolvable scale formulation, which yields a stabilized variational formula with only a resolvable scale. In this method, the SGS is identified as the residual of the flux, which represents the indistinguishable high-frequency component. This method is characterized by a residual equation proposed on the subgrid, thus reflecting the relationship between the residual of the flux and the residual of the source. A simple assumption is proposed that the residual of the flux is the scaling of the residual of the source. The scaling parameter is identified as a stabilization parameter, and it takes the inverse of the norm of the transport operator. This method has been verified by various benchmark problems, and the numerical results show that it has high accuracy, stability, and void applicability.