Apply doubly stochastic matrix M to a probability vector, then entropy increases? Consider a vector $p =(p_1,...p_n)$, $p_i>0$, $\sum p_i = 1$
and a matrix $M_{ij}$, which is doubly stochastic: $\sum_i M_{ij} = 1,  \sum_j M_{ij} = 1, M_{ij} > 0$.
Question 1  Just apply matrix M to a vector $p$ , i.e.  $q = Mp$ (i.e. $q_i = \sum_j M_{ij} p_j$)  is it true that entropy of new vector $q$ is greater  then of original vector $p$ ? I.e.
$$ H(q) = -\sum_i q_i ln(q_i) > -\sum_i p_i ln(p_i) = H(p) $$ 
Question 2 Is there a simple proof of it ? (It might follow from the Gibb's inequality, but it does not seem obvious for me). 
Question 3 What are generalizations,
in view of meta-principle "entropy always grows" that might be an example of some more general phenomena ? 

Motivation: There is a meta-principle that entropy grows is some natural systems, matrix applied to a vector is probably the most simple system one can consider (Markov chain), so the question above arises. After some thinking one can restrict from all matrices to doubly stochastic, because $M^n v$ tends to a uniform distribution (which has maximal entropy of all) only for doubly stochastic $M$ , so it cannot be true for general $M$. 
 A: This is a particular case of the following general principle. On the one hand, a bistochastic matrix $M$ has the property that for every non-negative vector (say, a probability vector), $Mp\succ p$. On another hand, the order $\succ$ can be defined by $x\succ y$ iff $\sum_if(x_i)\le\sum_if(y_i)$ for every convex function $f$.
Now, use the fact that $f=-H$ is convex.
For proofs, see my book Matrices (Springer-Verlag GTM #216, second edition) at sections 6.5 (Proposition 6.4) and 8.5 (Theorem 8.5).
Edit (at Alexander's request.) Let $\alpha\ne1$ be a positive parameter. The Renyi entropy is
$$\frac1{1-\alpha}\log\sum_jp_j^\alpha.$$
The map $p\mapsto Mp$ does increase Renyi entropy. Proof: if $\alpha>1$, the map $t\mapsto t^\alpha$ is convex, and $s\mapsto\frac1{1-\alpha}\log s$ is increasing. If instead $\alpha<1$, the map $t\mapsto t^\alpha$ is concave, and $s\mapsto\frac1{1-\alpha}\log s$ is decreasing.
A: The answer to Question 1 is yes as long as $p$ is not the uniform distribution $(\frac{1}{n},\frac{1}{n},\ldots,\frac{1}{n})$.
The proof (Question 2) is quite simple:
Recall from the Birkhoff–von Neumann theorem that every doubly stochastic matrix $M$ is a convex combination of permutation matrices.  We can interpret this as saying that there is a distribution $\theta$ on the set of all permutations of $A:=\{1,2,\ldots,n\}$ such that whenever $\mathbf{x}$ is a random variable from $A$ and $\pmb{\pi}$ is a random permutation of $A$ with distribution $\theta$ independent of $\mathbf{x}$, and we set $\mathbf{y}:=\pmb{\pi}(\mathbf{x})$, then we have
\begin{align}
   \mathbb{P}(\mathbf{y}=i\,|\,\mathbf{x}=j) &= M_{i,j} \;.
\end{align}
Note that if $\mathbf{x}$ is distributed as $p$, then $\mathbf{y}$ is distributed as $q:=Mp$, so $H(\mathbf{x})=H(p)$ and $H(\mathbf{y})=H(q)$.
Now,
\begin{align}
   H(\mathbf{y},\pmb{\pi}) &= H(\mathbf{y}) + H(\pmb{\pi}\,|\,\mathbf{y}) \;, \\
   H(\mathbf{y},\pmb{\pi}) &= H(\pmb{\pi}) + \underbrace{H(\mathbf{y}\,|\,\pmb{\pi})}_{H(\mathbf{x})} \;,
\end{align}
which implies
\begin{align}
   H(\mathbf{y}) &= H(\mathbf{x}) + H(\pmb{\pi}) - H(\pmb{\pi}\,|\,\mathbf{y}) \\
   &= H(\mathbf{x}) + I(\mathbf{y};\pmb{\pi})
\end{align}
where $I(\mathbf{y};\pmb{\pi})$ is the mutual information between $\mathbf{y}$ and $\pmb{\pi}$.  We know that $I(\mathbf{y};\pmb{\pi})\geq 0$ with equality if and only if $\mathbf{y}$ and $\pmb{\pi}$ are independent, which is the case if and only if $\mathbf{x}$ has the uniform distribution. Q.E.D.
For Question 3, suppose that $M$ is merely a positive stochastic matrix (not doubly stochastic).  Let $r$ denote the unique stationary distribution of $M$.  Then we have a similar ``entropy increase'' principle if we replace entropy with (minus) the Kullback-Leibler divergence relative to $r$.  There is a nice discussion on this in the book of Cover and Thomas.
A: This is a particular case of a well-known monotonicity property of the Kullback-Leibler divergence $D(\cdot\|\cdot)$. Namely, if $\alpha$ and $\beta$ are any two distributions on the same (say, finite) space $X$, and $P$ is a Markov operator on $X$, then 
$$
D(\alpha\|\beta) \ge D(\alpha P\|\beta P) \;.
$$
In your case double stochasticity of $P$ means that its stationary distribution is the uniform distribution $m_X$ on $X$, whereas 
$$
D(\alpha \| m_X) = \log \text{card} X - H(\alpha) \;.
$$
A: This is true, except that we have equality when $p$ is the uniform distribution, which is the limit distribution of your matrix. (Proof in here, 15.5, for example.)
Whenever $p_i\geq 1/n$, it is not so hard to show that $q_i\leq p_i$, and when $p_i\leq 1/n$, we have $q_i\geq p_i$. From there, it follows that $-\sum_i p_i \ln (q_i) \leq -\sum_i q_i \ln (q_i) $, and this gives the result combined with Gibbs'.
