How to deduce phase transitions from a phonon calculation?

The answer of @ProfM is already very complete, but I wanted to tackle your question from a more practical point of view.

The presence of imaginary frequencies indicate that there are atomic positions which are more energetically favorable at the ground state. So, the concept of "following" a mode means condensing it onto the reference structure, until you find the equilibrium positions.

To give an example, we can start with $ce{BaTiO_3}$ with $Pmoverline{3}m$ symmetry. By looking at the phonon frequencies you will notice the instability (imaginary frequency) at $Gamma$, which for this case in particular corresponds to the so-called ferroelectric mode.

Once you know the phonon of interest, you can read its eigendisplacements ($U_{FE}$), condense them onto the reference structure ($S_{ref}$) using different amplitudes $alpha$, $$S_{alpha} = S_{ref} + alpha U_{FE}$$ and then calculate the energy for each of the resulting structures ($S_{alpha}$). You will notice that, as you increase the value $|alpha|$, the total energy of the system decreases until it starts increasing again, creating a potential well. The structure with the minimum energy will be the one at equilibrium (according to the selected mode).

Coming back to the $ce{BaTiO_3}$ example, the new structure $S_{alpha, E_{min}}$ should show a $P4mm$ symmetry. However, you would need to relax the phase, as strain also needs to be included.

Note: $alpha$ can be any real value but it is case-specific (depending on the anharmonicities of the system and on the definition of the eigendisplacements), so would have to try with different options to determine which one is better for your case.

Background theory. In the harmonic approximation, the potential energy surface (PES) is expanded about an equilibrium point to second order, to obtain the Hamiltonian:

$$ hat{H}=sum_{p,alpha}-frac{1}{2m_{alpha}}nabla_{palpha}^2+frac{1}{2}sum_{p,alpha,i}sum_{p^{prime},alpha^{prime},i^{prime}}D_{ialpha;i^{prime}alpha^{prime}}(mathbf{R}_p,mathbf{R}_{p^{prime}})u_{pialpha}u_{p^{prime}i^{prime}alpha^{prime}}. $$

The basic quantity one builds when calculating phonons is the matrix of force constants:

$$ D_{ialpha;i^{prime}alpha^{prime}}(mathbf{R}_p,mathbf{R}_{p^{prime}})=frac{partial^2E}{partial u_{pialpha}u_{p^{prime}i^{prime}alpha^{prime}}}, $$

which is the expansion coefficient of the second order term in the potential energy surface $E$, with $i$ labelling the Cartesian direction, $alpha$ the atom in the basis, $mathbf{R}_p$ is the position of the cell $p$ in the crystal, and $u_{pialpha}$ is the amplitude of displacement of the corresponding atom. Using periodicity of the crystal, we can define the dynamical matrix at each $mathbf{q}$-point of the Brillouin zone as:

$$ D_{ialpha;i^{prime}alpha^{prime}}(mathbf{q})=frac{1}{N_psqrt{m_{alpha}m_{alpha^{prime}}}}sum_{mathbf{R}_p,mathbf{R}_{p^{prime}}}D_{ialpha;i^{prime}alpha^{prime}}(mathbf{R}_p,mathbf{R}_{p^{prime}})e^{imathbf{q}cdot(mathbf{R}_p-mathbf{R}_{p^{prime}})}, $$

where $N_p$ is the number of cells in the supercell over which periodic boundary conditions are applied, and $m_{alpha}$ is the mass of atom $alpha$. Diagonalizing the dynamical matrix gives eigenvalues $omega^2_{mathbf{q}nu}$ and eigenvectors $v_{mathbf{q}nu;ialpha}$. From these, it is possible to define a set of normal coordinates:

$$ u_{mathbf{q}nu}=frac{1}{sqrt{N_p}}sum_{mathbf{R}_p,i,alpha}sqrt{m_{alpha}}u_{pialpha}e^{-imathbf{q}cdot{mathbf{R}_p}}v_{-mathbf{q}nu;ialpha}, $$

in terms of which the Hamiltonian becomes a sum over uncoupled simple harmonic oscillators:

$$ hat{H}=sum_{mathbf{q},nu}-frac{1}{2}frac{partial^2}{partial u_{mathbf{q}nu}^2}+frac{1}{2}omega^2_{mathbf{q}nu}u_{mathbf{q}nu}^2. $$

The bosonic quasiparticles labelled by quantum numbers $(mathbf{q},nu)$ are called phonons, and have energy $omega_{mathbf{q}nu}$ and momentum $mathbf{q}$.

Dynamically stable structure. A dynamically stable structure is one whose equilibrium position is at a local minimum of the PES. As such, the eigenvalues $omega^2_{mathbf{q}nu}$ of the dynamical matrix (Hessian) are all positive numbers, and as a consequence the phonon frequencies $omega_{mathbf{q}nu}$ are all real.

Dynamically unstable structure. A dynamically unstable structure is one whose equilibrium position is at a saddle point of the PES. As such, some of the eigenvalues of the dynamical matrix are negative, and the corresponding phonon frequencies imaginary.

Physical interpretation. Phonons measure the curvature of the PES around the equilibrium position of the material. As we have seen, an imaginary frequency corresponds to a negative curvature, so it corresponds to a direction in the PES along which the energy decreases. This means that there is a lower energy configuration of the material, and we say that the structure is then dynamically unstable.

"Follow the imaginary modes". How can we find such lower-energy structure? The eigenvectors of the dynamical matrix associated with the imaginary phonons tell us the direction along which the energy decreases, so we can "follow those modes" to find the lower energy structure. This can be done by simply constructing a sequence of structures on which you displace the atoms by an amplitude $u_{mathbf{q}nu}$ (see equation above) of the imaginary phonon $(mathbf{q},nu)$, and calculating the total energy of each of the resulting structures. For a saddle point, the resulting curve will be something like a double well, and the minima of the double well correspond to your new lower-energy structure.

Finite temperature. The discussion up to this point concerns the potential energy surface, so temperature is neglected. If you are interested in a calculation at finite temperature, then you need the free energy surface. This is much harder to calculate and you need anharmonic terms in your Hamiltonian to describe it properly.

Perovskites. Perovskites typically have a cubic structure at high temperature, and then upon lowering the temperature undergo a number of phase transitions to lower-symmetry structures (tetragonal, orthorhombic, and so on). Imagine a perovskite that has only two phases, tetragonal at low temperature and cubic at high temperature (generalizing to more phases is trivial). Then if you calculate the phonons in the cubic structure (saddle point) you will find imaginary modes, and following them will take you to the tetragonal structure (minimum). If you calculate the phonons in the tetragonal structure, they will all have real frequencies. So why is the cubic phase stable at high temperatures? This is because, although the cubic phase corresponds to a saddle point of the potential energy surface, above some critical temperature it corresponds to a minimum of the free energy surface. As such, above that critical temperature the cubic phase becomes dynamically stable. As I mentioned above, to investigate this phase transition (e.g. to calculate the critical temperature), you need to include anharmonic terms (phonon-phonon interactions), which is much harder computationally.