1. Load the following Libraries

# Libraries
library(rgeos)

## Loading required package: sp

## rgeos version: 0.4-3, (SVN revision 595)
##  GEOS runtime version: 3.6.2-CAPI-1.10.2 
##  Linking to sp version: 1.3-1 
##  Polygon checking: TRUE

library(maptools)

## Checking rgeos availability: TRUE

library(raster)
library(rgdal)

## rgdal: version: 1.4-4, (SVN revision 833)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: /usr/share/gdal/2.2
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.3-1

library(parallel) 

The library parallel can be used to either run your script in parallel through parLapply or any similar of your preference. In this case, I use it to detect the number of cores in the system.

#2. Check the number of cores that you have

# Calculate the number of cores
no_cores <- detectCores() - 1

#3. Input & Stack your data

#List files (TIFFs or any other format)
Raster_Files <- list.files("/home/jorge/Windows Share/Final VHI", pattern = "*.tif$")
#Check the number of files in the raster list 
Number_Files <- length(Raster_Files)
#read-in the polygon shapefile
poly <- readShapePoly("/home/jorge/Windows Share/Final VHI/IDN_adm2.shp")

## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read

poly2 <- readShapePoly("/home/jorge/Windows Share/Final VHI/IDN_adm1.shp") 

## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read

#Check that the spatial points file has been sucessfully uploaded
plot(poly)

#create a raster stack
st_Files <- stack(paste0("/home/jorge/Windows Share/Final VHI/", Raster_Files))

#4 Process your Data

################################################################################
beginCluster(n=17)
#extract raster cell count (sum) or any statistic that you want (MEDIAN/MEAN/STDV/.../Etc) within each polygon area (poly)
for (i in 1:length(Raster_Files)){
  Extracted_Data <- extract(st_Files, poly2, fun=median, na.rm=TRUE, df=TRUE)
}

## Using cluster with 17 nodes
## Using cluster with 17 nodes
## Using cluster with 17 nodes
## Using cluster with 17 nodes

endCluster()
#write to a data frame
Stat_File <- data.frame(Extracted_Data)
#Check the Result
View(Stat_File)
################################################################################

#5 Export your resulting Data

########
#write to a CSV file
write.csv(Stat_File, file = "/home/jorge/Windows Share/Final VHI/MEDIAN.csv")