Given an undirected graph and a size parameter $k$, the Densest $k$-Subgraph (D$k$S) problem extracts the subgraph on $k$ vertices with the largest number of induced edges. While D$k$S is NP--hard and difficult to approximate, penalty-based continuous relaxations of the problem have recently enjoyed practical success for real-world instances of D$k$S. In this work, we propose a scalable and exact continuous penalization approach for D$k$S using the error bound principle, which enables the design of suitable penalty functions. Notably, we develop new theoretical guarantees ensuring that both the global and local optima of the penalized problem match those of the original problem. The proposed penalized reformulation enables the use of first-order continuous optimization methods. In particular, we develop a non-convex proximal gradient algorithm, where the non-convex proximal operator can be computed in closed form, resulting in low per-iteration complexity. We also provide convergence analysis of the algorithm. Experiments on large-scale instances of the D$k$S problem and one of its variants, the Densest ($k_1, k_2$) Bipartite Subgraph (D$k_1k_2$BS) problem, demonstrate that our method achieves a favorable balance between computation cost and solution quality.