code: def _M2_sparse_sym(Xvar, mask_X, Yvar, mask_Y, weights=None, column_selection=None):
      
      assert len(mask_X) == len(mask_Y), 'X and Y need to have equal sizes for symmetrization'
  
      if column_selection is None:
          mask_Xk = mask_X
          mask_Yk = mask_Y
          Xvark = Xvar
          Yvark = Yvar
      else:
          mask_Xk = mask_X[column_selection]
          mask_Yk = mask_Y[column_selection]
          Xvark = Xvar[:, _filter_variable_indices(mask_X, column_selection)]
          Yvark = Yvar[:, _filter_variable_indices(mask_Y, column_selection)]
  
      Cxxyy = np.zeros((len(mask_X), len(mask_Yk)))
      Cxxyy[np.ix_(mask_X, mask_Xk)] = _M2_dense(Xvar, Xvark, weights=weights)
      Cxxyy[np.ix_(mask_Y, mask_Yk)] += _M2_dense(Yvar, Yvark, weights=weights)
  
      Cxyyx = np.zeros((len(mask_X), len(mask_Yk)))
      Cxy = _M2_dense(Xvar, Yvark, weights=weights)
      Cyx = _M2_dense(Yvar, Xvark, weights=weights)
      Cxyyx[np.ix_(mask_X, mask_Yk)] = Cxy
      Cxyyx[np.ix_(mask_Y, mask_Xk)] += Cyx
  
      return Cxxyy, Cxyyx
docs: 2nd self-symmetric moment matrix exploiting zero input columns
  
      Computes X'X + Y'Y and X'Y + Y'X