Edit Now this answers the first question for the operator norm and the normalized Hilbert-Schmidt norm.
The answer depends on the norm you are considering. The answer is no for the operator norm, but is yes for the normalized Hilbert-Schmidt norm (at least if you replace $O(\varepsilon)$ by $o(1)$, see the answers to this question).
Here are some details on the counterexample for the operator norm.
- By a theorem of Lin (see here), for a pair of self-adjoint matrices of norm less than $1$, they approximately commute if and only if they can be approximated by commuting matrices.
- Voiculescu proved that the preceding does not hold for triples of self-adjoint matrices of norm less than $1$ (the link I gave here is not correct ; see for example the references in the paper by Exel and Loring given in the comments).
1+2 imply that there is a sequence of triples $A_1^n,A_2^n,A_3^n$ of matrices of norm less than $1$ which are pairwise close to (self-adjoint) commuting matrices, but whose distance to the triples of commuting matrices is bounded below.
- By continuity of the functional calculus and the fact that $t \in [-2,2] \mapsto e^{it}$ is a homeomorphism on its image, this implies that the unitary matrices $(e^{i A_1^n}, e^{i A_1^n},e^{i A_1^n})$ are pairwise close to pairs of commuting unitaries, but are at positive distance from triples of commuting unitaries. This is what you were looking for.