docs: def svd(a, compute_uv=True, rcond=None):
      """ svd decomposition of matrix ``a`` containing |GVar|\s.
  
      Args:
          a: Two-dimensional matrix/array of numbers
              and/or :class:`gvar.GVar`\s.
          compute_uv (bool): It ``True`` (default), returns
              tuple ``(u,s,vT)`` where matrix ``a = u @ np.diag(s) @ vT``
              where matrices ``u`` and ``vT`` satisfy ``u.T @ u = 1``
              and ``vT @ vT.T = 1``, and ``s`` is the list of singular
              values. Only ``s`` is returned if ``compute_uv=False``.
          rcond (float): Singular values whose difference is smaller than
              ``rcond`` times their sum are assumed to be degenerate for
              calculating variances for ``u`` and ``vT``.
              Default (``rcond=None``) is ``max(M,N)`` times machine precision.
  
      Returns:
          Tuple ``(u,s,vT)`` where matrix ``a = u @ np.diag(s) @ vT``
          where matrices ``u`` and ``vT`` satisfy ``u.T @ u = 1``
          and ``vT @ vT.T = 1``, and ``s`` is the list of singular
          values. If ``a.shape=(N,M)``, then ``u.shape=(N,K)``
          and ``vT.shape=(K,M)`` where ``K`` is the number of
          nonzero singular values (``len(s)==K``).
          If ``compute_uv==False`` only ``s`` is returned.
  
      Raises:
          ValueError: If matrix is not two-dimensional.
      """
      a = numpy.asarray(a)
      if a.dtype != object:
          return numpy.linalg.svd(a, compute_uv=compute_uv)
      amean = gvar.mean(a)
      if amean.ndim != 2:
          raise ValueError(
              'matrix must have dimension 2: actual shape = ' + str(a.shape)
              )
      if rcond is None:
          rcond = numpy.finfo(float).eps * max(a.shape)
      da = a - amean
      u0,s0,v0T = numpy.linalg.svd(amean, compute_uv=True, full_matrices=True)
      k = min(a.shape)
      s = s0 + [
          u0[:, i].dot(da.dot(v0T[i, :])) for i in range(k)
          ]
      if compute_uv:
          u = numpy.array(u0, dtype=object)
          vT = numpy.array(v0T, dtype=object)
          # u first
          daaT = da.dot(a.T) + a.dot(da.T)
          s02 = numpy.zeros(daaT.shape[0], float)
          s02[:len(s0)] = s0 ** 2
          for j in range(s02.shape[0]):
              for i in range(k):
                  if i == j:
                      continue
                  ds2 = s02[i]  - s02[j]
                  if abs(ds2) < rcond * abs(s02[i] + s02[j]) or ds2 == 0:
                      continue
                  u[:, i] +=  u0[:, j]  * u0[:, j].dot(daaT.dot(u0[:, i])) / ds2
          # v next
          daTa = da.T.dot(a) + a.T.dot(da)
          s02 = numpy.zeros(daTa.shape[0], float)
          s02[:len(s0)] = s0 ** 2
          for j in range(s02.shape[0]):
              for i in range(k):
                  if i == j:
                      continue
                  ds2 = s02[i]  - s02[j]
                  if abs(ds2) < rcond * abs(s02[i] + s02[j]) or ds2 == 0:
                      continue
                  vT[i, :] +=  v0T[j, :]  * v0T[j, :].dot(daTa.dot(v0T[i, :])) / ds2
          return u[:,:k], s, vT[:k, :]
      else:
          return s