Performance Optimization – Neighboring Surfaces in a Discrete Network

I have a mesh region R and I want to efficiently extract all pairs of adjacent surfaces. The way I calculated this is shown below:

R = Discretize region[Sphere[]];
Faces = MeshCells[R, 2, "Multicells" -> True][[1, 1]];
adj = Select[Subsets[faces, {2}],Length[Intersection[#[[1]], #[[2]]]]== 2 &];

Of course, this is not very efficient because I explicitly construct all the surface pairs and then filter them by assuming that both contain two equal vertices. Considerations on how to calculate the same efficiently?