Let $P\in \mathbb{R}^{n\times n}$ be the inverse of a positive definite M-matrix and $V\in \mathbb{R}^{n\times n}$ be any diagonal matrix. Prove (or disprove) that $PVPVP$ is elementwise nonnegative. 


I know of the following: 


$P$ is positive definite and elementwise nonnegative. Moreover, $p\_{jk}p\_{ii} \ge p\_{ji}p\_{ik}$ for any $i,j,k$. 


I can verify that the statement is true for $n=2$, but I don't know how to work with $n$ large. Playing around with randomly generated matrices in Matlab seems to suggest that the statement is true. Any hint or suggestion would be greatly appreciated. 


I've googled out that a very similar statement was put as a conjecture in this paper: *Optimization of an on-chip active cooling system based on thin-film thermoelectric coolers* (<http://dl.acm.org/citation.cfm?id=1870955>)


**Edit:** Perhaps someone can solve this easier question: Is there a positive semi-definite and elementwise nonnegative $P$ and diagonal $V$ such that $PVPVP$ is not elementwise nonnegative?



Product $PVPVP$ is elementwise nonnegative?
First, we repeat the arguments from this [stackexchange answer](https://math.stackexchange.com/a/989748). $P^{-1}$ is an $M$-matrix, and can thus be written as $s(I-A)$ for some positive $s$ and some $A$ with non-negative entries. As $P^{-1}$ is positive definite, the spectrum of $A$ lies to the left of $\{ z: \hbox{Re}(z) = 1 \}$, and hence by Perron-Frobenius the spectral radius of $A$ is less than $1$. Thus we have the absolutely convergent Neumann series
$$ P = s^{-1} (I + A + A^2 + \dots )$$
and hence
$$ PVPVP = s^{-3} \sum\_{i=0}^\infty \sum\_{j=0}^\infty \sum\_{k=0}^\infty A^i V A^j V A^k.$$
It thus suffices to show that
$$ \sum\_{i+j+k=m} A^i V A^j V A^k \quad (1)$$
has non-negative coefficients for each $m \geq 0$ (where $i,j,k$ are understood to be non-negative integers). By change of variables, this is
$$ \sum\_{0 \leq q \leq r \leq m} A^q V A^{r-q} V A^{m-r}.$$
Writing $A = (a\_{st})\_{1 \leq s,t \leq n}$ and $V = \hbox{diag}(v\_1,\dots,v\_n)$, the $st$ coefficient of (1) can be expanded as
$$ 
\sum\_{s=s\_0,s\_1,\dots,s\_m=t} a\_{s\_0 s\_1} \dots a\_{s\_{m-1} s\_m}\sum\_{0 \leq q \leq r \leq m} v\_{s\_q} v\_{s\_r}.$$
But the quadratic form
$$ \sum\_{0 \leq q \leq r \leq m} x\_q x\_r = \frac{1}{2}(x\_0+\dots+x\_m)^2 + \frac{1}{2} x\_0^2 + \dots + \frac{1}{2} x\_m^2$$
is positive definite, and the $a\_{st}$ are non-negative, and the claim follows.


[For the record, I found this argument while performing a perturbative analysis in the case where $P$ was close to $I$, or more precisely $P = (I-A)^{-1}$ for some $A$ with small non-negative entries.]


