Product $PVPVP$ is elementwise nonnegative? 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?
 A: Although you said you already proved the result for $n=2$, perhaps it's worth recording a proof here.
The statement is clear if both diagonal entries of $V$ have the same sign, so assume that $V = \left(\matrix{v_1&0\\ 0&-v_2}\right)$ with $v_1\ge0$ and $v_2\ge0$.  If $P = \left(\matrix{a&b\\ c&d}\right)$ then by direct computation,
$$PVPVP = \left(\matrix{a^3v_1^2 - bc(2av_1v_2 - dv_2^2)&b(a^2v_1^2-(ad+bc)v_1v_2 + d^2v_2^2)\\
 c(a^2v_1^2-(ad+bc)v_1v_2 + d^2v_2^2)& d^3v_2^2 - bc(2dv_1 v_2 - av_1^2 )}\right).$$
Since the determinant of $P$ is positive, $-bc \ge -ad$, so
$$a^2v_1^2-(ad+bc)v_1v_2 + d^2v_2^2 \ge a^2v_1^2-2adv_1v_2 + d^2v_2^2 = (av_1 - dv_2)^2 \ge 0.$$
Therefore the off-diagonal entries of $PVPVP$ are nonnegative.
If $2av_1v_2 - dv_2^2<0$ then the $(1,1)$ entry of $PVPVP$ is clearly nonnegative; otherwise,
$$a^3v_1^2 - bc(2av_1v_2 - dv_2^2)\ge a^3v_1^2 - ad(2av_1v_2 - dv_2^2) =a(av_1 -dv_2)^2 \ge 0. $$
Similarly, if $2dv_1 v_2 - av_1^2 < 0$ then the $(2,2)$ entry is clearly nonnegative; otherwise
$$d^3v_2^2 - bc(2dv_1 v_2 - av_1^2 )\ge d^3v_2^2 - ad(2dv_1 v_2 - av_1^2 ) = d(dv_2 - av_1)^2 \ge 0. $$
A: First, we repeat the arguments from this stackexchange answer.  $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.]
