The fundamental nature of flow instability in wall bounded flows changes with Mach number. The objectives of this study are two-fold, (i) compute the instability modes in high Mach number Poiseuille flows using linear stability analysis (LSA) and, (ii) perform direct numerical simulations (DNS) of the instability development using a solver based on gas kinetic method (GKM) for the purpose of code validation by comparison against LSA results. The LSA and DNS are performed for the case of Poiseuille flow over a range of Mach numbers – from moderately supersonic to hypersonic speeds. First, LSA is employed to identify the most unstable mode over the range of Mach numbers. We then perform two sets of GKM-DNS to corroborate the LSA results over the Mach number range. In the first set of simulations, the background field is initially perturbed with the most unstable mode identified by LSA and the evolution is monitored. It is shown that GKM-DNS accurately captures the exponential growth in kinetic energy for all Mach numbers. The second set of GKM-DNS simulations is performed by superposing the background pressure field with random initial perturbations. After an initial transient period, the modes predicted by LSA dominate the DNS flow field evolution. The wave-vector and mode shapes of the dominant instability are well replicated by GKM-DNS at each Mach number. These insights in the linear regime of high speed Poiseuille flow and validation of GKM are important for understanding and simulating wall bounded flows.