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)