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

Question

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

Answer 1

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} $$

Answer 2

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.