def mds(d, dimensions=3, method='cpu'):
    (n, n) = d.shape
    E = (-0.5 * d ** 2)

    Er = np.mat(np.mean(E, 1))
    Es = np.mat(np.mean(E, 0))

    # From Principles of Multivariate Analysis: A User's Perspective (page 107).
    F = np.array(E - np.transpose(Er) - Es + np.mean(E))

    [U, S, V] = svd_cpu(F)

    Y = U * np.sqrt(S)

    return (Y[:, 0:dimensions], S)