Hello MathOverflow community,

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

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

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 probability matrix $\boldsymbol{P}$ having stationary distribution $\boldsymbol{\pi}$ (column vector) 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}^\mathrm{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 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 inputing these in a "Newton's method-like" program, but all it outputed 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 "probability 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 1 (same dimensions as $\boldsymbol{P}$).

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.