It may be too difficult to characterize all matrices that are products of pairwise averaging matrices, so it is probably a good idea to try to solve an easier problem (at least computationally) in a quotient monoid in order to obtain and compute many counterexamples.
Let $[n]=\{1,\dots,n\}$. If $f\in S_{n}$, then let $\rho_{f}$ denote the permutation matrix $(\delta_{f(i),j})_{i,j}$ that corresponds to the permutation $f$. Suppose that $A\subseteq[n]\times[n]$. Then we say that $A$ is an $n\times n$-multipermutation if $A$ is the union of permutations. Let $M_{n}$ be the set of all $n\times n$-multipermutations.
Let $V_{n}=P([n]^{2})$, and define a composition operation in $V_{n}$ by $$A\circ B=\{(x,z)|\exists y,(x,y)\in B,(y,z)\in A\}.$$ Let $DS_{n}$ be the set of all doubly stochastic matrices. Define a mapping $\phi:DS_{n}\rightarrow M_{n}$ by letting $\phi((a_{i,j})_{i,j})=\{(i,j)\mid a_{i,j}>0\}$. Then $\phi$ is a surjective monoid anti-homomorphism (i.e. $\phi(BA)=\phi(A)\circ\phi(B)$ for doubly stochastic $A,B$, and $\phi$ preserves the identity).
Let $B_{n}$ be the submonoid of $M_{n}$ generated by elements of the form
$f\cup g$ for some $f,g\in S_{n}$. Then $R\in B_{n}$ if and only if $R=\phi(A)$ for some matrix $A$ that can be factored as a product of pairwise averaging matrices.
Therefore, if we can show that $\phi(A)\not\in B_{n}$, then we know that
$A$ cannot be factored as a product of pairwise averaging matrices.
Now, each element in $B_{n}$ can be put into the form $(1\cup f_{1})\circ\dots\circ(1\cup f_{r})\circ g$ or in the form
$u\circ (1\cup v_{1})\circ\dots\circ(1\cup v_{r})$. If we define $B_{n}^{\sharp}$ to be the submonoid of $B_{n}$ generated by the elements of the form $1\cup f$, then $B_{n}$ can be thought of as a sort of semidirect product of the permutation group $S_{n}$ and the monoid $B_{n}^{\sharp}$. In fact, if we define a monoid operation $*$ on $V_{n}\times S_{n}$ by letting $(R,f)*(S,g)=(RfSf^{-1},fg)$, and we define a monoid homomorphism $\Gamma:V_{n}\times S_{n}\rightarrow V_{n}$ by letting
$\Gamma(R,f)=Rf$, then $B_{n}=\Gamma[B_{n}^{\sharp}\times S_{n}]$.
Let $D_{n}=\{(R,(R\circ(1\cup f)))\mid R\in V,f\in S_{n}\}$. Then $G_{n}=(V_{n},D_{n})$ is an acyclic graph, so we can partially order $V_{n}$ by letting $f\leq g$ if and only if there is a directed path from $f$ to $g$. Furthermore, the partial ordering $\leq$ is compatible with $\subseteq$ since if $f\leq g$, then $f\subseteq g$ (this is useful since $\leq$ is difficult to compute, but $\subseteq$ is easy to compute). Observe that $R$ belongs to $B_{n}$ if and only if there exists some path from some permutation $f$ to $R$ in $(V_{n},D_{n})$. Define $\downarrow R=\{S\in V_{n}\mid S\subseteq R\}$. Since every path from $f$ to $R$ only goes through nodes in the set $\downarrow R$, we can restrict our path finding algorithm to induced subgraph whose nodes are precisely the elements in $\downarrow R$. Therefore, $R\in B_{n}$ if and only if there is some permutation $f$ and a path from $f$ to $R$ in the induced subgraph $G_{n}[\downarrow R]$.
Now, observe that $$N_{G_{n}[\downarrow S]}(R)=\{R\circ(1\cup f)\mid f\in S_{n},R\circ(1\cup f)\subseteq S\}=\{R\circ(1\cup f)\mid f\in S_{n},R\circ f\subseteq S\}.$$ Therefore, to calculate, $N_{G_{n}[\downarrow S]}(R)$, it suffices to find all permutations $f\in S_{n}$ such that $R\circ f\subseteq S$. However, if we define
$T_{y}=\{z\in[n]\mid\forall x,(x,y)\in R\rightarrow(x,z)\in S\}$, then
$R\circ f\subseteq S$ if and only if $f(y)\in T_{y}$ for all $y\in[n]$. Therefore, the problem of finding the elements in $N_{G_{n}[\downarrow S]}(R)$ reduces to the problem of finding all perfect matchings in bipartite graphs. However, according to this paper, the problem of finding all perfect matchings in a bipartite graph takes $O(c(n+m)+n^{5/2})$ computational effort and $O(nm)$ memory where $n$ is the number of vertices, $m$ is the number of edges, and $c$ is the number of perfect matchings in the bipartite graph. However, I am not convinced that finding all bijections $f$ with $R\circ f\subseteq S$ is the best strategy for
computing $N_{G_{n}[\downarrow S]}(R)$ since there may be many more permutations $f$
with $R\circ f\subseteq S$ than elements in $N_{G_{n}[\downarrow S]}(R)$.
An enhanced graph
Define a graph $G_{n,S}=(V_{n,S},E_{n,S})$ where
$V_{n,S}$ consists of tuples $(P,A,B,Q)$ such that $|A|=|B|$ and
$Q\subseteq R\cap(B^{c}\times[n])$.
$E_{n,S}$ consists of the edges of the following forms
i. $((P,[n],[n],\emptyset),(P,\emptyset,\emptyset,P))\in E_{n,S}$
ii. Suppose that $A\subseteq[n],B\subseteq[n]$, and $a\in A^{c},b\in B^{c}$. Then $$((P,A,B,Q),(P\cup(Q\circ\{(a,b)\}),A\cup\{a\},B\cup\{b\},Q\setminus(\{b\}\times[n])))\in E_{n,S}$$ such that
$P\cup(Q\circ\{(a,b)\})\subseteq S$.
Suppose that $f$ is a permutation. Then there is a directed path from
$(f,\emptyset,\emptyset,f)$ to $(P,A,B,Q)$ if and only if there are
$g_{1},\dots,g_{r}$ such that if $R=f\circ(1\cup g_{1})\circ\dots\circ(1\cup g_{r})$, then there is a some bijection
$h:A\rightarrow B$ such that $P=R\cup(R\circ h)$ and $Q=R\cap(B^{c}\times[n])$. Therefore, $S\in B_{n}$ if and only if there is a path from $(f,\emptyset,\emptyset,f)$ to $(S,\emptyset,\emptyset,S)$ for some permutation $f$ with $f\subseteq S$.
Therefore, one can determine whether $S\in B_{n}$ by using a path-finding algorithm. However, our path finding algorithm will probably not be very efficient at this point, so let us reduce the number of edges from $(f,\emptyset,\emptyset,f)$ to $(S,\emptyset,\emptyset,S)$ in order to improve efficiency. We say that a subset $E\subseteq E_{n,S}$ is traversable if $E$ contains each edge of the form
$((P,[n],[n],\emptyset),(P,\emptyset,\emptyset,P))$ and where if
$(P,A,B,Q)\in V_{n,S}$ and $A\neq[n]$ and $h:A^{c}\rightarrow B^{c}$ is bijective, then there are $a\in A^{c},b\in B^{c}$ with $b=h(a)$ and where
$$((P,A,B,Q),(P\cup(Q\circ\{(a,b)\}),A\cup\{a\},B\cup\{b\},Q\setminus(\{b\}\times[n])))\in E.$$
Observe that if $E$ is traverable, then $S\in B_{n}$ if and only if there is some path from $(f,\emptyset,\emptyset,f)$ to $(S,\emptyset,\emptyset,S)$ in the directed graph $(V_{n,S},E)$, so a path finding algorithm may be used to determine whether $S\in B_{n}$ or not.
Code
Here is the code in GAP for determining whether $R\in B_{n}$. I have made no attempt to optimize this code.
allperms:=function(x)
local list,n,table,i,newtable,j;
n:=Length(x); table:=[[]];
for i in [1..n] do
newtable:=[];
for list in table do
for j in [1..n] do
if j in list then continue; fi;
if x[i][j]=1 then Add(newtable,Concatenation(list,[j])); fi;
od; od;
table:=newtable; od;
return table; end;
implicator:=function(x,y)
local n,z,i,j,k;
n:=Length(x); z:=x*0;
for i in [1..n] do for j in [1..n] do
z[i][j]:=1;
for k in [1..n] do if y[k][j]=0 and x[k][i]=1 then z[i][j]:=0; break; fi;
od; od; od;
return z; end;
neighbors:=function(x,y)
local list,n,xx,z,perms,output,i,j;
n:=Length(x); z:=implicator(x,y); perms:=allperms(z); output:=[];
for list in perms do
xx:=StructuralCopy(x);
for i in [1..n] do for j in [1..n] do
if x[i][j]=1 then xx[i][list[j]]:=1; fi;
od; od;
Add(output,xx);
od;
return output; end;
size:=function(x) return Sum(List(x,Sum)); end;
test:=function(y)
local a,perms,xx,n,s,i,j,tar,table,har,bb,qq;
n:=Length(y); s:=size(y); perms:=allperms(y); tar:=NewDictionary(true,true); table:=[];
for a in perms do xx:=y*0;
for i in [1..n] do xx[i][a[i]]:=1; od;
Add(table,xx); AddDictionary(tar,xx,true);
od;
while Length(table)>0 do
qq:=Remove(table); har:=neighbors(qq,y);
for bb in har do if size(bb)=s then return true; fi; if not (LookupDictionary(tar,bb)=true) then AddDictionary(tar,bb,true); Add(table,bb); fi; od;
od;
return false;
end;
More code
Here is the code that finds a path in the graph $(V_{n,S},E)$ for some traverable set $E$.
newnode:=function(quadruple) local PP; PP:=quadruple[1]; return [PP,[],[],PP]; end;
testquad:=function(quadruple,pair,y) local a,b,P,Q,A,B,n,i; a:=pair[1]; b:=pair[2]; P:=quadruple[1]; A:=quadruple[2]; B:=quadruple[3]; Q:=quadruple[4]; n:=Length(P);
for i in [1..n] do if Q[b][i]=1 and y[a][i]=0 then return false; fi; od;
return true; end;
newpairzero:=function(n,zero,one) local m,l,lar,mar; m:=Length(zero); lar:=Difference([1..n],zero); mar:=Difference([1..n],one); l:=Minimum(lar); m:=Minimum(mar); return List(mar,v->[l,v]); end;
newpairone:=function(n,zero,one) local m,l,lar,mar; m:=Length(zero); lar:=Difference([1..n],zero); mar:=Difference([1..n],one); l:=Minimum(lar); m:=Minimum(mar); return List(lar,v->[v,m]); end;
extendpair:=function(quadruple,pair) local a,b,P,Q,A,B,n,i,AA,BB,PP,QQ; a:=pair[1]; b:=pair[2]; P:=quadruple[1]; A:=quadruple[2]; B:=quadruple[3]; Q:=quadruple[4]; n:=Length(P); QQ:=StructuralCopy(Q); AA:=Concatenation(A,[a]); BB:=Concatenation(B,[b]); Sort(AA); Sort(BB); PP:=StructuralCopy(P);
for i in [1..n] do PP[a][i]:=Maximum(QQ[b][i],PP[a][i]); QQ[b][i]:=0; od;
return [PP,AA,BB,QQ]; end;
test:=function(y) local P,Q,pair,list,n,tar,table,xx,i,j,A,B,quadruple; n:=Length(y); tar:=NewDictionary(true,true); table:=[]; xx:=[];
for i in [1..n] do xx[i]:=[]; for j in [1..n] do xx[i][j]:=0; od; xx[i][i]:=1; od;
Add(table,[xx,[],[],xx]);
while Length(table)>0 do quadruple:=Remove(table);
if not LookupDictionary(tar,quadruple)=fail then continue; fi;
AddDictionary(tar,quadruple,true); P:=quadruple[1]; A:=quadruple[2];
B:=quadruple[3]; Q:=quadruple[4];
if Length(A)=n then if P=y then return true; fi;
Add(table,newnode(quadruple)); else list:=newpairzero(n,A,B);
for pair in list do if testquad(quadruple,pair,y) then
Add(table,extendpair(quadruple,pair)); fi; od;
fi;
od;
return false; end;
Empirical observations
I ran the above algorithm to determine whether elements belonged to $B_{n}$, and most elements in $M_{n}$ that I have tested do not belong to $B_{n}$. In fact, in most cases, elements of the form $f\cup g\cup h$ where $f,g,h\in S_{n}$ do not belong to $B_{n}$, so most stochastic matrices of the form $\frac{1}{3}(\rho_{f}+\rho_{g}+\rho_{h})$ cannot be written as the product pairwise averaging matrices.
I will probably add details and improve this answer some time later in order to improve the algorithm for determining whether an element $R\in V_{n}$ belongs to $B_{n}$ or not.
Singular matrices
Suppose that $f$ is a permutation that can be written as the composition of $r$ disjoint cycles of lengths $k_{1},\dots,k_{r}$.
Then $$\text{Det}(t\cdot 1_{n}+(1-t)\rho_{f})=\prod_{j=1}^{r}[t^{k_{j}}+(t-1)^{k_{j}}].$$ In particular, $t\cdot 1_{n}+(1-t)\rho_{f}$ is non-singular if and only if $t=1/2$, and if $t=1/2$, then
$\text{Rank}(t\cdot 1_{n}+(1-t)\rho_{f})$ is the number of odd cycles in the permutation $f$.
If $R\in B_{n}$, then let $N(R)$ be the supremum of all sums $m_{1}+\dots+m_{r}$ such that $R=f(1+g_{1})\dots(1+g_{r})$ and when $g_{i}$ is written as a disjoint composition of cycles, $m_{i}$ of those cycles are odd. If $R\not\in B_{n}$, then set $N(R)=\infty$. We can adjust our algorithm for determining whether $R\in B_{n}$ or not so that it calculates the value $N(R)$ as well.
Then observe that if $A$ is a doubly stochastic matrix that can be written as a product of pairwise averaging matrices, then $\text{Null}(A)\leq N(\phi(A))$. Therefore, if $\text{Null}(A)>N(\phi(A))$, then we know that $A$ cannot be written as a product of pairwise averaging matrices.
Doubly stochastic matrices near a permutation matrix
The problem of determining whether $R\in B_{n}^{\sharp}$ or not is easier to solve than the problem of determining whether $R\in B_{n}$ since for $B_{n}^{\sharp}$, the path finding algorithm starts with a single vertex, namely $(1,\emptyset,\emptyset,1)$ instead of many vertices. Fortunately, the problem of determining whether a set belongs to $B_{n}^{\sharp}$ helps us prove that certain doubly stochastic matrices cannot be written as a product of pairwise averaging matrices.
If $A$ is a doubly stochastic matrix that is too close to a permutation matrix $\rho_{f}$, and $A$ can be written as a product of pairwise averaging matrices, then we know that $A$ can be put in the form
$A=\rho_{f}(t_{1}I_{n}+(1-t_{1})\rho_{g_{1}})\dots(t_{r}I_{n}+(1-t_{r})\rho_{g_{r}})$ for constants $t_{1},\dots,t_{r}\in(0,1)$ where the constants $t_{1},\dots,t_{r}$ are near $1$. Let us make this intuition precise.
Theorem: Let $\|\cdot\|$ be a norm on the vector space of all $n\times n$-matrices. Let $$M=\max_{f\in S_{n}}\|\rho_{f}\|,$$ and let
$$N=\min_{f\in S_{n}\setminus\{1\}}\|I_{n}-\rho_{f}\|.$$ Then whenever
$A$ is the product of pairwise averaging matrices, then
- $$A=(t_{1}I_{n}+(1-t_{1})\rho_{f_{1}})\cdots(t_{r}I_{n}+(1-t_{r})\rho_{f_{r}})$$
where $\min(t_{1},\dots,t_{r})\geq 1/2$
or
- $$\|A-I_{n}\|\geq N-2M(1-|\det(A)|).$$
Proof: Suppose that $f$ is a permutation that can be decomposed into disjoint cycles of lengths $k_{1},\dots,k_{r}$. Observe that
$$|\det(tI_{n}+(1-t)\rho_{f})|=\prod_{j=1}^{r}|t^{k_{j}}+(t-1)^{k_{j}}|
\leq|t^{k_{1}}+(t-1)^{k_{1}}|$$
$$\leq t^{k_{1}}+(1-t)^{k_{1}}\leq t^{2}+(1-t)^{2}\leq\frac{1}{2}+|t-\frac{1}{2}|.$$
Now, suppose that $A$ is the product of pairwise averaging matrices. Then without loss of generality, there are $t_{1},\dots,t_{r}$ such that $\min(t_{1},\dots,t_{r})\geq\frac{1}{2}$ along with permutations $f,g_{1},\dots,g_{r}$ where
$$A=\rho_{f}(t_{1}I_{n}+(1-t_{1})\rho_{g_{1}})\dots(t_{r}I_{n}+(1-t_{r})\rho_{g_{r}})$$.
In this case, we have $\det(A)\leq t_{1}\dots t_{r}$.
Now, $$\|A-\rho_{f}\|\leq\|A-\rho_{f}t_{1}\dots t_{r}I_{n}\|+\|\rho_{f}t_{1}\dots t_{r}I_{n}-\rho_{f}\|$$
$$\leq(1-t_{1}\dots t_{r})M+(1-t_{1}\dots t_{r})M\leq 2M(1-t_{1}\dots t_{r})\leq 2M(1-\det(A)).$$
Therefore, $$\|A-I_{n}\|\geq\|I-\rho_{f}\|-\|A-\rho_{f}\|\geq N-2M(1-\det(A)).$$ Q.E.D.
Observe that $M=\max\{\|A\|:\text{$A$ is doubly stochastic}\}$. Furthermore, $\|N\|\leq\|I_{n}\|+M\leq 2M$. If the norm $\|\cdot\|$ is submultiplicative ($\|AB\|\leq\|A\|\cdot\|B\|$ for all $A,B$), then $M=\|I\|=I_{n},\|N\|\leq 2$.
Theorem: Let $\|\cdot\|$ be a norm on the set of all $n\times
n$-matrices. Let $$M=\max_{f\in S_{n}}\|\rho_{f}\|,N=\min_{f\in
S_{n}\setminus\{1\}}\|I_{n}-\rho_{f}\|.$$
Suppose that $A$ is a doubly stochastic matrix. If $\phi(A)\not\in
B_{n}^{\sharp}$, and $$|A-I_{n}|<N-2M(1-\det(A)),$$ then $A$ cannot be
written as the product of pairwise averaging matrices.
Argument from dimensionality
Suppose that $R\in B_{n}$. Then let $T_{R}$ be the vector subspace of $M_{n}(\mathbb{R})$ consisting of all matrices $(a_{i,j})_{i,j}$ such that
$\{(i,j)\mid a_{i,j}\neq 0\}\subseteq R$ and where $\sum_{k}a_{i,k}=0$ for each $i\in[n]$ and where $\sum_{k}a_{k,j}=0$ for each $i\in[n]$. Observe that $\dim(T_{R})\geq|R|-2n$.
Then whenever $\phi(A)=R$, there is an open subset $U$ of $T_{R}$ such that
$A+B$ is doubly stochastic for each $B\in U$.
Proposition: Suppose that $R\subseteq[n]^{2}$, $R$ is the union of permutations, but $R$ cannot be put into the form $R=f(1\cup f_{1})\dots(1\cup f_{r})$ where $f_{1},\dots,f_{r}$ are non-identity permutations and $f$ is a permutation, and where $r\geq\dim(T_{R})$. Then there is a doubly stochastic matrix $A$ with $\phi(A)=R$ but where
$A$ cannot be written as a product of doubly stochastic matrices.
Proof: Simply observe that the dimension of the set of all doubly stochastic matrices $A$ with $\phi(A)=R$ is $\dim(T_{R})$ while the dimension of the set of all products of pairwise averaging matrices $A$ with $\phi(A)=R$ is at most $\dim(T_{R})-1$. Q.E.D.