docs: def opt_space(M_E, r=None, niter=50, tol=1e-6, print_out=False):
      '''
      Implementation of the OptSpace matrix completion algorithm.
      An algorithm for Matrix Reconstruction from a partially revealed set.
      Sparse treatment of matrices are removed because of indexing problems in Python.
      Args:
          M_E:    2D numpy array; The partially revealed matrix.
                  Matrix with zeroes at the unrevealed indices.
          r:      The rank to be used for reconstruction. If left empty, the rank is guessed in the program.
          niter:  Maximum number of iterations.
          tol:    Stop iterations if norm((XSY' - M_E) * E, 'fro') / sqrt(|E|) < tol, where
                  E_{ij} = 1 if M_{ij} is revealed and zero otherwise,
                  |E| is the size of the revealed set.
      Returns: The following
          X:      A M_E.shape[0]xr numpy array
          S:      An rxr numpy array
          Y:      A M_E.shape[1]xr numpy matrix such that M_hat = X*S*Y'
          errs:   A vector containing norm((XSY' - M_E) * E, 'fro') / sqrt(|E|) at each iteration.
      '''
      n, m = M_E.shape
      # construct the revealed set
      E = np.zeros(M_E.shape)
      E[np.nonzero(M_E)] = 1
      eps = np.count_nonzero(E) / np.sqrt(m * n)
  
      if r is None:
          print('Rank not specified. Trying to guess ...')
          r = guess_rank(M_E)
          print('Using Rank : %d' % r)
  
      m0 = 10000
      rho = 0
  
      rescal_param = np.sqrt(np.count_nonzero(
          E) * r / np.linalg.norm(M_E, ord='fro')**2)
      M_E = M_E * rescal_param
  
      if print_out:
          print('Trimming ...')
  
      M_Et = copy.deepcopy(M_E)
      d = E.sum(0)
      d_ = np.mean(d)
  
      for col in range(m):
          if E[:, col].sum() > 2 * d_:
              nonzero_ind_list = np.nonzero(E[:, col])
              p = np.random.permutation(len(nonzero_ind_list))
              M_Et[nonzero_ind_list[p[np.ceil(2 * d_):]], col] = 0
  
      d = E.sum(1)
      d_ = np.mean(d)
      for row in range(n):
          if E[row, :].sum() > 2 * d_:
              nonzero_ind_list = np.nonzero(E[row, :])
              p = np.random.permutation(len(nonzero_ind_list))
              M_Et[nonzero_ind_list[row, p[np.ceil(2 * d_):]]] = 0
  
      if print_out:
          print('Sparse SVD ...')
  
      X0, S0, Y0 = svds_descending(M_Et, r)
      del M_Et
  
      # Initial Guess
      X0 = X0 * np.sqrt(n)
      Y0 = Y0 * np.sqrt(m)
      S0 = S0 / eps
  
      if print_out:
          print('Iteration\tFit Error')
  
      # Gradient Descent
      X = copy.deepcopy(X0)
      Y = copy.deepcopy(Y0)
      S = getoptS(X, Y, M_E, E)
  
      errs = [None] * (niter + 1)
      errs[0] = np.linalg.norm(
          (M_E - np.dot(np.dot(X, S), Y.T)) * E, ord='fro') / np.sqrt(np.count_nonzero(E))
  
      if print_out:
          print('0\t\t\t%e' % errs[0])
  
      for i in range(niter):
          # Compute the Gradient
          W, Z = gradF_t(X, Y, S, M_E, E, m0, rho)
  
          # Line search for the optimum jump length
          t = getoptT(X, W, Y, Z, S, M_E, E, m0, rho)
          X = X + t * W
          Y = Y + t * Z
          S = getoptS(X, Y, M_E, E)
  
          # Compute the distortion
          errs[i + 1] = np.linalg.norm((M_E - np.dot(np.dot(X, S), Y.T))
                                       * E, ord='fro') / np.sqrt(np.count_nonzero(E))
          if print_out:
              print('%d\t\t\t%e' % (i + 1, errs[i + 1]))
  
          if abs(errs[i + 1] - errs[i]) < tol:
              break
      S = S / rescal_param
  
      return X, S, Y, errs