We study variational methods for rigorous approximation of invariant densities for a nonsingular map $T$ on a Borel measure space. The general method takes the form of a convergent sequence of optimization problems on $L^p$, $1 \leq p < \infty$ with a convex objective and finite moment constraints. Provided $T$ admits an invariant density in the appropriate $L^p$ space, weak convergence of the sequence of optimal solutions is observed; norm convergence can be obtained when the objective is a Kadec functional. No regularity or expansiveness assumptions on $T$ need to be made, and the method applies to maps on multidimensional domains. Objectives leading to norm convergence include Entropy, 'Energy' and 'Positively Constrained Energy'. Explicit solutions for the finite moment problems in the case of the 'Energy' functional are derived using duality - the optimality condition is then a linear algebra problem. Strong duality is obtained even though the dual functional may not be coercive and the set of moment test functions is not assumed to be pseudo-Haar. Finally, some numerical studies are presented for the case of moment test functions derived from a finite partition of the dynamical phase space and the results are compared with Ulam's method.