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
```