Investigate failures of block decomposition when computing steady states

The file solve_v8.mod fails computing the terminal steady state with trust region, while it works with fsolve. Part of the problem is that dogleg returns NaN due to 0*Inf. parameter.mat steadystate.mat

This seems to be caused by blocks computed from the Dulmage-Mendelsohn decomposition having singular Jacobians.

Isn’t that to be expected that some solvers fail on some problems while other succeed, and precisely the reason why we provide several solvers? Or is there a good reason to think that it’s a bug in the implementation of trust region and not a normal outcome for that specific algorithm?

I am not familiar enough with the workings of the algorithm to judge that. But for two reasons I think it's worth investigating.

1. Having a linear combination with weight 0 returning NaN due to putting 0 weight on an infinite value strikes me as a case one may want to trap. For example

    if alpha>0
        x = alpha*x + (1.0-alpha)*min(sgnorm, delta)*s;
    else
        x = (1.0-alpha)*min(sgnorm, delta)*s;
    end

2. The issue seems to be not so much with trust_region.m. solve_algo=9 (trust region on the full model) works. It's solve_algo=4's

Splits the model into recursive blocks and solves each block in turn using a trustregion solver with autoscaling

that fails due to a singular Jacobian. So maybe there is a problem in the decomposition.

changed the description

@sebastien The mex-replacement solve_algo=13 (#1818) fails with a cryptic Non-square block. Is that intended? solve_v8.mod steadystate.mat parameter.mat

added bug label

assigned to @sebastien

changed title from Investigate pathological case in trust_region.m to Investigate failures of block decomposition when computing steady states

changed the description

The file noffrc.mod also triggers the same issue with the MEX, but also crashes solve_algo=4 due to a zero Jacobian in one of the blocks.

mentioned in merge request !2035 (merged)

I had a look at noffrc.mod. The issue is that the Jacobian of the model is singular. In that case, the Dulmage-Mendelssohn decomposition has a rectangular block, and I don’t know how to use such a rectangular block when solving the system.

This case is properly trapped in the MEX for solve_algo=13 (though the message is a bit cryptic and should probably be improved). In solve_algo=4, this case is not properly trapped, and the code does computations that do not make sense (by the way, is the commit in !2035 (merged) related to that? if yes, it should probably reverted, because the code should not go that far).

So I think we should just improve error messages to make it clear that the Jacobian is singular, and that the user should either fix that, or use another algorithm.

1. !2035 (merged) is about issues that may arise within iterations in a block. Here, we would need to filter out a pathological case before entering the solutions algorithms.
2. Yes, we should improve the error message in the mex. @sebastien Could you do that?
3. What to do with solve_algo=4 depends on our decision for #1818. If we replace it with a mex after allowing for a jacobian_flag=false, then nothing needs to be done here.
4. Do I understand it correctly that a singular Jacobian may be due to the starting value? That would make it necessary to point to different starting values or a different optimizer.

4. Yes of course singular Jacobian may be due to bad starting value. In theory it could also arise in the course of the algorithm, if it jumps to a point where the Jacobian is singular. And of course it could also come from a structural problem of the model, though that case is more easy to avoid.

unassigned @sebastien

Ok, so we should first settle #1818 and then revisit this topic.

I think we can already trap the singular Jacobian case in solve_algo=4. It’s just about testing the squareness of blocks, and possibly print the same error message as in the MEX.

assigned to @sebastien

closed with commit fa71d471

mentioned in commit 5ccd5f26

I’ve reverted commit fa71d471, because it turns out that it breaks very simple .mod files (such as tests/arima/mod1a.mod; see https://git.dynare.org/Dynare/dynare/-/jobs/41662 for more).

For those files, the Jacobian is singular at the start of the trust-region+Dulmage-Mendelsohn algorithm, and the algorithm does something which is nonsense (it treats a rectangular block as if it were square). But it somehow works (maybe because the initial guesses are the steady state).

We need to find a better approach…

reopened

mentioned in commit JohannesPfeifer/dynare@5788f1bc

mentioned in merge request !2049 (merged)

In !2049 (merged) I try to

1. first check whether the initial guess already solves the problem
2. in case it does not, condition on the type of blocks. Solving an underdetermined block based on a smaller square system should be fine, while doing the same for the overdetermined system is wrong.

mentioned in commit 5670d393

The MATLAB version of the trust region solver has been dealt with. The same now needs to be done in the block_trust_region MEX.

closed with commit 16eabbbc

mentioned in commit 4d2c8c78