I did a benchmark to compare my own (improved) method with the ones proposed by [kglr](https://mathematica.stackexchange.com/users/125) and [Carl Woll](https://mathematica.stackexchange.com/users/45431). My initial thought was to treat this problem in a procedural, old-fashioned way. Here is the code:
```
ClearAll["Global`*"];
merge = Which[AnyTrue[#, Length@# < 2 &], #,
    #[[1, -1]] == #[[2, 1]], {#[[1]]~Join~Rest@#[[2]], {}},
    #[[1, 1]] == #[[2, -1]], {{}, #[[2]]~Join~Rest@#[[1]]},
    True, #] &;

pathsFind[li_List] := Module[{a = li},
  Do[
   Do[a[[{i, j}]] = merge@a[[{i, j}]], {j, 2, Length@a}, {i, 1, j - 1}];
   If[FreeQ[a, {}], Break[], a = a~DeleteCases~{}],
   Length@a];
  a]
```
Nothing fancy. It just loops through the elements multiple times until all paths are found. This code compares `allPaths` by kglr, `findPaths` by Carl Woll and `pathsFind`:
```
Module[{li, p, rand},
 li = DeleteCases[{a_, a_}]@DeleteDuplicatesBy[Sort]@RandomInteger[1000, {500, 2}];
 First /@ {
   RepeatedTiming[p = findPaths@li;],
   {Length@p},
   RepeatedTiming[p = allPaths@li;],
   {Length@p},
   RepeatedTiming[p = pathsFind@li;],
   {Length@p}
   }]
```
The result is `{0.417195, 372, 0.0131587, 372, 1.68831, 346}`. Obviously, `allPaths` is faster by some orders of magnitude. But it seems both faster methods are not so accurate and miss some paths along the way, since the number of groups in the last method is less than the other two.

By the way, thanks to all the nice ideas and efforts proposed here. I think I should better stick to my method for now.