  In this work, we study the circuit complexity for generalized coherent states
in thermal systems by adopting the covariance matrix approach. We focus on the
coherent thermal (CT) state, which is non-Gaussian and has a nonvanishing
one-point function. We find that even though the CT state cannot be fully
determined by the symmetric two-point function, the circuit complexity can
still be computed in the framework of the covariance matrix formalism by
properly enlarging the covariance matrix. Now the group generated by the
unitary is the semiproduct of translation and the symplectic group. If the
reference state is Gaussian, the optimal geodesic is still be generated by a
horizontal generator such that the circuit complexity can be read from the
generalized covariance matrix associated to the target state by taking the cost
function to be $F_2$. For a single harmonic oscillator, we discuss carefully
the complexity and its formation in the cases that the reference states are
Gaussian and the target space is excited by a single mode or double modes. We
show that the study can be extended to the free scalar field theory.
