I've written a function that returns the steady state given a regular, column-stochastic matrix. I want to use it to solve a larger problem, finding steady states that satisfy certain conditions.

The function seems to work fine on its own. But not as part of the larger problem.

Here's the function definition:

SteadyState[M_] :=
First[Eigenvectors[M]] / Total[First[Eigenvectors[M]]]


Used on its own it seems to work well, e.g.:

SteadyState[
{{1/2, 1/2},
{1/2, 1/2}}
]


returns {1/2, 1/2} as it should.

But when I try to use it inside of Solve, it doesn't behave as expected. For example:

Solve[
{{c11, c12},
{c21, c22}}
] == {q1, q2} &&
q1 >= 0 && q2 >= 0 &&
q1 + q2 == 1 &&
c11 >= 0 && c12 >= 0 && c21 >= 0 && c22 >= 0 &&
c11 + c21 == 1 &&
c12 + c22 == 1 &&
c11 == 1/2 && c12 == 1/2 && q1 == 1/2 && q2 == 1/2,
{q1, q2, c11, c12, c21, c22},
Reals
]


Returns the empty set of solutions {}, even though I've effectively handed it the solution in the last two constraints!

Could this be because my SteadyState function doesn't check that the input is regular and column stochastic? If so, why? And what would be a good way to rectify that?

• I suspect your system is under-determined, since you have 5 equations but 6 unknowns. Feb 22, 2019 at 17:08
• Also, if your problem is larger than 2x2, you might need to do this numerically. Feb 22, 2019 at 17:08
• This seems overly complicated. Which variables are known and which do you want to solve for? Feb 22, 2019 at 17:12
• What do you mean by "minimizes own expected error"? Feb 22, 2019 at 18:52
• Solve[SteadyState[{{c11, c12}, {1 - c11, 1 - c12}}] == {q1, 1 - q1}, {c11, c12}] gives an answer. Feb 22, 2019 at 18:57

I should have seen right away, for this purpose it's sufficient—and way more sensible—to just put the steady-state condition into Solve by hand, as one of the constraints.

To simplify the problem's expression I first defined some helpers:

ColumnStochastic[M_] := AllTrue[
Table[i, {i, Length[M]}],
Total[M[[All, #]]] == 1 &
] && AllTrue[Flatten[M], # >= 0 &]

pvec[v_] := AllTrue[v, # >= 0 &] && Total[v] == 1


The following code then runs no problem:

A = {{a11, a12}, {a21, a22}};
q = {q1, q2};

FindInstance[
A.q == q &&
ColumnStochastic[A] && pvec[q],
Join[q, Flatten[A]]
]


It finds the solution:

{{q1 -> 3/4, q2 -> 1/4, a11 -> 27/32, a12 -> 15/32, a21 -> 5/32, a22 -> 17/32}}


Using Solve or Reduce instead of FindInstance also works now to get the general class of solutions.

To find steady-states that meet additional constraints (as alluded in the OP), I can now just && them in alongside the others.