# Of all probability matrix $P$ having stationary distribution $\pi$, find the one having smallest diagonal

Hello MathOverflow community,

I am requesting your help today trying to solve a somewhat odd problem. Is there a way to find through some numerical algorithm such as Newton's method the stochastic matrix $$\boldsymbol{P}$$ having stationary distribution $$\boldsymbol{\pi}$$ (column vector, given as an input) and lowest diagonal (whose elements are closest to 0)?

I initially thought about minimizing the following loss function:

$$L(\boldsymbol{P}) = {(\mathrm{diag}(\boldsymbol{P}))}^\mathsf{T} \mathrm{diag}(\boldsymbol{P}) + \lambda_1 {||\boldsymbol{P}^\mathsf{T}\boldsymbol{\pi} - \boldsymbol{\pi}||}^2 + \lambda_2 {||\boldsymbol{P}\boldsymbol{e} - \boldsymbol{e}||}^2,$$

where $$\mathrm{diag}$$ returns the diagonal of its argument, $$\lambda_1$$ and $$\lambda_2$$ are both Lagrange multipliers, and $$\boldsymbol{e}$$ is a vector filled with 1 (same dimensions as $$\boldsymbol{\pi}$$).

As far as I know, this is a case of differentiating a scalar ($$L$$) with respect to a matrix ($$\boldsymbol{P}$$). I think it may involve something called tensors, but to be honest I have little to zero experience with this, and even less once you throw in Lagrange multipliers.

I did some calculations, and it would appear the differential is given by

$$\tfrac{\partial{}L}{\partial{}\boldsymbol{P}} = 2\mathrm{\boldsymbol{D}}(\boldsymbol{P}) - 2\lambda_1\boldsymbol{\pi}\boldsymbol{\pi}^\mathsf{T}(\boldsymbol{I} - \boldsymbol{P}) - 2\lambda_2 \boldsymbol{e} \boldsymbol{e}^\mathsf{T} (\boldsymbol{I} - \boldsymbol{P}^\mathsf{T}),$$

which at some point involved an "outer product" or "Kronecker product", but could be simplified to that. The $$D$$ function outputs a diagonal matrix having same diagonal as its argument. In turn, the Hessian matrix (matrix of second order derivatives) would be given by

$$\tfrac{\partial{}^2L}{\partial{}\boldsymbol{P}\partial{}\boldsymbol{P}^\mathsf{T}} = 2\boldsymbol{I} + 2\lambda_1\boldsymbol{\pi}\boldsymbol{\pi}^\mathsf{T} + 2\lambda_2 \boldsymbol{e} \boldsymbol{e}^\mathsf{T}.$$

I tried inputting these in a "Newton's method-like" program, but all it outputted was gibberish.

All of this is a bit out of my league, but I really tried to make it work by myself before running here. I would be so grateful if someone could help me out. I know a solution exists, because Excel's solver is able to find solutions (don't ask why I use Excel, in this case I don't have a choice).

Thanks,

RSMax

P.S. Just in case there would be multiple definitions around, by "stochastic matrix" I mean a square matrix whose elements are probabilities, and whose rows all sum to 1.

P.P.S. By stationary distribution, I am referring to "long terms odds", as given by:

$$\boldsymbol{\pi} = {(\boldsymbol{I} + \boldsymbol{E} - \boldsymbol{P}^\mathsf{T})}^{-1} \boldsymbol{e},$$

where $$\boldsymbol{E}$$ is a square matrix filled with ones (same dimensions as $$\boldsymbol{P}$$).

EDIT: Fixed some typos and added some clarifications regarding what is an input.

• Usually this sort of matrices is called "stochastic". – Dieter Kadelka Jun 11 at 7:01
• @Dieter Kadelka You are perfectly right, in fact I initially learned about these in a class named "Stochastic Processes", so in makes so much sense, but since I'm translating everything in my head as I write, I simply didn't make the connection right away. – RSMax Jun 11 at 17:16

If you make the objective to minimize the sum of the diagonal entries (i.e. the trace), your problem becomes a linear programming problem, solvable with readily available software (I think even Excel). In many cases the optimal solution will have all diagonal entries $$0$$.

EDIT: It may help to think of the problem this way. There are $$n$$ people, numbered from $$1$$ to $$n$$. Each person $$i$$ has a nonnegative amount $$\pi_i$$ of money, and we want to move all the money around so that everyone ends up with the same amount they started with. $$P_{ij}$$ is the fraction of person $$i$$'s money that goes to person $$j$$. Of course this will not be possible if one person has more than half the total amount of money. Otherwise it will be possible, by the following algorithm. We arrange the people in order clockwise from richest to poorest around a circular table, and everyone puts their money on the table in front of them (let's say in coins of a given denomination). The richest persion (#1) takes the amount $$\pi_1$$ off the table, starting with #2's coins, which come after his own coins. Each person then continues, taking the correct number of coins starting where the previous one left off. By induction, for $$k$$ from $$2$$ to $$n$$, when #k's turn comes, all his coins have already been taken.

EDIT: Here is a more explicit version of the algorithm I alluded to above.
We may assume that $$\pi_1 \ge \pi_2 \ge \ldots \pi_n > 0$$. If necessary, we reorder the rows and columns to make them decrease. If some $$\pi_i$$ are $$0$$, we make the corresponding columns all $$0$$ and put a $$1$$ arbitrarily (off the diagonal) in each of the corresponding rows, and follow the following algorithm for the rows and columns where $$\pi > 0$$.

If $$\pi_1 > 1/2$$, we let $$P_{11} = 2 - 1/\pi_1$$, $$P_{i1}=1$$ and $$P_{1i} = \pi_i/\pi_1$$ for $$i > 1$$, and all other $$P_{ij} = 0$$. It is easy to check that this works, and no solution can have $$P_{11} < 2 - 1/\pi_1$$.

Now suppose $$\pi_1 \le 1/2$$. Let $$S_k = \sum_{i=1}^k \pi_i$$ be the partial [[sums of $$\pi$$, so $$S_0 = 0$$ and $$S_n = 1$$. Let $$r_{ij}$$ be the length of the intersection of the intervals $$[S_{i-1}, S_i]$$ and $$[S_{j-1}+\pi_1, S_j + \pi_1] \mod 1$$. Then $$P_{ij} = r_{ij}/\pi_i$$.

For example, suppose $$\pi = [0.3, 0.24, 0.24, 0.16, 0.06]$$. The partial sums are $$[0, 0.3, 0.54, 0.78, 0.94, 1]$$. The shifted partial sums $$S_i + \pi_1$$ are $$[.3, .0.6, 0.84, 1.08, 1.24, 1.30]$$. Then, for example, the intersection of $$[S_{2} + \pi_1, S_3 + \pi_1] \mod 1 = [.84, 1] \cup [0, .08]$$ and $$[S_0, S_1]=[0,0.3]$$ has length $$0.08$$, so $$P_{13} = 0.08/0.3 = 4/15$$, while the intersection of $$[S_2+\pi_1, S_3 + \pi_1] \mod 1$$ with $$[S_3, S_4] =[.78, .94]$$ has length $$.1$$ so $$P_{43} = .1/.16 = 5/8$$ and the intersection of $$[S_2+\pi_1, S_3 + \pi_1] \mod 1$$ with $$[S_4, S_5] = [.94, 1]$$ has length $$.06$$ so $$P_{53} = .06/.06 = 1$$. The resulting matrix is $$\pmatrix{ 0 & 0 & 4/15 & 8/15 & 1/5\cr 1 & 0 & 0 & 0 & 0\cr 1/4 & 3/4 & 0 & 0 & 0\cr 0 & 3/8 & 5/8 & 0 & 0\cr 0 & 0 & 1 & 0 & 0\cr}$$

• Alternatively, an objective to minimize $\sum_i P_{i,i}^2$ yields a (convex) quadratic programming problem. – RobPratt Jun 11 at 15:04
• I totally agree with the EDIT section, this is something I have observed in practice using Excel's solver. What remains encouraging is that through this scheme, $P_{11}$ will nonetheless be less than $\pi_1$, given that $\pi_1 \geq 0.5$, ergo the diagonal has still been minimized in some way. – RSMax Jun 11 at 17:03
• Whether I minimize the trace or the sum the squares of all elements on the diagonal, my current issue (under the Newton's method approach) is that my "main" unknown is a matrix ($\boldsymbol{P}$, whereas the Lagrange multipliers are scalars, and I'm confused about how to wrap those 3 under a common structure so that I can recursively update it and converge toward a solution. Should I add rows and columns to the unknown $\boldsymbol{P}$ matrix so that it incorporates $\lambda_1$ and $\lambda_2$ on the main diagonal? – RSMax Jun 11 at 17:08
• If by linear programming you mean something else than Newton's method, I am very much interested. As I need to ultimately code the result, anything easier and/or more efficient than Newton's method will be received with open arms. Would you mind elaborating a bit? Thank you sir. – RSMax Jun 11 at 17:14
• Linear programming in Excel. But the algorithm I presented in the EDIT is simpler (though maybe not so easy to program in Excel). Maybe I should make it more explicit. – Robert Israel Jun 11 at 19:49

The following is too large for a comment and definitely not an answer.

Fix $$n \in \mathbb{N}$$ and let $$\cal{P}$$ be the set of stochastic $$n \times n$$-matrices with $$\text{diag}(P) = 0$$. Let further $$\mathcal{E} := \left\{\pi \in \mathbb{R}^n \colon \pi_i \geq 0, \sum \pi_i = 1, \exists P \in \mathcal{P} \text{ with } \pi \cdot P = \pi\right\},$$ each $$\pi$$ being a row vector. Then both $$\mathcal{P}$$ and $$\mathcal{E}$$ are compact convex subsets of $$\mathbb{R}^{n \times n}$$ resp. $$\mathbb{R}^n$$ and not empty if $$n \geq 2$$ and for $$n = 2$$ we have $$\mathcal{E} = \{(0.5,0.5\}$$. The problem is that $$(1,0,\ldots,0), \ldots, (0,0,\ldots,1) \not\in \mathcal{E}$$.

But for large $$n$$ it seems (I have no proof) it seems that almost any random stochastic $$\pi$$ is in $$\mathcal{E}$$. For this case your original problem is equivalent to "Find a $$P \in \mathcal{P}$$ with $$\pi \cdot P = \pi$$". Here the explicit construction of a homogeneous Markov chain with transition probability matrix $$P$$ (appropriately chosen) may be helpful to get a "feeling" for the solution. Of course this is impossible to formalize.

It seems to me that for $$n$$ large using numerical methods may be too heavy for this sort of problem. Edit (I've just seen the answer of Robert Israel) A reformulation as a linear programm with criterion $$\sum P_{ii} = \min!$$ and variables $$P_{ij}$$ is possible.

• But $\pi_i$ is given, right? – RobPratt Jun 11 at 15:01
• Right, I just edited my answer. – Dieter Kadelka Jun 11 at 15:04
• I tried with a $10\times{}10$ matrix using Excel's solver and it took about half an hour. That's much too slow for my purpose. I doubt the order of the $\boldsymbol{P}$ matrix will ever exceed $n = 10$ states though. In fact, I won't allow it too. – RSMax Jun 11 at 17:11
• If you have access to a unix/linux computer you should try the glpk package on mirrors.ocf.berkeley.edu/gnu/glpk It's very much faster than excel. After installing this package you have glpsol. With glpsol you can solve among others linear programs. There are examples where you can start – Dieter Kadelka Jun 11 at 17:23
• @RSMax: If you have a windows computer, there is a windows version of glpk: winglpk.sourceforge.net There is glpsol.exe in the package. I don't know anything about the quality of this solver. It should be much faster than excel, but you have to invest some time. – Dieter Kadelka Jun 12 at 11:50