Could fsolve cause trouble when using it to solve linear equations?

Hello,

I have used fsolve optimisation algorithm to solve linear integral equations using integral command. But I know fsolve works for non-linear equations. So can I use it to solve linear equation system? Thanks

Answer 1

[Mohammad Abouali]

Yes, you can, but if you have linear system of equations you have much better and optimized function for system of linear equations. So, why do you want to use fsolve and not use those?

By the way, fsolve is not for optimization. Solving a system for zero is different than optimizing an objective function. In some cases they work just the same but generally one shouldn't be used for the other.

[Meva]

I have given below and I supposed this is optimisation metho and uses fsolve.

F = @(x) [f1x(x(1), x(2), x(3)); f2x(x(1), x(2), x(3)); f3x(x(1), x(2), x(3))];
  x0 = [1.5; 1.5; 1.5];
  options = optimoptions('fsolve','Display','iter'); % Option to display output
  [x,~] = fsolve(F,x0,options); % Call solver

[Mohammad Abouali]

Once you are using an optimization algorithm such as Nelder-mead, your goal is to find an input set that minimizes the output of the objective function.

This is very different than using algorithms such as fsolve that their goal is to find the solution of a system. In this case you are looking for an input set that causes the output of the function to be zero.

For example, take y=x^2-1; It has it's solution/roots at x=+1 and x=-1 but the minimum of the function happens at x=0 where the function evaluates to -1.

If you use fsolve you get this:

[x,fval]=fsolve(@(x) (x^2-1),5)

x =
    1.0000
fval =
   3.1836e-07

and for the other root

[x,fval]=fsolve(@(x) (x^2-1),-5)
x =
   -1.0000
fval =
   3.1836e-07

you got x=+1 and x=-1 and the function output (fval) was very close to zero. Those are the roots to that function.

Now if you use an optimization algorithm for the same function, i.e. y=x^2-1 you get:

[x,fval]=fminsearch(@(x) (x^2-1),5)

x =
     0
fval =
    -1

Now you get x=0 as the result and the function value (fval=-1). so fminsearch told you where the function is minimized/optimized, and fsolve told you where the roots of the functions are.

You can see now that if you use an optimization function or if you use fsolve on the same function you get completely different answer and thats ok.

So, once again, optimization and solving a system for its root are not the same thing. and FSOLVE IS NOT for optimizing a function but finding its roots (although you get default options of fsolve from optimset, but still that doesn't make fsolve an optimiztaion algorithm.

In some cases these two method can be used interchangably and that is only if the minimum of the function is zero In that case solving for its root and optimizing the function practically should give the same answer but those are special cases and one shouldn't mix these two.

So now back to your question again: if you are looking for the roots of a system of linear system, there are much better functions available in matlab for that, fsolve afterall is for nonlinear systems. You can use fsolve but udsually at a higher cost.

But if your goal is trully optimizatoin then fsovle is the wrong choice.

[Matt J]

and FSOLVE IS NOT for optimizing a function but finding its roots

Actually, fsolve offers various choices of algorithms and they are all designed to optimize a least squares cost function except for the trust-region-dogleg algorithm. I had a similar conversation about this with MathWorks reps here,

http://www.mathworks.com/matlabcentral/answers/57929-levenberg-marquardt-in-lsqnonlin-vs-fsolve

The dogleg algorithm is the default, but even if it is not a rigorous least squares solver for nonlinear problems, it looks strongly like a variation of Gauss-Newton ( see Documentation ). Therefore for linear problems, as we're discussing here, it should reduce to a least squares solver. The agreement between resnorm1 and resnorm2 in the example below seems to confirm this,

    A=[1,1,1; 2 1 3].'; b=[5 5 7].'; x0=rand(size(A,2),1);

    [x,fval]=fsolve(@(x) A*x-b, x0);

    resnorm1=norm(fval)
    resnorm2=norm(A*(A\b)-b)

[Mohammad Abouali]

I read the post by MathWorks Reps that you mentioned and he is saying the same thing that I am saying here.

You are missing the point again. You are mistaking the similirity between the algorithms and their goals. Finding root of a function is not the same as finding the minimum or optimum point of the function. Although they both use a similar approach and that's why the default options of the fsolve are obtained by optimset (I think). Any way, you can always pose a root finding problem as an optimization problem (any way, that's what root solving algorithms such as fsovle do anyway). In fact you root finding is a special sort of optimziation in which case you are defining the optimum point to be zero. but the reverse is not always true. the optimum point of a function is not always the point the function becomes zero.

As I gave an example above, fsolve for (x^2-1) gives the root of that function; not the optimum point of that function which the function value is the lowest. This is a simple example and we know that roots and optimum points are not the same. But in real problem and higher dimension, it might not be possible to even plot the function and know how it looks like. If someone learns fsolve as optimization, and pass a function without knowing that it is giving the root of the function; the answer could be very wrong if the function could have negative values. Therefore, it is very crucial that people know what each function is exactly doing and what is the purpose of it.

When you are optimizing function F(x1,x2,...) and you pass that to optimization functinos such as fminsearch, the objective function is F(...). (By the way I know you don't like fminsearch). However, when you use something like fsolve() and you pass F(x1,x2,...) that's not the objective function. One way you can think of it is that you form a new objective function G(x1,x2,...)=|F(x1,x2,...)-0| and now it is an optimization problem. That's how these two problem are related and that's how you could practically use any optimization algorithm to find the root of a function. But you can use the same method to find what inputs are required to have a certain output value such as G(...)=|F(...)-myValue|. In these cases (root finding included) G is being optimized not F.

In your example above norm(vec) is always bigger than or equal to zero. so in these two cases the optimum point and the root are the same point. But if you look at the X^2-1, then you could see the difference.

As long as your definition of th optimum point is when the function value is zero, then you can use fsolve and call it optimization because the roots of a function in that case are the optimum point. but not always. in those cases that the root is not the optimum point then you should be carefull. After all, all the algorithms in fsolve are not optimizing F but minimizing the devition of F from zero. so pretty much optimizing another function constructed from F.

But not always the roots of a function are their optimum points and the minimum point of a function most of the time is not the root of the function. So you should distinguish between these two.

[Matt J]

I read the post by MathWorks Reps that you mentioned and he is saying the same thing that I am saying here.

No, I don't think you're saying the same thing as the MathWorks reps. Unless I've misuderstood, you are saying that none of the algorithms used by fsolve are optimization theoretic, and that is not what Shashank claimed in that thread.

You and I do agree that root-finding and optimization are not the same thing and that the function F(x) that you feed to fsolve is not something that is being minimized directly. However, it is not just a hypothetical possibility that the root-finding problem can be converted to an optimization problem. That is in fact what fsolve is doing internally and it is using known optimization principles and algorithms to extract a solution from this formulation. After all, there has to be some reason why fsolve is part of the Optimization Toolbox!

The exception to this is fsolve's default trust-region dogleg algorithm. Shashank claimed that this algorithm is not motivated by an optimization problem. Accordingly, for nonlinear equations, there is no apparent guarantee that a zero-gradient point of norm(F(x))^2 will be found. However, I think that in the special case where F(x) is linear, it is equivalent to a least algorithm, and my example appears to confirm that as well.

[Mohammad Abouali]

I did not say that none of the internal algorithm within fsolve is not optimization algorithm. I said, "By the way, fsolve is not for optimization." or "FSOLVE IS NOT for optimizing a function but finding its roots". These are copy paste from my replies and comments. I am saying fsolve shouldn't be used to optimize a function.

On the contrary, I actually said in my reply to you that fsolve and root finding algorithms turn the problem into an optimization and they solve it. (second paragraph of my reply to you). But they turn it into a problem that makes sure that the optimum point is where the function output is zero, i.e. the root.

again, if you have function F() and you want to optimize it; passing it to fsolve(F,...) does not necessarily provide you with the correct answer and my example above also clearly shows that. The only time that it works is that if the roots of the function are the minimum point of the function, as in your example.

Well I think we are saying the same thing anyway. "You and I do agree that root-finding and optimization are not the same thing and that the function F(x) that you feed to fsolve is not something that is being minimized directly." That's all I am trying to say, if you are looking for the minimum/optimum point of F() don't pass it to fsolve, but if you are looking for the roots or you already know that the minimum/optimum is the zero of the function then fsolve can be used; although if the goal is solely finding the root of the system of linear equations; then in this case, you have other options too such as:

• cgs (Conjugate gradients squared method),
• GMRES(Generalized minimum residual method (with restarts)),
• bicgstab (biconjugate gradients stabilized method)
• and my favorite one, AGMG (AGgregation-based algebraic MultiGrid) available on http://homepages.ulb.ac.be/~ynotay/AGMG/

These are for finding the root or solving a system of linear equations. There exist many more of these algorithms specific for solving a system of linear equations, and again, many of these algorithms do minimize the sum of squares error again (GMRES even has the "minimum residual" as part of its name), but their goal is not optimization.

By the way, if you look at Meva's profile it says "fluids". If (s)he is dealing with a system of linear equations arising from solving Navier-Stokes equation then above algorithms are much better suited for the job than fsolve. Actually AGMG (or other implementation of multigrid) is way more popular (and efficient) in solving the system of linear equations resulted from Poisson's equation for pressure than any other method. Used with preconditioning, they can do magic.

Answer 2

[John D'Errico]

Think of it like this. Fsolve CAN solve a linear system of equations, but do you want to use the wrong tool for the job? In general, many tools can solve problems they are not really targeted to solve. But for a linear system, use a tool designed to solve that particular problem VERY well and efficiently.

A nonlinear solver will assume that you have a nonlinear problem. That means it will need starting values. It will need convergence tolerances. Fsolve will first need to compute the derivatives of what it thinks are nonlinear equations. (Or you will need to supply them.) So all of this makes fsolve less efficient for the problem than it need be.

After all, while I can probably drive in a nail by banging it with a screwdriver, a hammer will do so much better. And I suppose I could find a way to drive a screw with a hammer, but a screwdriver will be the right tool for the job.

So instead, just use backslash (or linsolve, or pinv, etc.) to solve linear problems. These tools are designed to be quite efficient and accurate on linear problems.

[Mohammad Abouali]

agree. People do all sort of tricks to keep the model to a linear system as long as it doesn't invalidate the results and it is possible. Here, we do have a linear system to begin with, why should we use fsolve which is for non-linear systems anyway!?

[gautam vatsa]

Hi I have a similar difficulty. I totally admit your points - why at all should we use fsolve for solving linear problems. The issue in my case is that the my linear equations arrive after discretizing my partial differential equations (which are although linear but very complicated - coupled and 4th order). So for me it is much easier to form a function handle than to create the matrix vector form. Can you suggest some linear equation solver that admits function handle (consisting of set of linear equations in variables) just like fsolve does.

The function equationToMatrix requires me to represent my variables in symbolic variables, however I have too many variables and therefore, using symbolic variables makes my codes run too slow. Please suggest something.

[Torsten]

If your equations are linear in the vector x of unknowns, use the function handle to calculate the Jacobian of your system (with Jacobianest, e.g.) (Result is A) and evaluate your function handle with all variables set to 0 (Result is b).

Then your linear system to solve reads

A*x = -b

such that

x = A\(-b).