# Findroot apparently claiming to find singular Jacobians at all initial points for a system of 19 equations with 19 unknowns

The following code symbolically computes a system of equations and then applies FindRoot to it. When I run code the message 'Encountered a singular Jacobian...' comes up, for all initial points that I have tried.

k = 0.5;
fA = Max[0, Min[1, (1/h)*(Log[(pAC - k)/(pA - k)] + pAC/(pAC - k) - pA*(pA - k))]];
fB = Max[0, Min[1, (1/h)*(Log[(pBC - k)/(pB - k)] + pBC/(pBC - k) - pB*(pB - k))]];
zAB = Max[0, Min[1, 1/2 + (1/(2*g))*(Log[(pB - k)/(pA - k)] + pB/(pB - k) - pA/(pA - k))]];
zC = Max[0, Min[1, 1/2 + (1/(2*g))*(Log[(pBC - k)/(pAC - k)] + pBC/(pBC - k) - pAC/(pAC - k))]];
rAA = fA*zAB - (1/2)*Max[0, zAB - zC]*Max[0, fA - fB];
rBB = fB*(1 - zAB) - (1/2)*Max[0, (1 - zAB) - (1 - zC)]*Max[0, fB - fA];
rAC = (1 - fA)*zC - (1/2)*Max[0, zC - zAB]*Max[0, (1 - fA) - (1 - fB)];
rBC = (1 - fB)*(1 - zC) - (1/2)*Max[0, (1 - zC) - (1 - zAB)]*Max[0, (1 - fB) - (1 - fA)];
disutilityAA = fA*zAB*(h*(fA/2) + g*(zAB/2)) - (1/2)*Max[0, zAB - zC]*Max[0, fA - fB]*(h*(fA - (1/3)*(fA - fB)) + g*(zAB - (1/3)*(zAB - zC)));
disutilityBB = fB*(1 - zAB)*(h*(fB/2) + g*((1 - zAB)/2)) - (1/2)*Max[0, (1 - zAB) - (1 - zC)]*Max[0, fB - fA]*
(h*(fB - (1/3)*(fB - fA)) + g*((1 - zAB) - (1/3)*((1 - zAB) - (1 - zC))));
disutilityAC = (1 - fA)*zC*(h*((1 - fA)/2) + g*(zC/2)) - (1/2)*Max[0, zC - zAB]*Max[0, (1 - fA) - (1 - fB)]*
(h*((1 - fA) - (1/3)*((1 - fA) - (1 - fB)) + g*(zC - (1/3)*(zC - zAB))));
disutilityBC = (1 - fB)*(1 - zC)*(h*((1 - fB)/2) + g*((1 - zC)/2)) - (1/2)*Max[0, (1 - zC) - (1 - zAB)]*(h*((1 - fB) - (1/3)*((1 - fB) - (1 - fA))) +
g*((1 - zC) - (1/3)*((1 - zC) - (1 - zAB))));
everythingintermsofplayersactions = {xA -> 1/(pA - k), xAC -> 1/(pAC - k), xB -> 1/(pB - k), xBC -> 1/(pBC - k)};
\[Pi]A = (pA/(1 + tA) - c)*xA*(mA + rAA*mC) + (pAC/(1 + tC) - c)*xAC*rAC*mC //. everythingintermsofplayersactions;
\[Pi]B = (pB/(1 + tB) - c)*xB*(mB + rBB*mC) + (pBC/(1 + tC) - c)*xBC*rBC*mC //. everythingintermsofplayersactions;
focA = D[\[Pi]A, pA] == 0;
focAC = D[\[Pi]A, pAC] == 0;
focB = D[\[Pi]B, pB] == 0;
focBC = D[\[Pi]B, pBC] == 0;
dependenceontaxes = {pA -> pA[tA, tB, tC], pB -> pB[tA, tB, tC]};
focsp = {focA, focAC, focB, focBC} /. dependenceontaxes;
t = {tA, tB, tC};
dfocspbydt = D[focsp, {t}];
\[Pi]Aasafunctionoft = \[Pi]A /. dependenceontaxes;
\[Pi]Basafunctionoft = \[Pi]B /. dependenceontaxes;
sA = \[Pi]A + (tA/(1 + tA))*pA*xA*(mA + mC*rAA) + mA*(xA*k - xA*pA) //. everythingintermsofplayersactions /. dependenceontaxes;
sB = \[Pi]B + (tB/(1 + tB))*pB*xB*(mB + mC*rBB) + mB*(yB*k - xB*pB) //. everythingintermsofplayersactions /. dependenceontaxes;
sC = (tC/(1 + tC))*mC*(rBC*xBC*pBC + rAC*xAC*pAC) + mC*(k*(rAA*xA + rBB*xB + rAC*xAC + rBC*xBC) - rAA*xA*pA - rBB*xB*pB - rAC*xAC*pAC - rBC*xBC*pBC) -
disutilityAA - disutilityAC - disutilityBB - disutilityBC //. everythingintermsofplayersactions /. dependenceontaxes;
focsA = D[sA, tA] == 0;
focsB = D[sB, tB] == 0;
focsC = D[sC, tC] == 0;
allconditions = Join[Flatten[focsp], Flatten[dfocspbydt], {focsA}, {focsB}, {focsC}];
allconditionswithoutfunctions =
allconditions /. pAC[tA, tB, tC] -> pAC /. pA[tA, tB, tC] -> pA /. pB[tA, tB, tC] -> pB /. pBC[tA, tB, tC] -> pBC /. Derivative[1, 0, 0][pA][tA, tB, tC] ->
dpAdtA /. Derivative[0, 1, 0][pA][tA, tB, tC] -> dpAdtB /. Derivative[0, 0, 1][pA][tA, tB, tC] -> dpAdtC /. Derivative[1, 0, 0][pAC][tA, tB, tC] ->
dpACdtA /. Derivative[0, 1, 0][pAC][tA, tB, tC] -> dpACdtB /. Derivative[0, 0, 1][pAC][tA, tB, tC] -> dpACdtC /.
Derivative[1, 0, 0][pB][tA, tB, tC] -> dpBdtA /. Derivative[0, 1, 0][pB][tA, tB, tC] -> dpBdtB /. Derivative[0, 0, 1][pB][tA, tB, tC] -> dpBdtC /.
Derivative[1, 0, 0][pBC][tA, tB, tC] -> dpBCdtA /. Derivative[0, 1, 0][pBC][tA, tB, tC] -> dpBCdtB /. Derivative[0, 0, 1][pBC][tA, tB, tC] -> dpBCdtC;
parameters = {\[Rho] -> 1/2, c -> 1, mA -> 1/3, mB -> 1/3, mC -> 1/3, g -> 0.15, h -> 0.2};
allconditionswithoutfunctionswithnumbers = allconditionswithoutfunctions /. parameters;
FindRoot[{allconditionswithoutfunctionswithnumbers}, {{tA, 0.43}, {tB, 0.44}, {tC, 0.88}, {pAC, 6.1}, {pA, 4.}, {pBC, 7.1}, {pB, 3.2}, {dpAdtA, 1.1}, {dpAdtB, 1.2},
{dpAdtC, 1.05}, {dpACdtA, 1.06}, {dpACdtB, 1.03}, {dpACdtC, 1.02}, {dpBdtA, 1.04}, {dpBdtB, 1.03}, {dpBdtC, 1.01}, {dpBCdtA, 1.02}, {dpBCdtB, 1.04},
{dpBCdtC, 1.03}}]


The system consists of 19 equations and has 19 variables. Is it possible that this is due to a problem with FindRoot? (From the context of the problem it seems that there should be a solution and that the solutions should be locally unique.)

When I change the value of k to k=0 at the beginning, then I instead always get the error message 'The line search decreased the step size to within tolerance specified \ by AccuracyGoal and PrecisionGoal but was unable to find a sufficient \ decrease in the merit function. You may need more than \ MachinePrecision digits of working precision to meet these tolerances.'

• How does Sard's theorem collide with the fact that the mapping $f \colon \mathbb{R}^n \to \mathbb{R}^n$, $f(x)=0$ has noninvertible differential? Not at all. – Henrik Schumacher Mar 3 '18 at 14:27
• Sard's theorem only says that the set of all critical values is a null set. By definition, critical values are within the image of the mapping. – Henrik Schumacher Mar 3 '18 at 14:29
• Sorry, yes, I was mistaken. I have edited the question to correct it. – Lennart Mar 3 '18 at 14:41
• That does not suffice for Newton's algorithm (which is internally used in FindRoot). In patricular, I think the line search relies on the fact that the expression is comntinuously differentiable. There are algorithms that can deal with, e.g., semismooth expressions but I doubt that Mathematica employs such an algorithm; it just does not know how to compute a (generalized) derivative. However, you can help Mathematica by supplying a costumized Jacobian via the option Jacobian... – Henrik Schumacher Mar 3 '18 at 15:01
• The use of Min and Max means there may be "flat" regions, and this will manifest as regions where the Jacobian is rank deficient. – Daniel Lichtblau Mar 3 '18 at 16:23

In the meantime, I did some experiments. I repaced all occurences of Max[0,x] by Ramp[x] and Min[1,x] by 1 - Ramp[1 - x] since these can be easily differentiated by Mathematica.

Moreover, I computed the derivative of your expression

Off[Part::partd];
dpACdtB, dpACdtC, dpBdtA, dpBdtB, dpBdtC, dpBCdtA, dpBCdtB,
dpBCdtC};
XX = Table[X[[i]], {i, 1, Length[vars]}];
expr = allconditionswithoutfunctionswithnumbers[[All, 1]];
F = X \[Function] Evaluate[N[expr /. Thread[vars -> XX]]];
DF = X \[Function] Evaluate[N[D[F[XX], {XX, 1}]]];


Let's calculate the derivative at the initial point and check its rank:

X0 = {0.43, 0.44, 0.88, 6.1, 4., 7.1, 3.2, 1.1, 1.2, 1.05, 1.06, 1.03,
1.02, 1.04, 1.03, 1.01, 1.02, 1.04, 1.03};
Dimensions[DF[X0]]
MatrixRank[DF[X0]]


{19, 19}

12

Nope, it's not invertible, even with lots of good will.