docs: def semidefinite_relaxation(edm_missing, lamda, W=None, print_out=False, **kwargs):
      """Complete and denoise EDM using semidefinite relaxation.
  
      Returns solution to the relaxation of the following problem: 
          D = argmin || W * (D - edm_missing) || 
              s.t. D is EDM
      where edm_missing is measured matrix, W is a weight matrix, and * is pointwise multiplication. 
  
      Refer to paper "Euclidean Distance Matrices - Essential Theory, Algorithms and Applications", 
      Algorithm 5, for details. (https://www.doi.org/%2010.1109/MSP.2015.2398954)
  
      :param edm_missing: EDM with 0 where no measurement was taken. 
      :param lamda: Regularization parameter.  
      :param W: Optional mask. If no mask is given, a binary mask is created based on missing elements of edm_missing. 
                If mask is given 
      :param kwargs: more options passed to the solver. See cvxpy documentation for all options. 
      """ 
      from .algorithms import reconstruct_mds
  
      def kappa(gram):
          n = len(gram)
          e = np.ones(n)
          return np.outer(np.diag(gram), e) + np.outer(e, np.diag(gram).T) - 2 * gram
  
      def kappa_cvx(gram, n):
          e = np.ones((n, 1))
          return reshape(diag(gram), (n, 1)) * e.T + e * reshape(diag(gram), (1, n)) - 2 * gram
  
      method = kwargs.pop('method', 'maximize')
      options = {'solver': 'CVXOPT'}
      options.update(kwargs)
  
      if W is None:
          W = (edm_missing > 0)
      else:
          W[edm_missing == 0] = 0.0
  
      n = edm_missing.shape[0]
      V = np.c_[-np.ones((n - 1, 1)) / np.sqrt(n), np.eye(n - 1) -
                np.ones((n - 1, n - 1)) / (n + np.sqrt(n))].T
  
      H = Variable((n - 1, n - 1), PSD=True)
      G = V * H * V.T  # * is overloaded
      edm_optimize = kappa_cvx(G, n)
  
      if method == 'maximize':
          obj = Maximize(trace(H) - lamda *
                         norm(multiply(W, (edm_optimize - edm_missing)), p=1))
      # TODO: add a reference to paper where "minimize" is used instead of maximize. 
      elif method == 'minimize':
          obj = Minimize(trace(H) + lamda *
                         norm(multiply(W, (edm_optimize - edm_missing)), p=1))
  
      prob = Problem(obj)
  
      total = prob.solve(**options)
      if print_out:
          print('total cost:', total)
          print('SDP status:', prob.status)
  
      if H.value is not None:
          Gbest = V.dot(H.value).dot(V.T)
          if print_out:
              print('eigenvalues:', np.sum(np.linalg.eigvals(Gbest)[2:]))
          edm_complete = kappa(Gbest)
      else:
          edm_complete = edm_missing
  
      if (print_out):
          if H.value is not None:
              print('trace of H:', np.trace(H.value))
          print('other cost:', lamda *
                norm(multiply(W, (edm_complete - edm_missing)), p=1).value)
  
      return np.array(edm_complete)