title: Matrix Equation solver Python
body: <p>
  Here's the problem: I am trying to solve a second order matricial equation of the form :
  <a href="https://i.stack.imgur.com/spY3e.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/spY3e.png" alt="enter image description here"></a></p>
  
  <p>Where X(to find) and C (known) are of sizes [nxn]. (n is on the order of 1000).
  C is a symmetric covariance matrix known. (and X should be symmetric too)</p>
  
  <p>
  Here is my code:</p>
  
  <pre><code>from sympy import solve
  from sympy import Indexed, IndexedBase, Tuple
  import numpy as np
  
  X = IndexedBase('X',shape=(n,n))
  eqs = Tuple(np.dot(X,X)-np.dot(C,X)-np.eye(n))
  solve(eqs, X)
  </code></pre>
  
  <p>Is it the right way to do it? My code is taking ages.
  I am looking for any type of algorithm that could help me solve this kind of equation efficiently. </p>