For the example given, there are four equations, and 16 variables. So you will need to constrain the problem more if you want a unique solution.

If you treated it as a transformation using a 2 x 2 matrix, it could be done like so.

amat = Array[a, {2, 2}];
romat = Array[rho, {2, 2}, 0];
Solve[amat . 
   romat == {{rho[1, 1], rho[1, 0]}, {rho[0, 1], rho[0, 0]}}, 
 Flatten[amat]]

(* Out[11]= {{a[1, 1] -> -((-rho[1, 0]^2 + rho[1, 1]^2)/(
    rho[0, 1] rho[1, 0] - rho[0, 0] rho[1, 1])), 
  a[1, 2] -> -((rho[0, 0] rho[1, 0] - rho[0, 1] rho[1, 1])/(
    rho[0, 1] rho[1, 0] - rho[0, 0] rho[1, 1])), 
  a[2, 1] -> -((-rho[0, 0] rho[1, 0] + rho[0, 1] rho[1, 1])/(
    rho[0, 1] rho[1, 0] - rho[0, 0] rho[1, 1])), 
  a[2, 2] -> -((rho[0, 0]^2 - rho[0, 1]^2)/(
    rho[0, 1] rho[1, 0] - rho[0, 0] rho[1, 1]))}}*)