How to store all of the converged roots found by Matlab's fsolve algorithm?

Hi,

I'm practicing with Matlab's fsolve function, and am calling it with a few nested for loops that provide fsolve with various initial guesses.

How can I store all of the converged roots? Should I first store all of the roots, and then use another variable to store only the converged roots that satisfy the function tolerance?

For instance, I can write

% converged_roots = zeros(2, nguesses^2)      % Pre-allocate a 2 x nguesses^2 matrix of zeros to store the converged roots

nguesses = 20

for x_guess = linspace(-10, 10, nguesses)
    
    for y_guess = linspace(-5, 5, nguesses)
        
initial_guess = [ x_guess, y_guess ];

[ my_root, fval, exit_flag, output, Jacobian ] = fsolve( myfunction, initial_guess )

if norm( fval ) < 1e-5
    
converged_roots = my_root;

else
    
disp('not a converged root')

        end
    end
end

But, if I have, say, 5 converged roots found, the matrix 'converged_roots' only stores the last converged root -- so it's being overwritten each time in the loop.

How could I fix this?

If I pre-allocate a matrix of zeros (first line of my code, commented out), nothing happens either -- 'converged_roots' will still only return the last converged root. I think it might be because I'm not using a loop index i, j, anywhere in the code ...

Thanks,

Answer 1

[Walter Roberson]

Initialize

converged_roots = [];

and replace

converged_roots = my_root;

with

converged_roots{end+1} = my_root;

[Noob]

Hi Walter,

So, converged_roots = [ ] is an empty matrix stored in a variable?

What does {end + 1} do?

[Noob]

Also, I just tried this, and Matlab says consider pre-allocating for speed ...

[Noob]

It also doesn't work; I get this:

converged_roots =

1×100 cell array

Columns 1 through 8

[1×2 double] [1×2 double] [1×2 double] [1×2 double] [1×2 double] [1×2 double] [1×2 double] [1×2 double]

Columns 9 through 16

[1×2 double] [1×2 double] [1×2 double] [1×2 double] [1×2 double] [1×2 double] [1×2 double] [1×2 double]

Columns 17 through 24 ...

[Walter Roberson]

converged_roots{end+1} = my_root

means to create a new cell after the current end of converged_roots and store my_root in there.

It was a typing mistake on my part to use {end+1} instead of the (end+1) that I intended. However it turned out to be a fortunate mistake: without the mistake, we would not have discovered that your my_root is not a scalar value.

The fact that my_root is not a scalar value tells us that you are solving an equation of more than one variable, which is something you did not mention before.

The code is working. Each of your solutions is a vector, and the code is storing the vectors.

But you could modify the code to:

converged_roots(end+1,:) = my_root;

if you prefer to get out an array of roots.

Also, I just tried this, and Matlab says consider pre-allocating for speed

Any time you do not know in advance how many outputs you are going to get, you have two choices:

1. Expand the array by one slot every time you want to store another value. This will get you the caution to consider pre-allocating. It is not an error message, just a warning message that expanding arrays as you go is not the best for performance; OR
2. estimate the maximum size of the array and pre-allocate that. You are trying nguesses x nguesses combinations, and if every combination "converged", then you would have nguesses^2 different outputs. So you could pre-allocate converged_roots = zeros(nguesses^2, 2) and found_counter = 0; and then each time you found another, found_counter = found_counter + 1; converged_roots(found_counter,:) = my_root; . Then after all of the loops, trim away the parts not used: converged_roots(found_counter+1, :) = [];

[Noob]

Hi Walter,

I went to bed early last night -- thanks so much for this detailed comment.

In the last part of your comment, point no. 2, could I pre-allocate converged_roots = zeros(2, nguesses^2) instead, so that I store all of the converged roots as 2x1 column vectors, and so converged_roots has nguesses^2 columns to store the maximum possible number of converged roots? So for each converged root, I have converged_roots(: , found_counter) = my_root. Does this argument work?

And, why do we set converged_roots( : , found_counter+1 ) = [ ]? I understand the left-hand side, but why does the right-hand side consist of empty brackets?

Thanks,

[Noob]

Actually, the way you do it is much neater, in terms of spacing for printing / viewing solutions in the Command Window.

So, I just don't understand the last part: why you set the right-hand side to empty brackets [ ].

[Walter Roberson]

The technique allocates for the worst case where everything you try leads to a converged root. But you will usually not end up using all of the entries, so you will have an array in which the first count rows are used and the rest of the array is meaningless.

At that point you start wanting to do something with the results such as

disp(converged_roots)

but since only up to count is valid you would instead have to

disp(converged_roots(1:count,:))

This gets awkward, so it would be easier if we could get rid of the unused entries.

If we had initialized with nan instead of 0 then we could have used rmmissing(), but that function does more effort than we need, since we already know how much we want to save.

One way would be

converged_roots = converged_roots(1:count,:);

which explicitly extracts the part we want to keep and saves it. Some of the volunteers have demonstrated that this approach can have significant performance advantages, so you should certainly keep this in mind for the future. This goes along with a mental model of keeping a certain part of the array.

The other major approach is

converged_roots( count+1:end,:) = [];

Assigning literal [] to an element is a Mathworks syntax to indicate deletion. This goes along with the mental model of throwing away the parts you do not want.

There are some programming languages with dynamic memory in which deletion is more memory efficient because it can reuse the same memory space without moving it, returning any left-over memory to the memory pool. MATLAB does not work that way internally now but it could do that without notice (or just a mention of performance improvement). You should therefore not feel bad at all about writing the code in the way that matches your mental model.

[Noob]

Hi Walter,

Ok, thanks so much for your help - I really appreciate it.

[Walter Roberson]

By the way, you should expect that you might get a number of duplicates. Also, the algorithms used by fsolve do not guarantee that the root that will be found will be the "closest" to the original point.

You should consider using uniquetol() with 'byrows' to filter out the duplicates.

[Noob]

Hi Walter,

Good point about keeping the unique roots only -- thanks so much.

As for fsolve not guaranteeing that the root found will be "closest" to the minimizer, I am currently thinking that the main setting within fsolve to get the closest minimizer would be the 'First-Order Optimality' option. In the Mathworks documentation, it says that the lower the First-Order Optimality, the closer the function's gradient / Jacobian will be to zero -- analogous to finding 'critical points' in calculus, which may or may not be minimizers. So, I am thinking of setting the First-Order Optimality to a more stringent value, to get more accurate roots, e.g. if I set it to 1e-8, then, if I understand this option setting correctly, the root found by fsolve should be such that the infinity-norm of the gradient / Jacobian, evaluated at this root, will be less than 1e-8. What's not clear to me is the benefit of using the option 'FunctionTolerance', when we already have a function evaluation output argument from fsolve, when a root is found. I actually posted a separate question about this earlier -- trying to learn about fsolve's convergence analysis / metrics, and options that can be set. Please feel to comment / answer that posting, if you want.

Have a great day / night, and thanks again.

[Walter Roberson]

We are talking about different things for "closest". The factors you are talking about have to do with when the algorithm quits when it finds a point that "looks pretty good". What I was referring to is like this:

  S--
     \      +C
      \    /  \
  -----A--B----D-   0
        \/      \

You might be starting at S, and the closest zero might be at A. But the trial point might be at C, and since f(C ) < f(S) it would accept C -- and from there it might find the roots at B or D, neither of which are the closest to S.

With respect to a point "looking pretty good": because different functions are at different scales, f(x) = 1e-13 might be "better than can be expected considering round-off!" or could be "really quite bad" by comparison to how close it could get. Furthermore, f(x) = 1e-13 or even f(x) = 1e-50 might not be true roots and the true root might be quite a distance away. Consider for example (x^2 + 1e-50)*(x-20) which has a real root at 20 and at 0 has a value of -2e-49 which is not a root. The function tolerance always you to adjust to be more careful about roots -- or, on the flip side, to be less careful as there are cases where you only need the root to a couple of decimal places.

[Noob]

Hi Walter,

In your example (x^2 + 1e-50) * (x-20), we can find its roots by setting each factor equal to zero and solving, getting x_1 = 20 as one root, and the complex roots x_2 = -1e-50 * i , and x_3 = 1e-50 * i, which, as you said, are not zero. This is an illuminating example -- thanks so much.

Are the function evaluations, given as the second output argument by fsolve, sort of meaningless then? Is the Function Tolerance and the First-Order Optimality Tolerance the critical metrics to follow, when trying to find true roots of a function?

Thanks,

[Walter Roberson]

Remember that you are using a numeric solver with floating point numbers. x^2 + 1e-50 is indistinguishable from x^2 unless x^2 < 1e-34 and thus abs(x) < 1e-17 which is less than eps(1). The minimum step size is normally set larger than that. Numeric solvers have to assume that there is round-off error at the floating point level, so they often have an absolute tolerance below which they assume that a small enough absolute value is "really" just round off error on a calculation that yields 0.

[Noob]

Hi Walter,

I'm playing with your example in Matlab now -- so it appears that the choice of initial guesses is crucial; otherwise, it might appear that I have a root, when the solutions from fsolve are actually misleading. The first-order optimality and function tolerance doesn't seem to mitigate the issue either. I love this example. Thanks again.

[Noob]

Hi Walter,

I just saw your last comment, after I submitted an additional comment up there ^.

So if x^2 + 1e-50 is indistinguishable from x^2, then, numerically, x = 0 is considered a root of the function.

Is that correct to say?

And are these misleading "roots" still valuable to numerical analysts, say, from an optimization point of view?

Although they aren't true roots analytically, they would still give function evaluations that might satisfy the tolerances that are set.

Thanks,

[Walter Roberson]

x^2 +1e-50 is indistinguishable from x^2 for abs(x) > 1e-17 . But there are a lot of values between 0 and 1e-17 in double precision, so if fsolve happens to look in the area at all, it could potentially have some meaning.

The question is whether it will look in the area at all. For example at f = eps then f(x) is about 1e-31, so fsolve might have already decided that it had a good enough root. If the function tolerance were 1e-13 then the range +/- 0.00000007 satisfies that tolerance.

Thus, we should not be surprised if fsolve() decides that there is a root anywhere in +/- 0.00000007 unless we adjust the tolerances.

But no matter what absolute tolerance we use, fsolve() is going to have problems.

x^2 + sqrt(eps(realmin))

does not become distinguishable from x^2 until about 1e-72.

[Walter Roberson]

We could hypothesize that there might be different algorithms for looking at the tolerance of a root determined by jacobian (root has an even multiplicity so function does not cross 0 there), versus a root determined by looking for a zero crossing. But then we just substitute the function (x^2 - 1e-50) * (x-20) which does have 0 crossings (true roots) at +/- 1e-25 and ask ourselves whether we can realistically expect a solver that works by looking for zero crossings to be able to find such a root. Sure, a solver might deliberately try exactly 0, but we can obviously shift the root slightly ((x-epsilon)^2 - 1e-50)*(x-20) and have the same problem.