`FindRoot [ ]`::jsing: Encountered a singular Jacobian at a point...WHY

Question

I have a large system of linear and nonlinear equations. I know that the system is not singular at the point given because I calculated the Jacobian and evaluated it at a given point. The resulting matrix is invertible.

jacob=D[system,{vars}]; (*for the Jacobian matrix*)
Inverse[jacob]; (*inverse of the Jacobian*)
Det[Jacob]; (*determinant of Jacobian*)

Here is the message I get.

FindRoot::jsing: Encountered a singular Jacobian at the point`  {EG,EXR,FSAV,IADJ,PA1,PA2,PD1,PD2,PE1,PM2,PQ1,PQ2,PVA1,PVA2,PX1,PX2,QA1,QA2,QD1,QD2,QE1,QF11,QF12,QF21,QF22,QH11,QH12,QH21,QH22,QINT11,QINT12,QINT21,QINT22,QINV1,QINV2,QM2,QQ1,QQ2,QX1,QX2,WALRAS,WF1,YF11,YF12,YF21,YF22,YG,YH1,YH2} = {109.,1.,4.,1.,1.,1.,1.,1.,1.,1.,1.,1.,0.6,0.5,1.,1.,279.,394.,249.,394.,30.,72.,105.,73.,135.,28.,47.,157.,87.,84.,55.,50.,99.,28.,85.,105.,259.,538.,279.,394.,0.,1.,95.,125.,82.,83.,109.,285.,186.}. Try perturbing the initial point(s).`

To me, I miss something important to get a result from FindRoot[] (maybe a working precision level is too high). In case one wants to replicate the calculations, I provide the required data below:

system={PM2 - EXR pwm2 (1 + tm2), PE1 - EXR pwe1 (1 - te1), 
PQ1 QQ1 - PD1 QD1 (1 + tq1), 
PQ2 QQ2 - (PD2 QD2 + PM2 QM2) (1 + tq2), -PD1 QD1 - PE1 QE1 + 
PX1 QX1, -PD2 QD2 + PX2 QX2, PA1 - PX1 \[Theta]11 - PX2 \[Theta]12, 
PA2 - PX1 \[Theta]21 - PX2 \[Theta]22, -PA1 + ica11 PQ1 + ica21 PQ2 +
PVA1, -PA2 + ica12 PQ1 + ica22 PQ2 + PVA2, 
QA1 - ad1 QF11^\[Alpha]11 QF21^\[Alpha]21, 
QA2 - ad2 QF12^\[Alpha]12 QF22^\[Alpha]22, 
WF1 WFDIST11 - (PVA1 QA1 \[Alpha]11)/QF11, 
WF1 WFDIST12 - (PVA2 QA2 \[Alpha]12)/QF12, 
WF2 WFDIST21 - (PVA1 QA1 \[Alpha]21)/QF21, 
WF2 WFDIST22 - (PVA2 QA2 \[Alpha]22)/QF22, -ica11 QA1 + 
QINT11, -ica12 QA2 + QINT12, -ica21 QA1 + QINT21, -ica22 QA2 + 
QINT22, QX1 - QA1 \[Theta]11 - QA2 \[Theta]21, 
QX2 - QA1 \[Theta]12 - QA2 \[Theta]22, 
QQ2 - aq2 (QD2^-\[Rho]q2 (1 - \[Delta]q2) + 
QM2^-\[Rho]q2 \[Delta]q2)^(-1/\[Rho]q2), 
QM2/QD2 - ((PD2 \[Delta]q2)/(PM2 (1-\[Delta]q2)))^(1/(1 + \[Rho]q2)),
-QD1  + QQ1, 
QX1 - at1 (QD1^\[Rho]t1 (1 - \[Delta]t1) + QE1^\[Rho]t1 \[Delta]t1)^(1/  
\[Rho]t1), 
QE1/QD1 - ((PE1 (1 - \[Delta]t1))/(PD1 \[Delta]t1))^(
1/(-1 + \[Rho]t1)), -QD2 + 
QX2, -shry11 (QF11 WF1 WFDIST11 + QF12 WF1 WFDIST12) + 
YF11, -shry12 (QF21 WF2 WFDIST21 + QF22 WF2 WFDIST22) + 
YF12, -shry21 (QF11 WF1 WFDIST11 + QF12 WF1 WFDIST12) + 
YF21, -shry22 (QF21 WF2 WFDIST21 + QF22 WF2 WFDIST22) + 
YF22, -tr13 - EXR tr14 - YF11 - YF12 + YH1, -tr23 - EXR tr24 - 
YF21 - YF22 + YH2, QH11 - ((1 - mps1) (1 - ty1) YH1 \[Beta]11)/PQ1, 
QH12 - ((1 - mps2) (1 - ty2) YH2 \[Beta]12)/PQ1, 
QH21 - ((1 - mps1) (1 - ty1) YH1 \[Beta]21)/PQ2, 
QH22 - ((1 - mps2) (1 - ty2) YH2 \[Beta]22)/PQ2, -IADJ qinv1 + 
QINV1, -IADJ qinv2 + QINV2, -EXR pwe1 QE1 te1 - EXR pwm2 QM2 tm2 - 
PD1 QD1 tq1 - (PD2 QD2 + PM2 QM2) tq2 - EXR tr34 + YG - ty1 YH1 - 
ty2 YH2, EG - PQ1 qg1 - PQ2 qg2 - tr13 - tr23, QF11 + QF12 - QFS1, 
QF21 + QF22 - QFS2, -qg1 - QH11 - QH12 - QINT11 - QINT12 - QINV1 + 
QQ1, -qg2 - QH21 - QH22 - QINT21 - QINT22 - QINV2 + QQ2, 
FSAV + pwe1 QE1 - pwm2 QM2 + tr14 + tr24 + tr34, -EG + EXR FSAV - 
PQ1 QINV1 - PQ2 QINV2 - WALRAS + YG + mps1 (1 - ty1) YH1 + 
mps2 (1 - ty2) YH2, -cpi + cwts1 PQ1 + cwts2 PQ2};

vars = {EG, EXR, FSAV, IADJ, PA1, PA2, PD1, PD2, PE1, PM2, PQ1, PQ2, 
PVA1, PVA2, PX1, PX2, QA1, QA2, QD1, QD2, QE1, QF11, QF12, QF21, 
QF22, QH11, QH12, QH21, QH22, QINT11, QINT12, QINT21, QINT22, QINV1,
QINV2, QM2, QQ1, QQ2, QX1, QX2, WALRAS, WF1, YF11, YF12, YF21, 
YF22, YG, YH1, YH2};

parameters={\[Rho]t1 -> 0.96, \[Rho]t2 -> 1.1, \[Rho]q1 -> 0.43, 
\[Rho]q2 ->   0.43, \[Alpha]11 -> 0.496552, \[Alpha]21 -> 
0.503448, \[Alpha]12 -> 0.4375, \[Alpha]22 -> 0.5625, \[Beta]11 -> 
0.153846, \[Beta]21 -> 0.846154, \[Beta]12 -> 0.347518, \[Beta]22 ->
0.652482, \[Theta]11 -> 1., \[Theta]12 -> 0., \[Theta]21 -> 
0., \[Theta]22 -> 1., \[Delta]t1 -> 0.742333, \[Delta]q2 -> 
0.340471, shry11 -> 0.536723, shry21 -> 0.463277, 
shry12 -> 0.600962, shry22 -> 0.399038, ad1 -> 3.84818, 
ad2 -> 3.25772, at1 -> 2.56505, aq2 -> 2.42675, ica11 -> 0.301075, 
ica21 -> 0.179211, ica12 -> 0.139594, ica22 -> 0.251269, 
cwts1 -> 0.235119, cwts2 -> 0.764881, mps1 -> 0.264151, 
mps2 -> 0.220994, tq1 -> 0.0401606, tq2 -> 0.0371747, 
ty1 -> 0.0701754, ty2 -> 0.0268817, te1 -> 0, tm2 -> 0.371429, 
qg1 -> 13, qg2 -> 67, qinv1 -> 28, qinv2 -> 85, tr13 -> 25, 
tr23 -> 5, tr14 -> 40, tr24 -> 16, tr34 -> 15, QFS1 -> 177, 
QFS2 -> 208, WFDIST11 -> 0.711, WFDIST12 -> 1.579, WF2 -> 1, 
WFDIST21 -> 1, WFDIST22 -> 1, cpi -> 1.036, pwe1 -> 1, pwm2 -> 1};

givenPoint= {
{EG, 109}, {EXR, 1}, {FSAV, 4}, {IADJ, 1}, {PA1, 1}, {PA2, 
1}, {PD1, 1}, {PD2, 1}, {PE1, 1}, {PM2, 1}, {PQ1, 1}, {PQ2, 
1}, {PVA1, 0.6}, {PVA2, 0.5}, {PX1, 1}, {PX2, 1}, {QA1, 
279}, {QA2, 394}, {QD1, 249}, {QD2, 394}, {QE1, 30}, {QF11, 
72}, {QF12, 105}, {QF21, 73}, {QF22, 135}, {QH11, 28}, {QH12, 
47}, {QH21, 157}, {QH22, 87}, {QINT11, 84}, {QINT12, 55}, {QINT21,
 50}, {QINT22, 99}, {QINV1, 28}, {QINV2, 85}, {QM2, 105}, {QQ1, 
259}, {QQ2, 538}, {QX1, 279}, {QX2, 394}, {WALRAS, 0}, {WF1, 
1}, {YF11, 95}, {YF12, 125}, {YF21, 82}, {YF22, 83}, {YG, 
109}, {YH1, 285}, {YH2, 186}
};

**

• EDIT 1

** For those who want to replicate the calculations:

FindRoot[system/.parameters, givenPoint]

should be enough.

Answer 1

Fixing the typo 1/\[Rho] t1 --> 1/\[Rho]t1, I can reproduce the error. I think it is just that the Jacobian is too poorly conditioned at the initial points:

LinearSolve[
  jacob /. parameters /. Rule @@@ initialVal]["ConditionNumber"]

LinearSolve::luc: Result for LinearSolve of badly conditioned matrix {{0,-1.37143,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},<<48>>} may contain significant numerical errors.

(*  4.75451*10^16  <-- condition number > 10.^MachinePrecision *)

Use a working precision of roughly around $\log_{10}(\hbox{condition number})$ plus the PrecisionGoal:

sol = FindRoot[SetPrecision[system /. parameters, 24], 
   SetPrecision[givenPoint, 24], PrecisionGoal -> 8, 
   WorkingPrecision -> 24] // N
(*
{EG -> 116.757, EXR -> 1.57378, FSAV -> 10.052, IADJ -> 1.51676, 
 PA1 -> 0.796624, PA2 -> 1.09585, PD1 -> 0.505239, PD2 -> 1.09585, 
 PE1 -> 1.57378, PM2 -> 2.15832, PQ1 -> 0.52553, PQ2 -> 1.19291, 
 PVA1 -> 0.424617, PVA2 -> 0.722743, PX1 -> 0.796624, PX2 -> 1.09585, 
 QA1 -> 330.749, QA2 -> 337.713, QD1 -> 361.553, QD2 -> 337.713, 
 QE1 -> 51.3489, QF11 -> 104.764, QF12 -> 72.2356, QF21 -> 70.705, 
 QF22 -> 137.295, QH11 -> 61.6308, QH12 -> 97.7302, QH21 -> 149.331, 
 QH22 -> 80.8366, QINT11 -> 99.5802, QINT12 -> 47.1428, 
 QINT21 -> 59.2738, QINT22 -> 84.8569, QINV1 -> 42.4692, 
 QINV2 -> 128.924, QM2 -> 132.401, QQ1 -> 361.553, QQ2 -> 570.222, 
 QX1 -> 330.749, QX2 -> 337.713, WALRAS -> -5.72072*10^-14, 
 WF1 -> 0.936218, YF11 -> 94.7432, YF12 -> 125., YF21 -> 81.7784, 
 YF22 -> 82.9999, YG -> 159.551, YH1 -> 307.694, YH2 -> 194.959}
*)

The solution is accurate at machine precision:

system /. parameters /. sol
(*
{0., 0., 0., 0., 0., 0., 0., 0., 0., 2.22045*10^-16, -5.68434*10^-14, 5.68434*10^-14, 
 1.11022*10^-16, 0., 0., -2.22045*10^-16, 0., 0., 0., 
-1.42109*10^-14, 0., 0., 1.13687*10^-13, 1.11022*10^-16, 0., -5.68434*10^-14, 
 2.77556*10^-17, 0., -1.42109*10^-14, 0., -1.42109*10^-14, 0., 0., 
 1.42109*10^-14, 1.42109*10^-14, 0., 0., 0., 0., 0., 0., 0., 0., -2.84217*10^-14, 0., 
 5.68434*10^-14, 0., 0., 0.}
 *)

Numerical analysis remark: Why is the error message "singular" and not "badly conditioned"?

By Gastinel's Theorem (which may be found here), a matrix differs from a singular one by a relative error of $1/K$, where $K$ is the condition number (with the size of the error and condition number measured with respect to the same operator norm). When that relative error is less than the error to be expected from ordinary rounding, then the matrix might have been obtained by rounding a singular matrix to the working precision. I don't know what tolerance FindRoot uses to determine whether a matrix is effectively singular, but given that in this case $1/K$ is less than $MachineEpsilon, which is the smallest expected nonzero relative rounding error in a single entry, it is not surprising that the Jacobian is considered singular. Artificially increasing the precision of the data in parameters and initialVal does not really change this fact. But since FindRoot is an iterative method, a bad first step, whether or not there was one, could be corrected in subsequent iterations. Further, FindRoot uses a "damped" method that shortens extremely large steps.

Answer 2

When you fix the typo you can use the new in 12.0 affine covariant Newton solver:

sol = FindRoot[system /. parameters, givenPoint, 
   Method -> {"AffineCovariantNewton"}];
MinMax[system /. parameters /. sol]
{-6.821210263296962`*^-13, 2.50716114535976`*^-13}

More documentation is here.

sol
{EG -> 116.757, EXR -> 1.57378, FSAV -> 10.052, IADJ -> 1.51676, 
 PA1 -> 0.796624, PA2 -> 1.09585, PD1 -> 0.505239, PD2 -> 1.09585, 
 PE1 -> 1.57378, PM2 -> 2.15832, PQ1 -> 0.52553, PQ2 -> 1.19291, 
 PVA1 -> 0.424617, PVA2 -> 0.722743, PX1 -> 0.796624, PX2 -> 1.09585, 
 QA1 -> 330.749, QA2 -> 337.713, QD1 -> 361.553, QD2 -> 337.713, 
 QE1 -> 51.3489, QF11 -> 104.764, QF12 -> 72.2356, QF21 -> 70.705, 
 QF22 -> 137.295, QH11 -> 61.6308, QH12 -> 97.7302, QH21 -> 149.331, 
 QH22 -> 80.8366, QINT11 -> 99.5802, QINT12 -> 47.1428, 
 QINT21 -> 59.2738, QINT22 -> 84.8569, QINV1 -> 42.4692, 
 QINV2 -> 128.924, QM2 -> 132.401, QQ1 -> 361.553, QQ2 -> 570.222, 
 QX1 -> 330.749, QX2 -> 337.713, WALRAS -> -4.07116*10^-14, 
 WF1 -> 0.936218, YF11 -> 94.7432, YF12 -> 125., YF21 -> 81.7784, 
 YF22 -> 82.9999, YG -> 159.551, YH1 -> 307.694, YH2 -> 194.959}