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_{p\alpha}^2+\frac{1}{2}\sum_{p,\alpha,i}\sum_{p^{\prime},\alpha^{\prime},i^{\prime}}D_{i\alpha;i^{\prime}\alpha^{\prime}}(\mathbf{R}_p,\mathbf{R}_{p^{\prime}})u_{pi\alpha}u_{p^{\prime}i^{\prime}\alpha^{\prime}}.
$$
The basic quantity one builds when calculating phonons is the matrix of force constants:
$$
D_{i\alpha;i^{\prime}\alpha^{\prime}}(\mathbf{R}_p,\mathbf{R}_{p^{\prime}})=\frac{\partial^2E}{\partial u_{pi\alpha}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_{pi\alpha}$ 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_{i\alpha;i^{\prime}\alpha^{\prime}}(\mathbf{q})=\frac{1}{N_p\sqrt{m_{\alpha}m_{\alpha^{\prime}}}}\sum_{\mathbf{R}_p,\mathbf{R}_{p^{\prime}}}D_{i\alpha;i^{\prime}\alpha^{\prime}}(\mathbf{R}_p,\mathbf{R}_{p^{\prime}})e^{i\mathbf{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;i\alpha}$. 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_{pi\alpha}e^{-i\mathbf{q}\cdot{\mathbf{R}_p}}v_{-\mathbf{q}\nu;i\alpha},
$$
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.
