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.