4
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ How do you call FindRoot? $\endgroup$ – user21 Apr 23 at 12:06
  • 1
    $\begingroup$ You're missing assignments for the variables {t1, \[Rho]} in the parameters & initial values. Something's missing from your code? $\endgroup$ – Michael E2 Apr 23 at 12:26
  • $\begingroup$ @Michael E2: No, I do not think I miss something. All the parameters you refer to are in the list of parameters. $\endgroup$ – Tugrul Temel Apr 23 at 12:45
  • 3
    $\begingroup$ @TugrulTemel, you should add this the post and not in a comment where people have to search for it. Please remember to post complete working examples. Having incomplete code in posts demotivates looking at and helping you. $\endgroup$ – user21 Apr 23 at 12:51
  • 2
    $\begingroup$ No, you are mistaken: i.stack.imgur.com/CSnwN.png -- I think there is (at least) one typo is the posted code. $\endgroup$ – Michael E2 Apr 23 at 12:52
7
$\begingroup$

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.

$\endgroup$
7
$\begingroup$

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}
$\endgroup$
  • $\begingroup$ A nice way to announce the new method. :) $\endgroup$ – Michael E2 Apr 23 at 13:20
  • 1
    $\begingroup$ If you had not been so darn quick .... ;-) $\endgroup$ – user21 Apr 23 at 13:21
  • 1
    $\begingroup$ @MichaelE2 and thanks for finding that typo....! $\endgroup$ – user21 Apr 23 at 13:21
  • 1
    $\begingroup$ Well, that is an impediment for using Version12.0 functionality ;-) $\endgroup$ – user21 Apr 23 at 13:54
  • 1
    $\begingroup$ @TugrulTemel I fixed the typo for you and added the solution. $\endgroup$ – user21 Apr 23 at 15:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.