# Response of stationary distribution to perturbation of a stochastic matrix

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.

• It's better to built such parametrization that guarantees stochastic properties of $P(x)$ for any $x$ in some domain. Commented Jul 19, 2014 at 3:34
• @AlexanderVigodner that's what I mean, I think. $P(x)$ is a differentiable function from (some interval on) $\mathbb{R}$ to the set of stochastic matrices. Commented Jul 19, 2014 at 4:51

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.

• $I-P(x)$ is not invertible by assumption, so this can't be quite right. Something like this is right, though.
– Ian
Commented Jul 19, 2014 at 3:14
• ^Good catch. I'll edit my solution. Commented Jul 19, 2014 at 3:15
• I follow the first 5 lines (+1) but I can't understand how the sentence "Hence the solutions ... will form a line in the direction of $q(x)$" can be true. (Surely those are the solutions to $(I-P(x))q'(x)=0$ rather than the equation stated?) There also seems to be quite a big jump from the penultimate to the last line. Maybe it will become more obvious if I read up on pseudoinverses, but could you write something to help my intuition along a bit? Commented Jul 19, 2014 at 5:58
• Jimmi, your first variant was correct, now you made mistake. The left eigenvector is not $q$ it is ALWAYS vector of ones. See my extending of your first variant. Commented Jul 19, 2014 at 7:12
• I'm using the convention that $P_{ij}$ is the probability of transitioning from state $j$ to state $i$. Thus, $q$ is the right eigenvector which represents the stationary distribution not the left eigenvector of all ones. (Else, $q$ would be constant w.r.t. $x$). However, I might still have made a mistake. Commented Jul 19, 2014 at 7:26

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)$.