The Explanation from @ScienceGeyser provides a good explanation to the phenomenon.
There are two more things to address the question on how is this phenomenon avoided.
The feedback read by the microphone is not identical to the audio sent to the speakers. There is the physical response from the speakers, the acoustic of the device and the environment, if microphone and speakers fixed as parts of the same rigid device like in a cell phone or a headset.
For a situation where there is wires or wireless connection between speaker and microphone, the position from one in relation to the other may change.
In any case the system can be modeled as follows, $k(z)$ is the speaker signal, $m(z)$ is the microphone signal, and $s(z)$ is the response from the sources other than the speaker to the microphone. We can say assume that the response of the microphone to the speaker sound is $F(z)\, k(z)$ in which $F(z)$ is a slowly varying system, so if we process short segments of the signal we may assume $F(z)$ to be time invariant.
Then in order to successfully cancel the response from the system you have to compute $F(z)$.
Particular case
To simplify the notation of the formulas in the frequency domain we $z = exp(j \omega T_s)$ where $T_s$ is the sampling period, and $\omega$ some angular frequency, under these conditions that discrete Fourier transform coincides with z transform .
If you assume $F(z)$ to have a short input response and to change slowly enough that you can use it as a time invariant filter in a block of $N$ samples.
We have $m(z) = F(z) k(z) + s(z)$, in your device you have both $k(z)$ and $m(z)$ (in the time domain).
You can compute the energy in the frequency domain as $m(z) m(z)^*$, I am using the $*$ suffix to denote complex conjugation.
$$ \begin{eqnarray}
m(z)^*m(z) &=& (F(z) k(z) + s(z))^*(F(z) k(z) + s(z)) \\
&=& \underbrace{|F(z)k(z)|^2}_{\textrm{echo energy}} + 2 Re(F(z)k(z)s(z)^*) + \underbrace{|s(z)|^2}_{\textrm{external sound}}
\end{eqnarray} $$
It is reasonable to assume that $k(z)$ and $s(z)$ must be uncorrelated, and one can see that the echo will increase the energy of the signal from the microphone.
The corrected signal will be given by
$$\hat s(z) = m(z) - \hat{F}(z) k(z) = s(z) + k(z)(F(z) - \hat{F}(z))$$
And $\hat{F}(z)$ is chosen in a way that minimizes the energy of $\hat s(z)$, that in the frequency domain consists in a linear regression problem. Let $\hat F_i$, $m_i$ and $k_i$ be the coefficients of $\hat F(z)$, $m(z)$ and $k(z)$, respectively
$$ \hat F_i = \frac{m_i k_i^*}{\sum |k_i|^2} $$
Using big number of samples will give better estimate of $\hat F(z)$ but the time invariability assumption is more strongly violated. What can be done is to take advantage of the previous estimates. To address this we can filter $\hat F(z)$ itself.
$$\tilde{F}_i \leftarrow \tilde{F}_i \alpha + (1 - \alpha) \frac{m_i k_i^*}{\sum |k_i|^2}$$
And you can increase alpha for a more stable estimate, or decrease it to have a faster estimate.
Notes on complex arithmetic
To answer the question about the term 2 Re(.) left on the comment. The short answer is this is the crossed terms in the expansion of the product in the previous.
I skipped some steps there to avoid making long arrays of equations since this is a very common result in complex manipulations.
$$\begin{eqnarray}
(x + y)^* (x + y) &=& (x^* + y^*)(x + y) \\
&=& (x^* x + y^* y + x^*y + y^*x \\
&=& x^* x + y^*y + (x^*y) + (x^*y)^*
\end{eqnarray}$$
Let $x = a+bj$, with $j^2=-1$, term $x^*x=|x|^2$,
$$\begin{eqnarray}x^*x &=& (a-bj)(a+bj) \\&=& (a^2 - b^2j^2 -baj +abj) \\&=& (a^2 -b^2(-1)) \\&=& (a^2 + b^2) \\&=& |x|^2\end{eqnarray}$$
Similarly, $y^*y = |y|^2$.
The terms $x^*y + y^*x$ are conjugate from each other, thus they could be written as $w^* + w$, since $w = Re(w) + j Im(w)$ and $w^* = Re(w) - j Im(w)$, $w + w^* = 2 Re(w)$, where $Re(w)$ denotes the real part, and $Im(w)$ denotes te imaginary part.
Why the energy is the square of the amplitude
I will explain this using a Mass-spring system, let a mass $m$ attached to a sprint of constant $k$ oscillating with amplitude $A$. The maximum potential energy of the spring is achieved when it is deformed by $A$, and can be calculated as $\int_{0}^{A} k \, x\, dx = k A^2 / 2$. Also, if there are multiple oscilation modes (different frequencies), each one conserves energy separately, this is why we can talk about power at specific frequencies.
