Let $\{Y\_i\}\_{i=1}^N\in\mathbb{R}^{n\times m}$ be a set of full column rank matrices ($\mathrm{rank}(Y\_i)=m$ for all $i$) and $\{P\_i\}\_{i=1}^N\in\mathbb{R}^{m\times m}$ be a set of *positive definite* matrices. Let $A\in\mathbb{R}^{m\times n}$, $B\in\mathbb{R}^{m\times n}$ and consider the following equation
$$\tag{1}\label{eq:1}
\sum\_{i=1}^N Y\_i Y\_i^\top \Delta\_i Y\_i Y\_i^\top=0\_{n},
$$
where
$$
\Delta\_i := A^\top P\_i B+B^\top P\_i A,
$$
$0\_{n}$ denotes the $n\times n$ zero matrix and $(\cdot)^\top$ denotes transposition.



> 
> **My question.** Suppose that either $A$ or $B$ is of full (row) rank. Does \eqref{eq:1} imply $Y\_i^\top \Delta\_i Y\_i=0\_{m}$ for all $i=1,2,\dots,N$?
> 
> 
> 




---


**Two remarks.** 


1. The answer is in the affirmative if $P\_i=p\_i M$, for all $i$, with $M\in\mathbb{R}^{m\times m}$ being a positive definite matrix and $\{p\_i\}\_{i=1}^N$ being a set of positive real numbers. This follows from the fact that \eqref{eq:1} is equivalent to 
$$\tag{2}\label{eq:2}
\sum\_{i=1}^N \mathrm{tr}(Y\_i^\top XY\_i Y\_i^\top \Delta\_i Y\_i) = 0, \ \ \text{for all symmetric } X\in\mathbb{R}^{n\times n},
$$
where $\mathrm{tr}(\cdot)$ denotes the trace operator. More precisely, if we pick $X=A^\top MB+B^\top M A$ in the previous expression, the LHS of \eqref{eq:2} turns into a sum positive numbers which implies that $Y\_i^\top \Delta\_i Y\_i=0\_{m}$ for all $i=1,2,\dots,N$, as desired.
2. The assumption that either $A$ or $B$ is of full (row) rank seems to be crucial. Indeed, if we drop this assumption, then the answer is in the negative. For a counterexample pick $N=2$ and
$$
Y\_1=Y\_2=I\_2,\ P\_1=\begin{bmatrix} 2&1\\1&2\end{bmatrix},\ P\_2=\begin{bmatrix} 2&-1\\-1&2\end{bmatrix},\ A=\begin{bmatrix} 0&0\\1&0 \end{bmatrix},\ B=\begin{bmatrix}1&0\\0&0\end{bmatrix}.
$$



On a condition for a matrix sum to be zero
EDIT: Now with an embarrassingly simple counterexample. 


Same as my other integer counterexample, but now with very small magnitude integer entries in all matrices. $A$ and $B$ full rank, all $Y\_i'$ s are identity matrices, $m = n = N = 2$.



```
>> disp(P1)
     1     0
     0     1
>> disp(P2)
     1     0
     0     2
>> disp(A)
     3     0
     0     2
>> disp(B)
     0     1
    -1     0
>> Delta1 = A'*P1*B+B'*P1*A
Delta1 =
     0     1
     1     0
>> Delta2 = A'*P2*B+B'*P2*A
Delta2 =
     0    -1
    -1     0

```

+++++++++++++++++++++++++++++++++++++++++++


**EDIT: New exact counterexample immediately below:**


At the request of the OP, I am providing an exact counterexample with all matrix elements being integer. As with my other counterexample,s $A$ and $B$ are both full rank, and as with my second counterexample, all $Y\_i's$ are identity matrices.


$m = n = N = 2$.



```
>> disp(P1)
     1     0
     0     1
>> disp(P2)
     1     0
     0     2
>> disp(A)
    60    42
    44    38
>> disp(60*38-42*44)
   432
>> disp(B)
    22    19
   -20   -14
>> disp(22*(-14)-(-20)*19)
    72
>> Delta1 = A'*P1*B+B'*P1*A
Delta1 =
   880   688
   688   532
>> Delta2 = A'*P2*B+B'*P2*A
Delta2 =
  -880  -688
  -688  -532

```

+++++++++++++++++++++++++++++++++++++++++++
I am editing this post to add another counterexample which shows counterexample exists even if the $Y\_i$ are restricted to identity matrices. My original counterexample still stands, and is at the end (I also added singular values of Y1 and Y2, to demonstrate full rank). So as to make these counterexamples reproducible, I use exactly the matrices, as displayed.


The new counterexample is with $m = n = N = 2$, is full rank for $A$ and $B$, as shown by the singular values output from svd, and uses all $Y\_i$ = identity matrix. 


MATLAB output:



```
>> disp(A)
   0.110169472772473   0.046128343729780
   0.046222313869941  -0.019099888931544
>> disp(svd(A))
   0.124969084324384
   0.033899452001928
>> disp(B)
   0.091136400531647  -0.038150840712963
  -0.224387878673384  -0.094158263094596
>> disp(svd(B))
   0.253801784511592
   0.067540231952775
>> disp(P1)
  10.912674303492269  -1.179600975604882
  -1.179600975604882   9.910545724120434
>> disp(eig(P1))
   9.130000133988611
  11.693219893624091
>> disp(P2)
  10.173059165015179   1.124128224582275
   1.124128224582275  10.611386694428143
>> disp(eig(P2))
   9.246929508338818
  11.537516351104504
>> Delta1=A'*P1*B+B'*P1*A
Delta1 =
   0.061940244090696   0.027931132691275
   0.027931132691275   0.005765260330663
>> Delta2=A'*P2*B+B'*P2*A
Delta2 =
  -0.061940244090696  -0.027931132691276
  -0.027931132691276  -0.005765260330665
>> Delta1+Delta2
ans =
   1.0e-14 *
   0.037470027081099  -0.134614541735800
  -0.131838984174237  -0.230718222304915

```

As can be seem, Delta1 and Delta are non-zero, and their sum is the zero matrix to within numerical tolerance




---


Original counterexample:


Here i a counterexample with $n = m = N = 2$. So as to make it reproducible, I use exactly the matrices, as displayed. As a result, the required sum, `Y1*Y1'*Delta1*Y1*Y1'+Y2*Y2'*Delta2*Y2*Y2'` , has max element deviation from zero of 1e-13. Using my "internal" counterexample (i.e., not rounded to 16 digits displayed), the sum has a max element deviation of 2e-16, so I believe it is legitimate.


Output from MATLAB follows:


svd displays the singular values, demonstrating that $A$ and $B$ are both of full rank.



```
>> disp(A)
   0.041926790622772   0.023648187924399
   0.027723496832564   0.006701524620265
>> disp(svd(A))
   0.055543689989045
   0.006744907360037
>> disp(B)
  -0.034363733684070   0.000684366074426
   0.046234868075085   0.017527803003490
>> disp(svd(B))
   0.059260236050568
   0.010697937974866
>> disp(P1)
   5.545624992234104   1.061540046345759
   1.061540046345759   6.241963078368660
>> disp(eig(P1))
   4.776615019858404
   7.010973050744360
>> disp(P2)
   6.237907041799323  -0.887038870822917
  -0.887038870822917   5.003187529900254
>> disp(eig(P2))
   4.539820203457168
   6.701274368242409
>> disp(Y1)
  -3.152948589228529   1.160433382904003
   1.308458703988197  -3.430528144820526
>> disp(svd(Y1))
   4.534793991971979
   2.050346659786896
>> disp(Y2)
  -0.864647428071400  -0.271153534587496
  -4.073950010112803  -4.647399186185172
>> disp(svd(Y2))
   6.228780751866293
   0.467779477983433
>> Delta1=A'*P1*B+B'*P1*A
Delta1 =
   0.002114883742232   0.002336175771742
   0.002336175771742   0.002535655871380
>> Delta2=A'*P2*B+B'*P2*A
Delta2 =
  -0.006897453634519  -0.002343057693250
  -0.002343057693250   0.000633794608936
>> Y1*Y1'*Delta1*Y1*Y1'+Y2*Y2'*Delta2*Y2*Y2'
ans =
   1.0e-12 *
  -0.013097162243625  -0.023654689318420
  -0.023522850334246  -0.132810429320784
>> Y1'*Delta1*Y1
ans =
   0.006089643801744   0.009696210846174
   0.009696210846174   0.014088678339379
>> Y2'*Delta2*Y2
ans =
  -0.011144492247550  -0.001620866985435
  -0.001620866985435   0.007276519088282

```

When I use my "internal" $A$, which has deviation from the above provided A of 



```
   1.0e-15 *
   0.228983498828939   0.409394740330526
  -0.423272528138341  -0.196023752785379

```

I get 



```
>> Y1*Y1'*Delta1*Y1*Y1'+Y2*Y2'*Delta2*Y2*Y2'
ans =
   1.0e-15 *
  -0.117961196366423   0.242861286636753
   0.159594559789866  -0.027755575615629

```

