Why can we disentangle the squeezing operator $\exp[\frac12(\xi a^{\dagger 2}-\xi^* a^2)]$ via the $\mathfrak{su}(1,1)$ algebra? In the context of quantum mechanics, we define the squeezing operator $S(\xi)$ as:
$$S(\xi)\equiv \exp[\frac12(\xi a^{\dagger 2}-\xi^* a^2)],$$
where $a^\dagger$ and $a$ are the so-called creation and annihilation operators, which satisfy $[a,a^\dagger]=1$.
It is often convenient to consider the disentangled form of $S(\xi)$, that is, to split the operator into an exponential involving only $a^{\dagger 2}$ and one involving only $a^2$. This can be done as following:
$$S(r e^{i\theta})=
\exp\left[\frac12e^{i\theta}\tanh(r) a^{\dagger 2}\right]
\exp\left[-\ln\cosh(r) \left(a^\dagger a+\frac12\right)\right]
\exp\left[-\frac12e^{-i\theta}\tanh(r) a^{2}\right].
$$
This result is obtained via the following observations (see also here for more details):


*

*Define $K_+\equiv \frac12 a^{\dagger 2}$ and $K_-\equiv K_+^\dagger=\frac12 a^2$, and observe that
$$[K_-,K_+]=\frac12(1+2 a^\dagger a)\equiv K_0 \quad\text{ and }\quad
[K_0,K_\pm]=\pm2 K_\pm.$$

*Observe that these define the $\mathfrak{su}(1,1)$ Lie algebra, a faithful representation for which is given in terms of the Pauli matrices:
$$i\sigma_-\equiv\begin{pmatrix}0 & 0 \\ i & 0\end{pmatrix}=K_-,
\qquad i\sigma_+\equiv\begin{pmatrix}0 & i \\ 0 & 0\end{pmatrix}=K_+, \\
\sigma_3 = \begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}= K_0.$$

*Using this representation in terms of $2\times2$ matrices, we compute the matrix exponential $\exp[i(\xi \sigma_+ - \xi^* \sigma_-)]$ using standard techniques.

*We use Gauss decomposition to write the resulting matrix as product of matrices which correspond to exponentials of the original operators.


My question relates to the mathematical justification behind this procedure.
I understand that the representation of the operators in terms of $2\times2$ matrices is faithful, but that only relates to the commutators, the operators themselves don't behave the same way. So why should we expect to be able to compute the exponential in the representation and then convert the results back into the original operators?
 A: Consider the operators $$K_{+}=\frac{e^{i\theta}}{2}{a^{\dagger}} ^2~,~ K_{-}=\frac{e^{-i\theta}}{2}a^2~,~ K_0=a^{\dagger}a+\frac{1}{2}$$
These obey the exact same commutation relations quoted above. The algebra of these operators is linear and closes amongst the operators (hence is a Lie algebra). This is a really important fact for the considerations that follow below.
Consider now a one-dimensional flow along an arbitrary curve in the Lie group space that these generators span:
$$f(\lambda)=\exp(\alpha(\lambda)K_+)\exp(\beta(\lambda)K_0)\exp(\gamma(\lambda)K_-)$$ 
with $$\alpha(0)=\beta(0)=\gamma(0)=0$$
Denote the finite set of generators of the group as $\bar{G}=\{G_i, i=1,..., \dim(G)\}$.
Since the algebra closes and is linear, we can find functions $\Delta_{ijk}(s)$ such that
$$e^{sG_i}G_je^{-sG_i}=\sum_{k=1}^{dim(G)}\Delta_{ijk}(s)G_k$$
We specialize to the case considered here, even though the theorem holds with very general assumptions as outlined above. Take the derivative of $f(\lambda)$ and using properties of the matrix exponential defined above and deduce that:
$$\frac{df}{d\lambda}=\Big[\alpha'(\lambda)K_++\beta'(\lambda)(e^{\alpha K_+}K_0e^{-\alpha K_+})+\gamma'(\lambda)(e^{\alpha K_+}e^{\beta K_0}K_- e^{-\beta K_0}e^{-\alpha K_+})\Big]f(\lambda)$$
It is sufficient to note that the function in brackets $\rho(\lambda)$ is linear in the generators $\rho(\lambda)=a(\lambda)K_++b(\lambda)K_0+c(\lambda)K_-$. Then the general solution to this equation is given by the "time"-ordered exponential:
$$f(\lambda)=T_{\lambda}\exp\Big(\int_{0}^{\lambda}\rho(t)dt\Big)$$
I believe that it can be shown for finite-dimensional groups again under general assumptions and convergence properties, that there must exist functions such that
$$f(\lambda)=\exp(A(\lambda)K_++B(\lambda)K_0+C(\lambda)K_-)$$
This can be perhaps justified by looking at the time-ordered exponential as a product of infinitesimally small group flows, consider the BCH formula to put all the exponentials together and use the fact that all commutators in the BCH formula close under the generator algebra and are linear in the generators. Of course, there is no guarantee that BCH converges, so the above statement should be taken with a grain of salt.
The purpose of the above analysis however is not to prove, but rather motivate the fact that the disentangled form of the squeezing operator, should it exist, it must be INDEPENDENT of the representation, and only depending on the algebra of commutators.
Fortunately, to justify the squeezing operator decomposition, a general answer to the questions posed above is not necessary. It can be shown directly that defining the flow
$$g(\lambda)=\exp(-\mu(\lambda)K_+)\exp(\lambda(K_+-K_-))\exp(\mu(\lambda)K_-)$$
it's derivative can be directly evaluated similarly to the above sketch as
$$\frac{dg}{d\lambda}=\Bigg[[(1-\mu^2-\frac{d\mu}{d\lambda})+\frac{d\mu}{d\lambda}(\mu\cosh\lambda-\sinh\lambda)^2)]K_++(\cosh^2\lambda\frac{d\mu}{d\lambda}-1)K_-+(\cosh\lambda\frac{d\mu}{d\lambda}(\mu\cosh\lambda-\sinh\lambda)-\mu)K_0\Bigg]g(\lambda)$$ 
and miraculously, if we set the coefficient of $K_-$ to zero- which happens for $\mu(\lambda)=\tanh(\lambda)$- we very simply get
$$\frac{dg}{d\lambda}=-\mu(\lambda)K_0g(\lambda)$$
which we can trivially integrate for the quoted result, with evidently using only commutators to reach it:
$$g(\lambda)=\exp\Big(-K_0\int_{0}^{\lambda}dt~\mu(t)\Big)=\exp\Big(-K_0\ln\cosh\lambda\Big)$$
This has been a long answer, but the tl;dr is: 
1) Most statements of the form found in the OP can be generally computed by using commutators alone,
2) The above nice disentanglement formula corresponds to finding an integrable flow within the group. 
I hope this is helpful, but still I would love to see an answer that puts all the above into perspective.
