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

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 No1: Fixed some typos and added some clarifications regarding what is an input.

EDIT No2: Substituting the sum of the squares of the diagonal elements for the trace indeed made this a linear problem which I was then able to solve using the simplex's "Big M" method. Works like a charm, even on Excel! My VBA code:

Public Const MAXVALUE As Double = 1.79769313486231E+308 'Arbitrary large value
Public Const EPSILON As Double = 0.0001 'Tolerance parameter for the Simplex method
Public Const BIGM As Double = 100# 'Big M constant for the Simplex method

Public Function Simplex(odds As Range) As Variant
'///////////////////////////////////////////////////////////////////////
'/// :function:    Simplex
'/// :scope:       Public
'/// :description: For odds.Cells.Count > 10, there can be instability issues.
'/// :return:      Wrapper (data validation) for the Simplex_ function.
'/// :odds:        An array of long-term odds (Range).
'///////////////////////////////////////////////////////////////////////
Dim n As Long
n = odds.Cells.Count

With WorksheetFunction
If .Count(odds) <> n Then
Simplex = CVErr(xlErrNum)
ElseIf .CountIf(odds, ">1") > 0 Then
Simplex = CVErr(xlErrValue)
ElseIf .CountIf(odds, "<0") > 0 Then
Simplex = CVErr(xlErrValue)
ElseIf Abs(.Sum(odds) - 1) > EPSILON Then
Simplex = CVErr(xlErrValue)
Else
Simplex = Simplex_(odds)
End If
End With
End Function

Private Function Simplex_(pi As Range) As Double()
'///////////////////////////////////////////////////////////////////////
'/// :function:    Simplex_
'/// :scope:       Public
'/// :description: Implements the "Big M" variant of the Simplex method.
'/// :return:      The lowest trace probability matrix having stationary distribution Odds.
'/// :pi:          An array of long-term odds.
'///////////////////////////////////////////////////////////////////////
Dim StopFlag As Boolean
Dim output() As Double
Dim MinRatio As Double
Dim MaxCost As Double
Dim factor As Double
Dim Ratio As Double
Dim Pivot As Double
Dim x() As Variant
Dim Row As Long
Dim Col As Long
Dim n As Long
Dim i As Long
Dim j As Long
Dim k As Long

n = pi.Cells.Count
'As we're looking for a matrix, the number of unknowns is n^2 + others for the constraints
ReDim x(0 To 2 * n + 1, 0 To n * n + 2 * n + 2) As Variant
ReDim output(1 To n, 1 To n) As Double
x(0, 0) = "BASIS"
x(0, 1) = "Z"

For i = 1 To n * n
Row = 1 + (i - 1) \ n
Col = 1 + (i - 1) Mod n
x(0, i + 1) = "p" & Row & Col
Next

For i = 1 To n
Col = 1 + (i - 1) Mod n
x(0, i + n * n + 1) = "c.r" & Col
x(0, i + n * n + n + 1) = "c.pi" & Col
x(i, 0) = "c.r" & Col
x(i + n, 0) = "c.pi" & Col
Next

x(0, n * n + 2 * n + 2) = "RHS" 'Right-hand side
x(2 * n + 1, 0) = "COST"

For i = 1 To n
For j = 1 To n * n
Row = 1 + (j - 1) \ n
x(i, j + 1) = IIf(Row = i, 1, 0)
Next

x(i, n * n + 2 * n + 2) = 1

For j = 1 To n * n
Row = 1 + (j - 1) \ n
Col = 1 + (j - 1) Mod n
x(i + n, j + 1) = IIf(Col = i, pi(Row).Value, 0)
Next

x(i + n, n * n + 2 * n + 2) = pi(i).Value
Next

For i = 1 To n * n
Row = 1 + (i - 1) \ n
Col = 1 + (i - 1) Mod n
x(2 * n + 1, i + 1) = IIf(Row = Col, -1, 0)
Next

For i = 1 To 2 * n
For j = i + 1 To 2 * n
x(i, n * n + 1 + j) = 0
x(j, n * n + 1 + i) = 0
Next

x(i, 1) = 0
x(i, n * n + 1 + i) = 1
x(2 * n + 1, n * n + 1 + i) = -BIGM
Next

x(2 * n + 1, 1) = 1
x(2 * n + 1, n * n + 2 * n + 2) = 0

For k = 1 To 2 * n
Col = n * n + 1 + k
MinRatio = MAXVALUE

For i = 1 To 2 * n
If x(i, Col) > 0 Then
Ratio = x(i, n * n + 2 * n + 2) / x(i, Col)

If Ratio < MinRatio Then
MinRatio = Ratio
Row = i
End If
End If
Next

Pivot = x(Row, Col)

For i = 1 To 2 * n + 1
If i <> Row Then
factor = -x(i, Col) / Pivot

For j = 1 To n * n + 2 * n + 2
x(i, j) = x(i, j) + factor * x(Row, j)
Next
End If
Next

For j = 1 To n * n + 2 * n + 2
x(Row, j) = x(Row, j) / Pivot
Next
Next

Do
MaxCost = -MAXVALUE

For i = 1 To n * n
If x(2 * n + 1, i + 1) > MaxCost Then
MaxCost = x(2 * n + 1, i + 1)
Col = i + 1
End If
Next

If MaxCost <= 0 Then
Exit Do
Else
MinRatio = MAXVALUE

For i = 1 To 2 * n
If x(i, Col) > 0 Then
Ratio = x(i, n * n + 2 * n + 2) / x(i, Col)

If Ratio < MinRatio Then
MinRatio = Ratio
Row = i
End If
End If
Next

Pivot = x(Row, Col)

For i = 1 To 2 * n + 1
If i <> Row Then
factor = -x(i, Col) / Pivot

For j = 1 To n * n + 2 * n + 2
x(i, j) = x(i, j) + factor * x(Row, j)
Next
End If
Next

For j = 1 To n * n + 2 * n + 2
x(Row, j) = x(Row, j) / Pivot
Next

x(Row, 0) = x(0, Col)
End If
Loop

For i = 1 To n * n
Row = 1 + (i - 1) \ n
Col = 1 + (i - 1) Mod n

For j = 1 To 2 * n
If x(j, 0) = "p" & Row & Col Then
output(Row, Col) = x(j, n * n + 2 * n + 2)
Exit For
End If
Next
Next

Simplex_ = output
End Function

• Usually this sort of matrices is called "stochastic". Jun 11, 2020 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. Jun 11, 2020 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. Jun 11, 2020 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. Jun 11, 2020 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? Jun 11, 2020 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. Jun 11, 2020 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. Jun 11, 2020 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? Jun 11, 2020 at 15:01
• Right, I just edited my answer. Jun 11, 2020 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. Jun 11, 2020 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 Jun 11, 2020 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. Jun 12, 2020 at 11:50