title: Improve precision in solution for Sylvester equation with python (scipy)
body: <p>I am trying to solve numerically the following Sylvester equation (eq. (11) from <a href="https://arxiv.org/abs/1906.08808" rel="nofollow noreferrer">arXiv:1906.08808</a>)</p>
  
  <p><a href="https://i.stack.imgur.com/vruVg.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/vruVg.png" alt="Sylvester equation"></a></p>
  
  <p>where W+ is defined as exp(K*t). I have defined the following function to compute V(t) at a given time t.</p>
  
  <p><strong>Edited</strong> (after comment from <a href="https://stackoverflow.com/questions/59376329/improve-precision-in-solution-for-sylvester-equation-with-python-scipy#comment104950528_59376329">Warren Weckesser</a>)</p>
  
  <pre class="lang-py prettyprint-override"><code>import numpy as np
  from scipy.linalg import expm, solve_sylvester
  
  def covariance_matrix(K, V0, D, t, precision=precision):
      """
      It computes the covariance matrix at time t.
  
      See eq. (11) in arXiv:1906.08808
      """
      Wp = expm(K*t)
      WpT = Wp.transpose()
      KT = K.transpose()
      C = -D + K @ Wp @ V0 @ WpT + Wp @ V0 @ WpT @ KT + Wp @ D @ WpT
      return solve_sylvester(K, KT, C)
  </code></pre>
  
  <p>However, it is not working as expected. For example, if I call the function with t=0, I don't recover the initial V0 as it should (see screenshot below)</p>
  
  <p><a href="https://i.stack.imgur.com/QnF4M.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/QnF4M.png" alt="Screenshot Jupyter notebook"></a></p>
  
  <p>This looks like some kind of numerical error. I have tried different precisions (<code>np.float64</code> and <code>np.float128</code>) but I don't see any difference. My guess is that the precision is lost when using expm from <code>scipy.linalg</code>.</p>
  
  <p>Any ideas how to change this code to get better precision?</p>
  
  <p><strong>Expanded</strong> (also after edition)</p>
  
  <p>The problem with these small errors (order of 1e-13 at time 0) is that later I am using this matrix to compute some quantity that is very sensible to this little differences.</p>
  
  <p>I leave here the exact code that I am using later with this matrix.</p>
  
  <pre class="lang-py prettyprint-override"><code>def entanglement(V):
      """
      It gives the logarithmic negativity computed from the covariance matrix.
      See: Phys. Rev. A 65, 032314 (2002) for details.
      """
      Ia = V[:2, :2]
      Ib = V[2:, 2:]
      L  = V[:2, 2:]
      s = det(Ia) + det(Ib) - 2 * det(L)
      v = np.sqrt((s - np.sqrt(s**2 - 4 * det(V))) / 2)
      ent = -np.log2(2*v)
      if ent &gt; 0:
          return ent
      else:
          return 0
  </code></pre>
  
  <p><strong>Minimal, reproducible example</strong></p>
  
  <pre class="lang-py prettyprint-override"><code>#!/usr/bin/python3
  import numpy as np
  from scipy.linalg import expm, det, solve_sylvester
  
  precision = np.float128
  
  def drift_matrix(w, eta, Q):
      return np.array([[0,           w,    0,        0],
                       [-w*(1-eta), -w/Q, -w*eta,      0],
                       [0,           0,    0,        w],
                       [-w*eta,      0, -w*(1-eta), -w/Q]], dtype=precision)
  
  def init_covariant_matrix(n, s=0):
      return np.diag(np.ones(4, dtype=precision) * (n + 1/2) * np.cosh(2*s))
  
  def diagonal_matrix(n, w, Q):
      return np.diag(np.array([0,1,0,1], dtype=precision) * w/Q * (1 + 2*n))
  
  def covariance_matrix(K, V0, D, t):
      Wp = expm(K*t)
      WpT = Wp.transpose()
      KT = K.transpose()
      C = -D + K @ Wp @ V0 @ WpT + Wp @ V0 @ WpT @ KT + Wp @ D @ WpT
      return solve_sylvester(K, KT, C)
  
  w = precision(0.1)
  eta = precision(1.36e-6) / w**2
  n = precision(0)
  Q = precision(1e4)
  K = drift_matrix(w, eta, Q)
  V0 = init_covariant_matrix(n, 1.73)
  D = diagonal_matrix(n, w, Q)
  V0_1 = covariance_matrix(K, V0, D, 0)
  
  # What I want (ideally)
  np.all(np.abs(V0 - V0_1) &lt; np.finfo(np.float128).resolution)
  
  # I would still be happy if...
  np.all(np.abs(V0 - V0_1) &lt; np.finfo(np.float64).resolution)
  </code></pre>