We present a multi-level scheme to solve the linear system that arises when one attempts to find the steady state solution of the kinetic model equations using implicit time stepping in the discrete velocity method (DVM). At each physical node, we observe that the relaxation state is parameterized by macroscopic quantities derived from the local velocity distribution function (VDF). Since the number of macroscopic quantities is usually much smaller compared to the number of velocity nodes, the gradient associated with the relaxed VDF can be written as the product of low-rank matrices. Using the Woodbury matrix identity, all the velocity components of the local VDF can be updated simultaneously. This means at a local scale, the evolution of VDF can be fully implicit, and free from the restrictions in time step size even if the Knudsen number is relatively small. From a global perspective, we notice that the linear system is sparse in block sense. The diagonal blocks are dense due to relaxation, and the off-diagonal blocks are diagonal, if not zero, due to convection. Since the diagonal blocks can be inverted using the Woodbury identity, the global implicit system can be solved using a block-relaxation method (such as block Gauss-Seidel). We upgrade this block-relaxation method to a multi-level one, reducing the number of blocks by half for each inner iteration. As a result, both high and low frequency error modes can be effectively reduced. We observe that using singular value decomposition, the low-rank property in the diagonal blocks can be preserved. Therefore, the numerical complexity of one multi-level iteration is of the same order compared to one semi-implicit sweep in conventional DVM. Numerical examples show that the rate of convergence is much higher and almost independent of the Knudsen number.