No solution found using fsolve, what should I do?

I'm trying to solve a nonlinear system with six equations and six unknowns (Basically it is trying to find the fixed point of a mapping from to ).

I wrote the following small experiment:

tempx=1; 
p0=0.5*ones(1,6);
fsub=@(p)ComputeCCP(tempx,p);
[trueccp,fv,ext]=fsolve(fsub,p0);

It takes a long time to run this code, and in the end, it returns "No solution found". The function I feed into fsolve is called fsub, and it is a function of p. It is constructed using a function called ComputeCCP when its first parameter tempx is fixed at 1. I attached ComputeCCP in this question. It is slow to evaluate because I used 1 million simulations inside this function to compute double integrals (which gives probabilities) with high precision.

Why fsolve is not working here? What shall I do to find a solution to this system? Thanks in advance!

Answer 1

[John D'Errico]

A simple rule is fsolve assumes the function you feed it is well behaved. That means things like continuity. Differentiable. At the very least, it means that if you were to evaluate the function at the same set of values twice in a row, it would generate exactly the same result.

format long g
fsub(.5*ones(1,6))
ans =
  Columns 1 through 3
                  0.096333                  0.302968                   0.30198
  Columns 4 through 6
                   0.10278                  0.319573                  0.318503
                   
fsub(.5*ones(1,6))
ans =
  Columns 1 through 3
                  0.096683                   0.30191                  0.301975
  Columns 4 through 6
                  0.102587                  0.319008                  0.318855

Do you see anything interesting there? That when evaluated at exactly the same point twice in a row, your function is not even consistent, returning the same result. That means continuity and especially differentiability are completely out of the question.

If I look at your code:

z11=normrnd(0,3);
z12=normrnd(0,3);
z21=normrnd(0,3);
z22=normrnd(0,3);
   
e11=normrnd(0,3);
e12=normrnd(0,3);
e21=normrnd(0,3);
e22=normrnd(0,3);

it appears to be a simulation of some sort. Is that going to produce the well-defined function I said fsolve REQUIRES? No.

I'm sorry. It does not matter how badly you want to use fsolve, you cannot do so.

[Tmat]

Hi John, this is extremely helpful! Thank you so much!

Indeed, I used simulations inside the function to calculate double integral with irregular regions. What you said is very insightful. I guess at least one candidate solution is to generate these simulated values in advance and outside this function and set as parameters of the function ComputeCCP. This will make it producing the same value each time when it is called. I will try this.

What else could I do? Do you know any other solvers in Matlab or other methods that might fit my problem better? Thanks again!

[Torsten]

I don't know the background of your problem, but most probably, it is possible to solve it by Monte-Carlo simulations instead of solving deterministic nonlinear equations.

[Tmat]

@Torsten Thanks, Torsten! My problem is a deterministic problem. The Monte-Carlo simulation is used to approximate double integrals. I want to find the fixed point of a mapping whose calculation involves complicated double integrals on irregular regions. Monte-Carlo simulation is used to approximate this mapping.

[Torsten]

A Monte Carlo Simulation to approximate integrals within a root solver with Newton's method cannot work. The Monte Carlo Simulation can never be exact enough such that the root solver could approximate the derivatives of the functions which is necessary for the Newton step.

[Tmat]

@Torsten Thanks, Torsten! What other solvers or method do you suggest Also, what other method do you suggest to approximate the integrals?

[Torsten]

Could you explain your problem in more detail ? It's difficult to extract it from your code.

[John D'Errico]

I just saw your response. Yes, you COULD use one set of values for thiose parameters, then using fsolve. But that seems a poor idea. An important question is why you are using Monte Carlo to compute the integral.

d1(i,1)=1*(U1_00>=U1_10)*(U1_00>=U1_01)*(U1_00>=U1_11)+2*(U1_10>=U1_00)*(U1_10>=U1_01)*(U1_10>=U1_11)+3*(U1_01>=U1_00)*(U1_01>=U1_10)*(U1_01>=U1_11)+4*(U1_11>=U1_00)*(U1_11>=U1_10)*(U1_11>=U1_01);
d2(i,1)=1*(U2_00>=U2_10)*(U2_00>=U2_01)*(U2_00>=U2_11)+2*(U2_10>=U2_00)*(U2_10>=U2_01)*(U2_10>=U2_11)+3*(U2_01>=U2_00)*(U2_01>=U2_10)*(U2_01>=U2_11)+4*(U2_11>=U2_00)*(U2_11>=U2_10)*(U2_11>=U2_01);

That is just one part of your code, so the domain of interest seems rather complicated. But in fact, it does not seem really that nasty, and is the integration kernel even linear on that domain? I think it is, but reading your code is literally impossible. I don't really want to spend the time to unravel your code to learn exactly what you are doing, as there are no comments at all in it. It is not even clear to me how many dimensions the domain lives in.

So what is the domain of interest? Can tools like integral or integral2 be used to compute the integral instead? Have you tried?

[Tmat]

@Torsten Right. In my problem, I want to find the p such that p=fsub(p), where p is a probability vector of dimension 6. For any candidate p, the calculation of fsub(p) involves the calculation of double integrals, which I used Monte Carlo simulation to approximate.

[Tmat]

@John D'Errico Thanks, John. The integration kernals are highly nonlinear as they are the joint probability density function of of (e11+z11,e12+z12) and the joint probability density function of (e21+z21,e22+z22). In the current code I used a very simple case, that is eij and zij are both independent standard normal. But in the final most general code, eij and zij could be any continuous distributions that are arbitrarily correlated (ei1 and ei2 could also be correlated), and I want to be able to compute the fixed point, i.e., p=fsub(p) quickly without the need to find the analytical expression for the pdf of (ei1+zi1,ei2+zi2). Also, finally, I want to be able to easily scale this code to higher dimensions, for example, when the current U1_00 becomes U1_000 and others. Monte Carlo integral seems the only way that satisfy these needs because simulating from any joint distributions is easy. To do integral2, I need the analytical expression for the joint pdfs, right?

[Torsten]

@Tmat

Imagine you have an ignorant audience and you are to explain your mathematical problem. Do you think the audience would understand what you are talking about ? Or even to help you with your problem ?

[Tmat]

@Torsten Hi Torsten, sorry for the delayed response and thanks for your help. I solved my problem mentioned in this question by passing simulated values as parameters with extremely large simulation sample size. The solutions are also verified.