If you can live with some slight "inaccuracies" on the borders of the polygons, you can achieve a very fast processing by rasterizing the shapefile, and then extracting the value of the points on the raster using raster::extract. Something like this would work, and outputs a SpatialPointsDataFrame.

library(rgdal)
library(gdalUtils)
library(sf)
library(raster)
library(sf)
library(lazyeval)
library(tidyr)
download.file("http://gis.ices.dk/shapefiles/ICES_ecoregions.zip",
              destfile = "ICES_ecoregions.zip")
unzip("ICES_ecoregions.zip")

# read eco region shapefiles
ices_eco <- sf::read_sf("ICES_ecoregions_20150113_no_land.shp")
## Make a large data.frame (361,722 rows) with positions in the North Sea:
lon <- seq(-18.025, 32.025, by=0.05)
lat <- seq(48.025, 66.025, by=0.05)
c <- expand.grid(lon=lon, lat=lat)

# Get the Ecoregion for each position
pings  <- SpatialPointsDataFrame(c, 
                                 data = data.frame(id = 1:dim(c)[1]), 
                                 proj4string = CRS(st_crs(ices_eco)$proj4string))

extr_data = list()
temprast  <- tempfile(fileext = ".tif")
tempshape <- tempfile(fileext = ".shp")
counter = 1
for (poly_n in seq_along(along = ices_eco$Ecoregion))  {

  message("Working on polygon: ", poly_n)
  # extract one polygon andrasterize it with 0.1 deg resolution
  subshape    <- ices_eco[poly_n,]
  st_write(subshape, tempshape, update = TRUE, quiet = TRUE)

  # crop the points dataframe on the extent of the rasterized polygon
  # to save time
  subpoints   <- raster::crop(pings, extent(as.numeric(st_bbox(subshape))[c(1,3,2,4)]))

  # if there are points "left", extract the value of the rasterized polygon, for each point
  if (!is.null(subpoints)) {
    rast <- gdal_rasterize(tempshape, temprast, burn = poly_n,
                           tr = c(0.01,0.01), a_nodata = -999,
                           output_Raster = TRUE)
    extr_points <- raster::extract(rast, subpoints, sp = TRUE)

    extr_points <- extr_points[which(extr_points@data[ , 2] == poly_n), ] %>%
      as("sf")

    # If any of the points fell in the polygon, extract the value and return them in
    # the list. otherwise, return null

    if ((dim(extr_points)[1]) != 0 ) {
      extr_data[[counter]] <- extr_points
      counter <- counter + 1
    }
  }
  unlink(temprast)
  unlink(tempshape)
}

# Do some juggling to remove duplicates and go back to a SpatialPointsDataFrame
out <- data.table::rbindlist(extr_data) %>%
  tibble::as_tibble() %>% 
  # Following two lines needed to remove few "duplicated" points. I think you get those
  # for points falling on the boundary between two cells. Increasing resolution in rasterization 
  # could avoid this
  dplyr::group_by(id) %>% 
  dplyr::summarize_(ecoregion = interp(~first(var), var = as.name(names(extr_data[[1]])[2])),
                    geometry  = interp(~vargeo[1],  vargeo = as.name("geometry"))) %>% 
  # transform bakck to a spatialpointsdataframe
  sf::st_as_sf() %>%
  as("Spatial") 

# Add the extracted info to the original SpatialPointsDataFrame
pings@data <- dplyr::left_join(pings@data, out@data) 
head(pings@data)
summary(pings)

>     head(pings@data)
  id ecoregion
1  1        17
2  2        17
3  3        17
4  4        17
5  5        17
6  6        17
>     summary(pings)
Object of class SpatialPointsDataFrame
Coordinates:
        min    max
lon -18.025 32.025
lat  48.025 66.025
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Number of points: 361722
Data attributes:
       id           ecoregion     
 Min.   :     1   Min.   : 9.00   
 1st Qu.: 90431   1st Qu.: 9.00   
 Median :180862   Median :11.00   
 Mean   :180862   Mean   :11.86   
 3rd Qu.:271292   3rd Qu.:15.00   
 Max.   :361722   Max.   :17.00   
                  NA's   :180720  

On my PC, this completes in about 1 minute. Note that in the code above, Im also processing one "polygon" at a time to further increase speed and avoid creating a huge temporary raster file. Additionally, Im using the sf functions for reading and writing shapefiles, which are much faster than the rgdal ones.

Also, I'm rasterizing the shapefile to a 0.01x0.01 degrees resolution'. You can reduce the possibility of "errors" by increasing the resolution of the rasterization, by changin the tr parameter in the instruction rast <- gdal_rasterize(tempshape, temprast, burn = poly_n, tr = c(0.01,0.01), a_nodata = -999, output_Raster = TRUE)

HTH !