We address the numerical approximation of Mean Field Games with local couplings. For powerlike Hamiltonians, we consider both the stationary system introduced in [51,53] and also a similar system involving density constraints in order to model hard congestion effects [65,57]. For finite difference discretizations of the Mean Field Game system as in [3], we follow a variational approach. We prove that the aforementioned schemes can be obtained as the optimality system of suitably defined optimization problems. In order to prove the existence of solutions of the scheme with a variational argument, the monotonicity of the coupling term is not used, which allow us to recover general existence results proved in [3]. Next, assuming next that the coupling term is monotone, the variational problem is cast as a convex optimization problem for which we study and compare several proximal type methods. These algorithms have several interesting features, such as global convergence and stability with respect to the viscosity parameter, which can eventually be zero. We assess the performance of the methods via numerical experiments.