Suppose I have three mutually independent non-Gaussian noise $E_A$, $E_B$, $E_C$. There are two variables generated by linear combinations of these noise components: $M=pE_A+qE_B$, $N=rE_B$. By linear regress $M$ on $N$ (with least square error), we will get regression residual ($qE_B$) independent to $N$, so we know the common noise components shared by variables $M$ and $N$ is exactly the noise components that $N$ has ($E_B$).
However, consider another case where $M=pE_A+qE_B$, $N=rE_B+sE_C$. Now for either linear regression of $M$ on $N$ or $N$ on $M$, the residual will not be independent to the regressor. So in this case, is there anyway that one could separate out the common noise that $M$ and $N$ shares (i.e. to obtain a variable that has only noise component $E_B$)? I expect the method to be as simple as possible, and preventing from using ICA.
I'm a little confused. I don't have enough reputation for commenting, so I decide to write a partial answer.
We have mutually independent random variables (noises) $E_A,E_B$, and $E_C$. I suppose you mean that $M=pE_A+qE_B$ and $N = rE_B$ are both observable, and that $p,q,r$ are some fixed but unknown scalar parameters. Say we collect $n$ independent observations, denoted $$\mathbf{y}=(M_1,...,M_n)'=(pe_{A1}+qe_{B1},...,pe_{An}+qe_{Bn})'=p\mathbf{e}_A+q\mathbf{e}_B,$$ and $$\mathbf{x}=(N_1,...,N_n)=(re_{B1},...,re_{Bn})'=r\mathbf{e}_B.$$ The least square method is to orthogonally project $\mathbf{y}$ onto $\mathcal{C}(\mathbf{x})$, the column space of $\mathbf{x}$.
Formally, we write a model $\mathbf{y}=\mathbf{x}\theta+{\epsilon}$, where $\theta$ is a scalar-valued parameter, $\mathbf{x}\theta$ is regarded as $\mathbb{E}(\mathbf{y}|\mathbf{x})$, and anything remaining is absorbed into the error vector $\mathbf{\epsilon}$. Assuming $\mathbf{x}=r\mathbf{e}_B$ is not a null vecor, we have the orthogonal projector onto $\mathcal{C}(\mathbf{x})$
$$\mathbf{P_x}\overset{\mathrm{def}}{=}\mathbf{x}(\mathbf{x}'\mathbf{x})^{-1}\mathbf{x}'=\frac{\mathbf{x}\mathbf{x}'}{\mathbf{x}'\mathbf{x}}=\frac{(r\mathbf{e}_B)(r\mathbf{e}_B)'}{(r\mathbf{e}_B)'(r\mathbf{e}_B)}=\frac{\mathbf{e}_B\mathbf{e}_B'}{\mathbf{e}_B'\mathbf{e}_B}\overset{\mathrm{def}}{=}\mathbf{P}_{\mathbf{e}_B}.$$
Hence, the residual, which is obtained by projecting $\mathbf{y}$ onto the orthogonal space of $\mathcal{C}(\mathbf{x})$, is $$\hat{\epsilon}=(\mathbf{I}-\mathbf{P_x})\mathbf{y}=(\mathbf{I}-\mathbf{P}_{\mathbf{e}_B})\mathbf{y} = (\mathbf{I}-\mathbf{P}_{\mathbf{e}_B})(p\mathbf{e}_A+q\mathbf{e}_B)=(\mathbf{I}-\mathbf{P}_{\mathbf{e}_B})p\mathbf{e}_A,$$ because $(\mathbf{I}-\mathbf{P}_{\mathbf{e}_B})\mathbf{e}_B=\mathbf{0}$, where $\mathbf{I}$ is the identity matrix. The residual, viewed as an observed vector, is linear independent to the observed regressor $\mathbf{x}=(N_1,...,N_n)'$, but not independent in the probability/distribution sense, because of the projection matrix $\mathbf{P}_{\mathbf{e}_B}$. And the fitted value is $$\hat{\mathbf{y}}=\mathbf{P_x}\mathbf{y}=\mathbf{P}_{\mathbf{e}_B}(p\mathbf{e}_A+q\mathbf{e}_B)=p\mathbf{P}_{\mathbf{e}_B}\mathbf{e}_A+q\mathbf{e}_B.$$
I don't see where $qE_B$ is coming from in the residual given by you.
In your second example, the residual is $$(\mathbf{I}-\mathbf{P}_{\mathbf{e}_B+\mathbf{e}_C})\mathbf{y}=(\mathbf{I}-\mathbf{P}_{\mathbf{e}_B+\mathbf{e}_C})(p\mathbf{e}_A+q\mathbf{e}_B).$$ There is nothing we can do here, because $\mathbf{e}_B+\mathbf{e}_C$ can be any vector of length $n$ and so $\mathbf{e}_B$ is not necessarily in $\mathcal{C}(\mathbf{e}_B+\mathbf{e}_C)$. Your best bet is to regress $N$ on something related to $\mathbf{e}_C$ only, and then the residual would be related to $\mathbf{e}_B$ only; and then regress $M$ on the residual to get another residual that related with $\mathbf{e}_A$ only.