  The optimized linear $\delta$-expansion is applied to multi-field $O(N_1)
\times O(N_2)$ scalar theories at high temperatures. Using the imaginary time
formalism the thermal masses are evaluated perturbatively up to order
$\delta^2$ which considers consistently all two-loop contributions. A
variational procedure associated with the method generates nonperturbative
results which are used to search for parameters values for inverse symmetry
breaking (or symmetry nonrestoration) at high temperatures. Our results are
compared with the ones obtained by the one-loop perturbative approximation, the
gap equation solutions and the renormalization group approach, showing good
agreement with the latter method. Apart from strongly supporting inverse
symmetry breaking (or symmetry nonrestoration), our results reveal the
possibility of other high temperature symmetry breaking patterns for which the
last term in the breaking sequence is $O(N_1-1) \times O(N_2-1)$.
