2
$\begingroup$

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
$\endgroup$
2
  • 2
    $\begingroup$ Usually this sort of matrices is called "stochastic". $\endgroup$ Jun 11, 2020 at 7:01
  • $\begingroup$ @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. $\endgroup$
    – RSMax
    Jun 11, 2020 at 17:16

2 Answers 2

3
$\begingroup$

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

$\endgroup$
5
  • $\begingroup$ Alternatively, an objective to minimize $\sum_i P_{i,i}^2$ yields a (convex) quadratic programming problem. $\endgroup$
    – RobPratt
    Jun 11, 2020 at 15:04
  • $\begingroup$ 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. $\endgroup$
    – RSMax
    Jun 11, 2020 at 17:03
  • $\begingroup$ 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? $\endgroup$
    – RSMax
    Jun 11, 2020 at 17:08
  • $\begingroup$ 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. $\endgroup$
    – RSMax
    Jun 11, 2020 at 17:14
  • $\begingroup$ 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. $\endgroup$ Jun 11, 2020 at 19:49
0
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ But $\pi_i$ is given, right? $\endgroup$
    – RobPratt
    Jun 11, 2020 at 15:01
  • $\begingroup$ Right, I just edited my answer. $\endgroup$ Jun 11, 2020 at 15:04
  • $\begingroup$ 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. $\endgroup$
    – RSMax
    Jun 11, 2020 at 17:11
  • $\begingroup$ 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 $\endgroup$ Jun 11, 2020 at 17:23
  • $\begingroup$ @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. $\endgroup$ Jun 12, 2020 at 11:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.