Solving $A^*A=B$ There are two 4×4 matrices $A$ and $B$ with $B$ symmetrical and
$$A^*A=B$$
Given that $A$ is completely unknown and that $B$ is known but is a function of a variable $k$, I want to solve for each element of $A$ as a function of $k$.
Progress so far
I believe that given the equation, we have 10 equations and 10 unknowns, meaning $A$ is also symmetrical. 
I first tried without the conjugate. To do so, I tried using Matlab and Maple solve functions to process this system of equations but they either never finish or give up with an error warning that I cannot make much use of.
I also tried using conjugate transpose on Matlab, but it doesn't enjoy the conjugate of a symbol.
I have the feeling I am doing something very wrong here. Are there perhaps any mathematical properties I can leverage to solve this?
 A: The conjugate transpose of $A$ is just the adjoint, commonly denoted $A^*$. So you're trying to solve $$(A^*)A=B.$$
Yes, this requires that $B^*=B$. If in fact $B^*=B$ you can start by diagonalizing $B$. You get $$B=PDP^{-1}$$where $D$ is diagonal and $P$ is unitary. Now if $E$ is another diagonal matrix and $A=PEP^{-1}=PEP^*$ then $A^*=PE^*P^*$, so $$A^*A=PE^*EP^{-1};$$you just need to find a diagonal matrix $E$ with $E^*E=D$, which is easy.
A: You can't do what you want unless $B$ has only non-negative eigenvalues because $A^*A$ has only non-negative eigenvalues, assuming $A^*$ denotes conjugate transpose. In that case, you can choose $A$ to be the unique positive square root $\sqrt{B}$ of $B$, so that $A=A^*$. That results in a well-posed system with a unique solution. If you normalize $B$ so that the eigenvalues are in $[0,1]$, then you might try a binomial series expansion for $(1-x)^{1/2}$:
$$
         A=\sqrt{B}=(I-(I-B))^{1/2}=I-\sum_{n=1}^{\infty}c_n (I-B)^n.
$$
The coefficients $c_n$ are positive and sum to $1$, which makes the convergence monotone in an operator sense. The convergence rate is determined by how far the non-zero eigenvalues of $B$ are from $0$. This method should be numerically stable, and will result in a positive $A$. If $B=I$, the sum is $0$ and $A=I$ is the correct solution.
A: A few things that aren't a complete answer, but are too long for a comment:


*

*$A$ need not be symmetric, and it's actually quite non-unique, since if $U$ is unitary, then $(UA)^{*}(UA)=A^{*}U^{*}UA=A^{*}A=B$ (I'm using $A^{*}$ to denote the conjugate transpose of a matrix $A$).

*The Cholesky factorization gives such a decomposition, and may be suitable depending on the matrix $B$. Note that you may end up with some algebraic functions of $k$ for the entries of $A$ (square roots of rational functions in $k$, assuming the entries of $B$ are polynomials in $k$) using this approach.

*Another approach is to find the eigendecomposition of $B=UDU^{*},$ where $U$ is unitary and $D$ has real, nonnegative entries (otherwise the representation $B=A^{*}A$ is impossible). Then $A=D^{1/2}U^{*}$ works, though again this may have algebraic functions of $k$ for the entries of $A$.

*This may be clear already, but unless you have some extra structure on $B$, even if its entries are polynomials in $k,$ the entries of $A$ will not be.
