Prove that with some conditions,$\sum\limits_{i=1}^mA_iBA_i^T=B$ iff $A_iB=BA_i(1\leq i\leq m)$. $A_i\in\mathbb{R}^{n\times n}(1\leq i\leq m)$ ,$\sum\limits_{i=1}^mA_iA_i^T=I$,$B\in\mathbb{R}^{n\times n}$ and $\text{tr}(B)=1$,$B$ is positive.
Prove that $\sum\limits_{i=1}^mA_iBA_i^T=B$ iff $A_iB=BA_i(1\leq i\leq m)$.
One direction is simple:if $A_iB=BA_i(1\leq i\leq m)$,then  $\sum\limits_{i=1}^mA_iBA_i^T=\sum\limits_{i=1}^mBA_iA_i^T=B$,but what about the other direction?
Update:Thank @mechanodroid for answer and I realize that the condition:$\sum\limits_{i=1}^mA_i^TA_i\preceq I$ should be added after reviewing the background of this problem.The form $\sum A_iBA_i^T$is called the Kraus representation and the property $\sum\limits_{i=1}^mA_iBA_i^T=B$ is called "unital".And for Kraus representation there is a default rule $\sum\limits_{i=1}^mA_i^TA_i\preceq I$ (which is called trace-preserving) as mentioned by @mechanodroid, but I omitted that,my fault.
 A: Consider the square of the Frobenius norm of $A_iB-BA_i$, which is the trace of $$\begin{aligned}&(A_iB-BA_i)(A_iB-BA_i)^t \\=&(A_iB-BA_i)(BA_i^t - A_i^tB)\\ = &A_iB^2A_i^t-A_iBA_i^tB-BA_iBA_i^t+BA_iA_i^tB\end{aligned}$$
And we sum up over $i$, $$(\sum_{i} A_iB^2A_i^t)-B^2-B^2+B^2= (\sum_{i} A_iB^2A_i^t)-B^2$$ After taking the trace, $$\text{Tr}((\sum_{i} A_iB^2A_i^t)-B^2)=\text{Tr}(\sum_i B^2A_i^tA_i-B^2)=\text{Tr}(B^2(\sum_i A_i^tA_i-I))$$
Now since $B$ and hence $B^2$ is positive definite, and if we further assume $\sum_i A_i^tA_i\le I$, then $B^2(I-\sum_i A_i^tA_i)$ must be
positive semi-definite. Hence $$\text{Tr}(B^2(\sum_i A_i^tA_i-I))=-\text{Tr}(B^2(I-\sum_i A_i^tA_i))\le 0$$
But $\text{Tr}(B^2(\sum_i A_i^tA_i -I))$ is a sum of Frobenius norm squared, hence nonnegative, so it must be $0$, and each individual norm must be $0$.
A: The statement is false. Here is a counterexample:
$$
A_1=\pmatrix{\frac{1}{\sqrt{2}}&0&0\\ \frac{1}{\sqrt{6}}&\frac{1}{\sqrt{6}}&\frac{1}{\sqrt{6}}\\ 0&0&\frac{1}{\sqrt{2}}},
\,A_2=\pmatrix{\frac{1}{\sqrt{2}}&0&0\\ -\frac{1}{\sqrt{6}}&\frac{1}{\sqrt{6}}&-\frac{1}{\sqrt{6}}\\ 0&0&\frac{1}{\sqrt{2}}},
\,B=\pmatrix{0\\ &\frac{1}{3}\\ &&\frac{2}{3}}.
$$
One may directly verify that $A_1A_1^T+A_2A_2^T=I$ and $A_1BA_1^T+A_2BA_2^T=B$. However, since diagonal matrices with distinct diagonal entries only commute with diagonal matrices, $B$ does not commute with any of the $A_i$s.
In the counterexample above, $B$ has three distinct eigenvalues. The statement in question is true if $B$ has only two distinct eigenvalues. In this case, by a change of orthonormal basis we may assume that $B=\lambda_1I_{n_1}\oplus\lambda_2I_{n_2}$ (where $\lambda_1>\lambda_2$) and
$$
A_i=\pmatrix{A_{11}^{(i)}&A_{12}^{(i)}\\ A_{21}^{(i)}&A_{22}^{(i)}}.
$$
By the given conditions, we have
$$
\sum_i(tI-B)^{-1/2}A_i(tI-B)A_i^T(tI-B)^{-1/2}=I
$$
whenever $t>\lambda_1$. Thus
$$
\begin{align}
n&=\operatorname{tr}\left(\sum_i(tI-B)^{-1/2}A_i(tI-B)A_i^T(tI-B)^{-1/2}\right)\\
&=\sum_i\left\|(tI-B)^{-1/2}A_i(tI-B)^{1/2}\right\|_F^2\\
&=\sum_{i=1}^m\left(
\left\|A_{11}^{(i)}\right\|_F^2
+\left\|A_{12}^{(i)}\right\|_F^2\frac{t-\lambda_2}{t-\lambda_1}
+\left\|A_{21}^{(i)}\right\|_F^2\frac{t-\lambda_1}{t-\lambda_2}
+\left\|A_{22}^{(i)}\right\|_F^2
\right)\\
&=\sum_{i=1}^m\left(
\left\|A_{11}^{(i)}\right\|_F^2
+\left\|A_{12}^{(i)}\right\|_F^2
+\left\|A_{21}^{(i)}\right\|_F^2
+\left\|A_{22}^{(i)}\right\|_F^2
\right)\\
&\quad+\left(\sum_{i=1}^m\left\|A_{12}^{(i)}\right\|_F^2\right)\frac{\lambda_1-\lambda_2}{t-\lambda_1}
+\left(\sum_i\left\|A_{21}^{(i)}\right\|_F^2
\right)\frac{\lambda_2-\lambda_1}{t-\lambda_2}\\
\end{align}
$$
on the interval $(\lambda_1,\infty)$. Hence $\sum_{i=1}^m\left\|A_{12}^{(i)}\right\|_F^2$ and $\sum_{i=1}^m\left\|A_{21}^{(i)}\right\|_F^2$ must be zero, i.e., $A_{12}^{(i)}$ and $A_{21}^{(i)}$ must be zero for each $i$. This means each $A_i$ is a block-diagonal matrix that commutes with $B$.
A: Denote by $\phi : \Bbb{R}^{n\times n} \to \Bbb{R}^{n\times n}$ the map $\phi(X) = \sum_{i=1}^m A_iXA_i^t$ where $\sum_{i=1}^m A_iA_i^t = I$.
Assume $\phi(B) = B$ for some matrix $B \ge 0$ (the assumption $\operatorname{Tr}B = 1$ is not necessary, as is evident from the fact that we can scale both sides of the equality by an arbitrary nonnegative scalar). The proof will be presented as a long series of claims:

*

*$\phi$ decreases trace, i.e. for all $X \in \Bbb{R}^{n\times n}$ we have $\operatorname{Tr}\phi(X) \le \operatorname{Tr} X$. (This works only if $\sum_{i=1}^n A_i^tA_i \le I$ which doesn't follow from $\sum_{i=1}^n A_iA_i^t = I$. Therefore this proof is incomplete.).


*For $X \ge 0$ we have $\phi(X) \ge 0$. Furthermore, if $\phi(X) \ge X$ then $\phi(X) = X$.
The latter follows from the fact that $\phi(X)-X$ is a positive matrix with trace $\le 0$, hence it is equal to $0$.

*

*As Ben Grossman suggests, if we define the block matrix $A = [A_1 \cdots A_n] \in \Bbb{R}^{n\times mn}$, we have $AA^t = I_{n\times n}$ and $$\phi(X) = A\operatorname{diag}(X,\ldots, X)A^t, \qquad \text{ for all }X \in \Bbb{R}^{n\times n}.$$
Indeed, we have
\begin{align}
A\operatorname{diag}(X,\ldots, X)A^t &= \begin{bmatrix} A_1 & \cdots & A_m\end{bmatrix}\begin{bmatrix} X & \cdots & 0 \\ \vdots & \ddots & \vdots \\
0 & \cdots & X\end{bmatrix} \begin{bmatrix} A_1^t \\ \vdots \\ A_m^t\end{bmatrix}\\
&= \begin{bmatrix} A_1X & \cdots & A_mX \end{bmatrix}\begin{bmatrix} A_1^t \\ \vdots \\ A_m^t\end{bmatrix}\\
&= \sum_{i=1}^m A_iXA_i^t\\
&= \phi(X).
\end{align}

*

*We have $A^tA \le I_{mn \times mn}$.

Indeed, $A^tA \ge 0$ and $\sigma(A^t A) = \sigma(AA^t) \cup \{0\} = \{0,1\}$ so the claim follows.

*

*For all $k \in \Bbb{N}_0$ we have $\phi(B^k) = B^k$.

Indeed, we have
\begin{align}
B^k &= \phi(B)^k\\
&= (A\operatorname{diag}(B,\ldots,B)A^t)^k\\
&= A\operatorname{diag}(B,\ldots, B)\underbrace{A^tA}_{\le I}\operatorname{diag}(B,\ldots, B)\underbrace{A^tA}_{\le I}\operatorname{diag}(B,\ldots, B)A^t \cdots A\operatorname{diag}(B,\ldots, B)A^t\\
&\le A\operatorname{diag}(B,\ldots, B)^kA^t\\
&= A\operatorname{diag}(B^k\ldots, B^k)A^t\\
&= \phi(B^k).
\end{align}
so in fact $\phi(B^n)= B^k$.

*

*By linearity of $\phi$ it follows that $\phi(p(B)) = p(B)$ for every polynomial $p \in \Bbb{R}[z]$.


*For every orthogonal eigenprojection $P$ of $B$ we have $\phi(P) = P$.
Indeed, $P$ can be expressed as a real polynomial in $B$.

*

*For every orthogonal eigenprojection $P$ of $B$ we have $PA_i = A_iP$ for all $1 \le i \le m$.

Indeed, we have
$$0 = (I-P)P(I-P) = (I-P)\phi(P)(I-P) = \sum_{i=1}^m \underbrace{((I-P)A_i) P ((I-P)A_i)^t}_{ \ge 0}$$
so for all $1 \le i \le m$ it follows
$$((I-P)A_i) P ((I-P)A_i)^t = 0 \implies ((I-P)A_iP)  ((I-P)A_iP)^t = 0 \implies (I-P)A_iP = 0.$$
The same calculation for the eigenprojection $I-P$ shows $PA_i(I-P) = 0$ so we finally get $PA_i = A_iP$.

*

*We have $BA_i = A_iB$ for all $1 \le i \le m$.

This follows from the fact that $B$ is a real linear combination of its eigenprojections.
A: With the additional condition that $\sum_iA_i^TA_i\preceq I$, the statement is true. It is trivial when $B$ is a scalar multiple of $I$. Now suppose $B$ has at least two different eigenvalues. By a change of orthonormal basis, we may assume that
$$
B=\pmatrix{\lambda I\\ &D}
\quad\text{and}\quad
A_i=\pmatrix{X_i&Y_i\\ Z_i&W_i}
$$
where $D\prec\lambda I$. By the given conditions, we have $\sum_iA_i(\lambda I-B)A_i^T=\lambda I-B$. Hence the diagonal sub-blocks on one side are equal to those on the other, i.e.,
$$
\begin{align}
&\sum_iY_i(\lambda I-D)Y_i^T=0,\tag{1}\\
&\sum_iW_i(\lambda I-D)W_i^T=\lambda I-D,\tag{2}\\
\end{align}
$$
As $\lambda I-D\succ0$, $(1)$ implies that $Y_i=0$, while $(2)$ implies that
$$
\operatorname{tr}\left((\lambda I-D)^{1/2}\left(I-\sum_iW_i^TW_i\right)(\lambda I-D)^{1/2}\right)=0.\tag{3}
$$
However, by the additional condition $\sum_iA_i^TA_i\preceq I$, we see that $I-\sum_iW_i^TW_i$ is positive semidefinite. Therefore $(3)$ implies that $I-\sum_iW_i^TW_i=0$. In turn, the condition $\sum_iA_i^TA_i\preceq I$ implies that $Z_i=0$. Hence each $A_i$ is a block-diagonal matrix. Proceed recursively, we see that if we write $B=\lambda_1 I\oplus\lambda_2I\oplus\cdots\oplus\lambda_k I$ where the $\lambda_i$s are the distinct eigenvalues of $B$, then each $A_i$ must be a block-diagonal matrix with the same diagonal block sizes as $B$. Hence $A_i$ commutes with $B$.
The above proof only requires that $B$ is symmetric. Positive semidefiniteness of $B$ and the trace condition $\operatorname{tr}(B)=1$ are not needed.
A: I assume $B$ "is positive" means symmetric positive definite.  For clarity we may assume each $A_k \neq \mathbf 0$ (otherwise just remove the zero matrices from the sum and re-index).
Note: OP's subsequently added constraint that
$I-\sum_{k=1}^m A_k^T A_k\succeq \mathbf 0$
$\implies  I= \sum_{k=1}^m A_k^T A_k $
proof: we need to check that the trace of this symmetric PSD matrix is zero (since the trace gives the Frobenius norm of said matrix's square root) and
$\text{trace}\big(I-\sum_{k=1}^m A_k^T A_k\big) $
$= \text{trace}\big(I\big)-\sum_{k=1}^m \text{trace}\big(A_k^T A_k\big)$
$=n-\sum_{k=1}^m \text{trace}\big(A_k A_k^T\big)$
$=n-\text{trace}\big(\sum_{k=1}^m A_k A_k^T\big)$
$= 0$
since we are told $\sum_{k=1}^m A_k A_k^T=I$

Now for the proof:
$B=\sum_{k=1}^m A_k B A_k^T \implies A_kB =BA_k$
The technique is to compute the squared Frobenius norm of the difference, then find it implies Cauchy-Schwarz is met with equality and exploit the associated equality conditions.
$0= \big\Vert B-\sum_{k=1}^m A_k D A_k^T \big\Vert_F^2 =  2 \cdot \text{trace}\Big(B^2\Big) -\text{trace}\Big(B\big(\sum_{k=1}^m A_k BA_k^T\big) \Big)-\text{trace}\Big(\big(\sum_{k=1}^m A_k BA_k^T\big)B \Big)$
$\implies \text{trace}\big(B^2\big) =\text{trace}\big(B\sum_{k=1}^m A_k BA_k^T \big) = \sum_{k=1}^m \text{trace}\big(BA_k BA_k^T \big)$

Now consider individual $A_k$
$\text{trace}\big(BA_k BA_k^T \big)=\text{trace}\big((BA_k)( BA_k^T) \big)=\text{trace}\big( (BA_k^T)(BA_k) \big)=\text{trace}\big( (A_kB)^T(BA_k) \big)$
$\leq \big\Vert A_k B\big\Vert_F\cdot \big\Vert BA_k\big\Vert_F$
$\leq \frac{1}{2}\Big(\big\Vert A_k B\big\Vert_F^2+ \big\Vert BA_k\big\Vert_F^2\Big)$
$= \frac{1}{2}\Big( \text{trace}\big( B^2 A_k^T A_k \big) +  \text{trace}\big(   B^2 A_kA_k^T \big)\Big)$
$= \frac{1}{2}\cdot \text{trace}\Big( B^2 \big( A_k^T A_k + A_kA_k^T \big)\Big)$
where the inequalities are (i) Cauchy-Schwarz and (ii) $\text{GM}\leq \text{AM}$.
Summing over $k$:
$\sum_{k=1}^m \text{trace}\big(BA_k BA_k^T \big)$
$\leq \sum_{k=1}^m  \frac{1}{2}\cdot \text{trace}\Big( B^2 \big( A_k^T A_k + A_kA_k^T \big)\Big)$
$=  \frac{1}{2}\cdot \text{trace}\Big(\sum_{k=1}^m B^2  \big( A_k^T A_k + A_kA_k^T \big)\Big)$
$=  \frac{1}{2}\cdot \text{trace}\Big( B^2 \sum_{k=1}^m \big( A_k^T A_k + A_kA_k^T \big)\Big)$
$=  \frac{1}{2}\cdot \text{trace}\Big( B^2 \big( 2 I \big)\Big)$
$=\text{trace}\big( B^2\big)$
And we know $\text{trace}\big(B^2\big) = \sum_{k=1}^m \text{trace}\big(BA_k BA_k^T \big)$ from the earlier squared Frobenius norm calculation, i.e. the upper bound is met with equality .
The equality conditions of (i) tell us
$ A_k B = \sigma_k\cdot  BA_k$
left multiplying by $A_k^T$ and taking the trace
$0\lt \text{trace}\big(A_k^T A_k B\big) = \sigma_k\cdot \text{trace}\big(A_k^T BA_k\big) = \sigma_k\cdot \text{trace}\big(A_k A_k^T B\big)\implies \sigma_k \gt 0$
since a symmetric PSD matrix times a symmetric PD matrix has positive trace (which is equal to the Frobenius norm of the product of their square roots and the PD matrix square root is invertible so said product is non-zero)
Combining this with the equality conditions of (ii) we have
$0\lt \big\Vert BA_k\big\Vert_F=\big\Vert A_k B\big\Vert_F=  \big\Vert\sigma_k\cdot  BA_k\big\Vert_F =\sigma\cdot \big\Vert BA_k\big\Vert_F\implies \sigma =1$
conclude $A_k B = 1 \cdot B A_k$
i.e. the desired commuting relationship holds
