Glauber dynamics on spin systems are well known to suffer exponential slowdowns at low temperatures due to the emergence of multiple metastable phases, separated by narrow bottlenecks that are hard for the dynamics to cross. It is a folklore belief that if the dynamics is initialized from an appropriate random mixture of ground states, one for each phase, then convergence to the Gibbs distribution should be much faster. However, such phenomena have largely evaded rigorous analysis, as most tools in the study of Markov chain mixing times are tailored to worst-case initializations. In this paper we develop a general framework towards establishing this conjectured behavior for the Ising model. In the classical setting of the Ising model on an N-vertex torus in g.,Currency signd, our framework implies that the mixing time for the Glauber dynamics, initialized from a 1/2-1/2 mixture of the all-plus and all-minus configurations, is N1+o(1) in dimension d=2, and at most quasi-polynomial in all dimensions d≥ 3, at all temperatures below the critical one. The key innovation in our analysis is the introduction of the notion of "weak spatial mixing within a phase", a low-temperature adaptation of the classical concept of weak spatial mixing. We show both that this new notion is strong enough to control the mixing time from the above random initialization (by relating it to the mixing time with plus boundary condition at O(logN) scales), and that it holds at all low temperatures in all dimensions. This framework naturally extends to more general families of graphs. To illustrate this, we use the same approach to establish optimal O(NlogN) mixing for the Ising Glauber dynamics on random regular graphs at sufficiently low temperatures, when initialized from the same random mixture.