# Solving a linear system of symbolic variables and symbolic parameters

I am trying to run the following code to solve a system of 4 equations.

    Reduce[ftft FT + ptft PT + uft U + nft OLF == FT &&
ftpt FT + ptpt PT + upt U + npt OLF == PT &&
ftu FT + ptu PT + uu U + nu OLF ==
U && (1 - ftft - ftpt - ftu) FT + (1 - ptft - ptpt -
ptu) PT + (1 - uft - upt - uu)  U + (1 - nft - npt -
nu) OLF == OLF, {FT, PT, U, OLF}]


However, when I run this, it takes an incredibly long time (I left it to run overnight and did not complete by the morning). Do you have any advice as to what the problem is? I'm new to Mathematica and it seems like this is a relatively simple system of equations.

Any help is appreciated.

• Because this is a system of homogeneous linear equations, the four variables can only be zero, unless the determinant is zero. – bbgodfrey Dec 8 '15 at 2:52
• Hmm. Perhaps I should explain a little more. I am essentially trying to get a closed form solution for the steady state of a discrete markov chain and then hoping to take the closed form to the data. The parameters (e.g. ftft, ftpt) are all nonzero and strictly less than 1. So, 0 is a solution but there should be a non-zero solution as well. – A. Rambachan Dec 8 '15 at 2:59
• (1) It is slow in part due to a phenomenon known as "coefficient swell" wherein symbolic terms get quite large. Consider what the result would look like if computed with Cramer's rule. You will have 17 determinants involving 4! sums of products of 4 matrix and vector elements. – Daniel Lichtblau Dec 8 '15 at 15:10
• (2) Moreover Reduce will look for nongeneric solutions wherein special relations between parameters change the form of solution set. This may be what you want but it is also going to require serious computation. – Daniel Lichtblau Dec 8 '15 at 15:12
• (3) On the topic of what you want, I would suspect it is not a symbolic solution per se. My guess is you mostly want a recipe for computing a solution given specific parameter values. For that I would advocate a numeric approach. This has the advantage of not suffering from truncation or cancellation error, each of which can make a mess out of a numeric substitution into a symbolic solution. Also numeric formulations can be solved quite fast. – Daniel Lichtblau Dec 8 '15 at 15:14

It sounds like what you are looking for is the eigenvalue corresponding to the eigenvalue 1. Why not approach this directly?

amat = {{a, b, c, d}, {e, f, g, h}, {i, j, k, l},
{(1 - a - e - i), (1 - b - f - j), (1 - c - g - k), (1 - d - h - l)}};
Eigenvalues[amat]
Eigenvectors[amat]


As you can see, there is an eigenvalue at 1, and the corresponding eigenvector is about 3 lines of symbols.

• [Why didn't I think of that?] And an upvote of course. – Daniel Lichtblau Dec 8 '15 at 17:38

I do not know why Reduce is so slow, but Solve works well, recognizing that the determinant of the coefficients is zero. (Edit: result below simplified)

Solve[{ftft FT + ptft PT + uft U + nft OLF == FT,
ftpt FT + ptpt PT + upt U + npt OLF == PT,
ftu FT + ptu PT + uu U + nu OLF == U,
(1 - ftft - ftpt - ftu) FT + (1 - ptft - ptpt - ptu) PT + (1 - uft -
upt - uu) U + (1 - nft - npt - nu) OLF == OLF},
{FT, PT, U, OLF}] // Simplify

(* {{PT -> (FT (-(ftu nft + nu - ftft nu) upt + npt (-1 + ftft + ftu uft + uu - ftft uu) -
ftpt (nft + nu uft - nft uu)))/(nu ((-1 + ptpt) uft - ptft upt) -
npt (ptft + ptu uft - ptft uu) + nft (-1 + ptpt + ptu upt + uu - ptpt uu)),
U -> (FT (-nu (-1 + ftft + ftpt ptft + ptpt - ftft ptpt) + ftu (nft + npt ptft -
nft ptpt) + (ftpt nft + npt - ftft npt) ptu))/(nu (uft - ptpt uft + ptft upt) +
npt (ptft + ptu uft - ptft uu) - nft (-1 + ptpt + ptu upt + uu - ptpt uu)),
OLF -> -((FT (-1 + ptpt + ftu uft - ftu ptpt uft + ftu ptft upt + ptu upt + uu -
ptpt uu + ftpt (ptft + ptu uft - ptft uu) - ftft (-1 + ptpt + ptu upt + uu -
ptpt uu)))/(nu (uft - ptpt uft + ptft upt) + npt (ptft + ptu uft - ptft uu) -
nft (-1 + ptpt + ptu upt + uu - ptpt uu)))}}*)


which expresses {PT, U, OLF} in terms of FT. Of course, one could equally well express any three of the variables in terms of the fourth.

Incidentally, the determinant of the system of equations can be shown to be zero by

{b, m} = CoefficientArrays[{ftft FT + ptft PT + uft U + nft OLF == FT,
ftpt FT + ptpt PT + upt U + npt OLF == PT,
ftu FT + ptu PT + uu U + nu OLF == U,
(1 - ftft - ftpt - ftu) FT + (1 - ptft - ptpt - ptu) PT + (1 -
uft - upt - uu) U + (1 - nft - npt - nu) OLF == OLF},
{FT, PT, U, OLF}];
Det[m // Normal]
(* 0 *)


b is, of course, identically zero.

• @A.Rambachan By the way, Solve finds general solutions only, while Reduce tries to find special solutions too. Perhaps, it tried too hard. – bbgodfrey Dec 8 '15 at 3:28
• The det need not be zero. The point is that 1 is an eigenvalue for the matrix. To see this, note that {1,1,1,1} is a left eigenvector for that eigenvalue. The goal in the post is to find a right eigenvector. – Daniel Lichtblau Dec 8 '15 at 15:18
• I had in mind the original version of the matrix: {{a, b, c, d},{e, f, g, h}, {i, j, k, l}, {(1 - a - e - i), (1 - b - f - j), (1 - c - g - k), (1 - d - h - l)}}. I gather you incorporated the -b part so that 0 replaces 1 as the important eigenvalue. (Why did you change the naming of parameters? Was there something stated in comments I have missed?) – Daniel Lichtblau Dec 8 '15 at 15:56
• Ah. Your answer and my edit crossed paths. This makes things more clear. Or at least less murky. Whichever better describes the state of affairs. – Daniel Lichtblau Dec 8 '15 at 16:00