Correlations of two flow harmonics vn and vm via three- and four-particle cumulants are measured in 13 TeV pp, 5.02 TeV p+Pb, and 2.76 TeV peripheral Pb+Pb collisions with the ATLAS detector at the LHC. The goal is to understand the multi-particle nature of the long-range collective phenomenon in these collision systems. The large non-flow background from dijet production present in the standard cumulant method is suppressed using a method of subevent cumulants involving two, three and four subevents separated in pseudorapidity. The results show a negative correlation between v2 and v3 and a positive correlation between v2 and v4 for all collision systems and over the full multiplicity range. However, the magnitudes of the correlations are found to depend on the event multiplicity, the choice of transverse momentum range and collision system. The relative correlation strength, obtained by normalisation of the cumulants with the 〈vn2〉 from a two-particle correlation analysis, is similar in the three collision systems and depends weakly on the event multiplicity and transverse momentum. These results based on the subevent methods provide strong evidence of a similar long-range multi-particle collectivity in pp, p+Pb and peripheral Pb+Pb collisions.