Spatial 1 to 1 join by proximity in sf and R?

I'm trying to join a Simple Feature of POINT geometry type to another Simple Feature of POLYGON geometry type using st_join and nngeo::st_nn as the person asking this other question also did. The difference is that I want a 1:1 match from both sides (i.e. one point per polygon).

Above is a picture of how my data looks like when mapped. My code looks like:

matched<-st_join(buildings, addresses, join=nngeo::st_nn, maxdist=50,k=1)

Instead of a 1 to 1 match, I keep getting more than one point merged to the same polygon or more than one polygon merged to the same point depending on how order the arguments to st_join, with repeated points or polygons accordingly.

What I want this to do is to match each polygon to the closest point which also has that polygon as the closest to itself (i.e., to match polygons to points when they are mutually the closest to each other), and to randomise matching when more than one point or polygon meet those criteria.

Is there any way of doing that with sf or sp in r? I have thought about extracting the distance matrices and doing a two-sided matching between the two sets of features (like a Gale- Shapley algorithm), but I want to give it a shot with what is already available before putting many hours into coding that.