Solve large number of independent systems with fsolve, where some systems do not have a solution

Suppose I have N nonlinear systems of equations for , where and are each vectors of length m. For each n, I want to feed in the vector and obtain the solution using fsolve.

I know that if the system has a solution for every , then I can do this by feeding in an -length vector of appended 's and solving for the vector of appended 's using a single call of fsolve, where I provide an appropriate structure of the Jacobian to specify that the N systems are independent.

I want to find the solution for each n where the solution exists. However, when there are even a few n where the solution does not exist, fsolve often will not find the solution even for those n where it does exist. I think this is because the few n without a solution throw off the sum of squares, so that fsolve "gives up" quickly on finding a solution for the entire -length system. Is there some way to "vectorize" fsolve so that it finds the exact solution for each n where the solution exists, or will I have to resort to looping over all n and using fsolve on each individual system?

Apologies if anything is unclear; this is my first time asking a question here.

Answer 1

[Matt J]

I see no reason why fsolve would "give up" if you've set MaxIterations and MaxFunctionEvaluations to Inf, regardless of whether a solution to all equations exists.

However, it is often not a good idea to try to combine smaller systems into an Nm-dimensional systems. For one thing some sub-systems could take a lot longer to converge than others, even if a solution to every subsytem does exist. With the Nm-dimensional formulation, you lose the ability to apply separate stopping criteria to each of the subsystems, and will be bottle-necked by the slowest system.

Also, although I don't know this to be true for your case, often the solution of one system can serve as a good initial x0 for the next equation system in the sequence, reducing the iterations that the next equation system will require. You lose the ability to exploit that if you try to solve all N systems simultaneously.

[Matthew Schwartzman]

Thanks, after examining further I realize that fsolve is not "giving up" early, but trying further iterations where it continues to change even for those n that have solutions, and actually increasing the error for these n by doing so.

In any case, I think you are right that I should probably just solve each m-dimensional system sequentially. I was taught at some point to avoid writing for-loops in Matlab if possible, but it looks like that is my best option.

[Matt J]

I just wanted to add that you could also consider a hybrid strategy where you divide the N systems into small batches of say 10, which you would solve as a 10*m-dimensional system. That way, the bottlenecking may not be as bad, and you would reduce the outer loop to N/10 iterations.