Q:

Many points in many polygons, ID in R

I have a shapefile with about ~5k polygons.  I have a database of about 3.8M points.  
I want to know, for every point, which polygon it's in. 
For example, Point A may be in Polygon A, Polygon B, and Polygon C. 
The problem I have is two fold: 
1) How do I tie a specific point match to a specific polygon?  Right now, I can only get the sum total of the number of points within  a polygon using the following code: 
colSums(gContains(all_polygons, spdf, byid = TRUE)) -> res

2) How do I make this computationally easier?  I am less familiar with GIS analysis in R, but in traditional SQL language, we may do a cross join of all polygons with all points, and then throw a flag 1,0 if the point is in the polygon.  This is incredibly computationally intense, as it generates a table that is 5k* 3.8M. Is there another way around this?  

A:

If you use the sf package, then st_intersection can get you all the intersections, with none of the non-intersections...
Sample data - some points and two polygons which overlap.  The polygons have a "Region" attribute which is A or B:

Then if I do:
inters = st_intersection(points,polys)

I get:
       Name Region                       geom
4         D      A POINT (-15.73582 13.59678)
5         R      A POINT (-15.68296 13.56941)
6         H      A POINT (-15.63011 13.57446)
7        KK      A POINT (-15.58938 13.56983)
8       IIJ      A  POINT (-15.53783 13.5972)
10        Y      A POINT (-15.46808 13.56646)
10.1      Y      B POINT (-15.46808 13.56646)
11   ksjdkj      B POINT (-15.46981 13.64015)
12   kdjihj      B POINT (-15.39616 13.61362)
14   ksjicc      B POINT (-15.27443 13.62752)

And the columns are the attributes from each overlapping point-polygon combination - Name and geometry comes from  points, and Region comes from polys. 
The points that aren't in polygons do not appear in this result. 
If speed is an issue, your problem is trivially parallel over polygons, so you can divide your polygons into N sets and run on N processors for a potential N-fold speed increase. 
For example, suppose you have 1200 polygons in polys.  Suppose computing st_intersection(points, polys) takes 10 seconds.  Computing the intersection with half the polygons with st_intersection(points, polys[1:600,]) should take five seconds.  If you can do that at the same time as st_intersection(points, polys[601:1200,]) then you can do all 1200 polygon tests in five seconds - and then a short time to join the two results. 
If you have an 8-core processor then you could divide your polygons into 8 groups, and get an 8-fold speedup, using R's various parallel packages.  
If you have access to some cloud resources, like Amazon EC2 or Google Cloud or Azure, you can "spin up" 600 virtual machines, have each one test 2 polygons, and get a potential massive speedup.