Is the Lowdin Orthogonalization used in diagonalizing the atomic orbitals really a basis change?
Consider a nonsingular, linear transformation $\mathbf A$ of a set of vectors arranged into a matrix $\mathbf V$ $$\mathbf V'=\mathbf {VA}$$

We can say the matrix $\mathbf V'$ is orthogonal (and thus that we have linearly transformed the vectors to form an orthogonal basis) if $$\mathbf {V'^{\dagger}V'}=\mathbf {A^{\dagger}V^{\dagger}VA}=\mathbf {A^{\dagger}SA}=\mathbf{1}$$ 

In general, this occurs when $\mathbf{A}=\mathbf{S}^{-1/2}\mathbf{B}$, where $\mathbf{B}$ can be any unitary matrix. The Lowdin orthogonalization commonly used is just the case where $\mathbf{B=1}$, so we can see that we are properly changing the basis. 

But your issue is specifically with transformations of the density matrix, where the transformation matrix is not unitary in general. The standard change of basis formalism would suggest that to go to AO basis from MO, we would do a transformation something like this:
$$\mathbf{C}\mathbf{D}^{MO}\mathbf{C}^{-1}=\mathbf{D}^{AO}\mathbf{S}$$
Ensuring orthogonality of the MOs requires that $\mathbf{C^{\dagger}\mathbf{S}\mathbf{C}}=\mathbf{1}$, which we can sub in to the equation above to give: $$\mathbf{C}\mathbf{D}^{MO}\mathbf{C^{\dagger}\mathbf{S}\mathbf{C}}\mathbf{C}^{-1}=\mathbf{C}\mathbf{D}^{MO}\mathbf{C}^{\dagger}\mathbf{S} =\mathbf{D}^{AO}\mathbf{S} \to \mathbf{C}\mathbf{D}^{MO}\mathbf{C}^{\dagger}=\mathbf{D}^{AO}$$

So the form in terms $\mathbf C$ and its adjoint is a direct result of the conventional similarity transform combined with the orthonormality condition on the MOs. 

This latter derivation was adapted from$^1$ which can be found here.

As to why we avoid using transformations involving the inverse of the transformation matrix, I suspect Zhe is correct that it is a matter of cost, since obtaining the transpose/adjoint of a matrix requires almost no work at all, while an inverse is relatively difficult to obtain.


T. Helgaker et. al Chemical Physics Letters 327 (2000). 397–403