Variational methods are a common approach for computing properties of ground states but have not yet found analogous success in finite temperature calculations. In this work we develop a new variational finite temperature algorithm (VAFT) which combines ideas from minimally entangled typical thermal states (METTS), variational Monte Carlo (VMC) optimization and path integral Monte Carlo (PIMC). This allows us to define an implicit variational density matrix to estimate finite temperature properties in two and three dimensions. We benchmark the algorithm on the bipartite Heisenberg model and compare to exact results.Strongly correlated fermionic and frustrated spin systems span many interesting physical systems. Computing properties of these systems is difficult, particularly in dimensions greater then one where density matrix renormalization group (DMRG) Of course, it is also possible to use the variational technique at finite temperature. Most naturally one might parameterize the finite temperature many body density, optimizing again over the parameters α. However, this requires specifying a set of FTDM over which to optimize, a task that is naturally more difficult than specifying a set of wavefunctions. Additionally, this optimization needs to minimize the free energy of the system, a quantity which is harder to evaluate, requiring techniques such as coupling constant integration. For one example of this approach, see ref. [19].Producing a good variational ansatz for very high or low temperature density matrices is straightforward.At high temperature (small β), one can approximate exp[−βĤ] by either Taylor expanding as (1 − βĤ) or using the Trotter break-up. At low temperature, one can instead expand in terms of the low-energy excitations, writingρ =and E i are variational estimates for the ith eigenvector and eigenvalue respectively and k is some cutoff. Unfortunately, this can be difficult, requiring orthogonal variational estimates of a number of excited states, and as one goes to higher temperature the number of eigenstates increases exponentially. Nonetheless, this framework is often applied in the context of density functional theory (DFT) to produce finite temperature functionals [20].At any temperature, one can approximate the density matrix as a purificationρ ≈ i λ i |Ψ i Ψ i | over a set |Ψ i of variational wavefunctions. However, one cannot efficiently enumerate more than a polynomial number of |Ψ i , nor is it obvious how to optimize over these wavefunctions.In spite of these difficulties, it would be useful if we could use known variational wave-functions to compute finite temperature properties. In this paper, we describe a new variational finite temperature approach (VAFT), in which we implicitly define a variational finite temperature density matrix. Our variational FTDM is a purification of an exponentially large number of variational wavefunctions. We compute properties of this FTDM via an algorithm which stochastically samples from its diagonal without ever needing to explicitly represent...