We present an algorithm for the optimization and thermal equilibration of spin glasses - or more generally, cost functions of the Ising form $H=sum_{langle i jrangle} J_{ij} s_i s_j + sum_i h_i s_i$, defined on graphs with arbitrary connectivity. The algorithm consists of two repeated steps: i) the optimized construction of a random tree of spin clusters on the input problem graph, and ii) the thermal sampling of the generated tree. The randomly generated trees are constructed so as to optimize a balance between the size of the tree and the complexity required to draw Boltzmann samples from it. We benchmark the algorithm on several classes of problems and demonstrate its advantages over existing approaches.