Response of stationary distribution to perturbation of a stochastic matrix

Question

Suppose I have a (left) stochastic matrix $P$, i.e. a non-negative matrix with column sums equal to 1. Its maximum eigenvalue will be equal to 1, and the corresponding eigenvector $\mathbf q$, if it's unique, can be interpreted the stationary distribution of a stochastic process with transition probabilities $P$.

Now suppose I have a family of stochastic matrices, $P(x)$, characterised by a real parameter $x$. I would like to know how the leading eigenvector $\mathbf{q}$ changes in response to a small change in $x$. In other words, given the values of $p_{ij}$ and $\frac{dp_{ij}}{dx}$, I would like to know $\frac{dq_i}{dx}$.

I'm aware that there is a general method for calculating the derivatives of eigenvectors. However, it's a little bit complicated, and I get the feeling it should simplify in this case, since we know that the leading eigenvalue must remain equal to 1 as the matrix changes.

Since $\mathbf{q}(x)$ is to be interpreted as a probability distribution, it's advantageous to assume it's normalised as $\sum_i q_i=1$, rather than $\mathbf{q}^T\mathbf{q}=1$.

I'm happy to accept answers that assume the leading eigenvector is unique (i.e. $P$ is primitive), but if there's an elegant way to handle the other cases, that would be great.

Answer 1

Assume that $P(x)$ has only one eigenvalue equal to $1$. Let $q(x)$ be the unique unit-norm leading eigenvector of $P(x)$ for all $x$. For convenience, I'll drop the "$(x)$" in the computations.

We require $Pq = q$, i.e. $(I-P)q = 0$ and $q^Tq = 0$ for all $x$.

Differentiation yields: $(I-P)\dfrac{dq}{dx} - \dfrac{dP}{dx}q = 0$ and $q^T\dfrac{dq}{dx} + \dfrac{dq}{dx}^Tq = 0$,

These can be rewritten as (1) $(I-P)\dfrac{dq}{dx} = \dfrac{dP}{dx}q$ and (2) $q^T\dfrac{dq}{dx} = 0$.

By assumption, $P$ has only one eigenvalue equal to $1$, so $I-P$ has only one $0$ eigenvalue. Thus, $I-P$ has a one dimensional nullspace, namely $\text{span}(q)$. So, if we can find one solution $\dfrac{dq}{dx} = v$ to (1), then all solutions to (1) will be of the form $\dfrac{dq}{dx} = v+tq$ for some $t \in \mathbb{R}$. Using the pseudoinverse, one solution is $\dfrac{dq}{dx} = (I-P)^+\dfrac{dP}{dx}q$. Therefore, the solutions to (1) are all in the form $\dfrac{dq}{dx} = (I-P)^+\dfrac{dP}{dx}q + tq$ for some $t \in \mathbb{R}$.

Condition (2) requires that $\dfrac{dq}{dx}$ is orthogonal to $q$. Since the solution $(I-P)^+\dfrac{dP}{dx}q$ is already orthogonal to the nullspace of $I-P$, which includes $q$, this is the only solution of (1) which also satisfies (2).

Therefore, the solution to (1) and (2) is: $\dfrac{dq}{dx} = (I-P)^+\dfrac{dP}{dx}q$.

If we instead wish to normalize $q$ such that $\sum_{i}q_i = 1$ instead of $q^Tq = 1$, then condition (2) becomes $\displaystyle\sum_{i}\dfrac{dq_i}{dx} = 0$. The solution will still be in the form $\dfrac{dq}{dx} = (I-P)^+\dfrac{dP}{dx}q + tq$ for some $t \in \mathbb{R}$. However, finding the value of $t$ to satisfy (2) might be harder.

Answer 2

Let me proceed the solution of JimmyK4542. He is absolutely correct and let me just prove that the solution indeed $$ q^\prime(x)=(I−P(x))^+P′(x)q(x) $$ Notice that the left eigenvector for $P(x)$ corresponding to 1 is the vector of ones. Let's denote such vector (normalized by 1) as $e$. Thus $$ e^T (I-P(x))q(x)^\prime= e^T P^\prime(x) q(x)=0 $$ Thus $$ P^\prime(x) q(x)=(I-ee^T)P^\prime(x) q(x) $$ But $$ (I-ee^T)=(I-P(x))(I-P(x))^+ $$ So $q^\prime(x)=(I−P(x))^+P′(x)q(x)$ indeed satisfies the equation $$ (I-P(x)) q'(x)=P′(x)q(x) $$ Moreover there is no other solution $q'(x)$ because otherwise we will get $$ (I-P(x))(q'(x)-q_2'(x))=0 $$ namely $q_2'$ would have a component along $q(x)$.