Chebop/linop convergence has gone awry

[asgeirbirkis]

It seems that recent changes in linop convergence decisions have destroyed a lot of the convergence behaviour for Newton iteration once we're close to the solution. This is definitely related to #898, and perhaps other pull requests. On e31a3c, which I think was a point I still was happy, running the following problem (Carrier eqn. with damping off):

dom = [-1 1];
N = chebop(@(x,u) 0.01.*diff(u,2)+2.*(1-x.^2).*u+u.^2-1, dom);
rhs = 0;
N.bc = @(x,u) [u(-1); u(1)];
x = chebfun(@(x) x, dom);
N.init =  2.*(x.^2-1).*(1-2./(1+20.*x.^2));
options = cheboppref();
options.display = 'iter';
options.errTol = 1e-10;
options.damped = false;
u = solvebvp(N, rhs, options);

gave me the following output:

Iter.   || du ||   Contraction    stepsize   len(du)   len(u)
--------------------------------------------------------------
 01     2.76e+00          NaN      1.0000     125       153 
 02     1.59e+00     5.77e-01      1.0000     165       165 
 03     1.53e+00     9.65e-01      1.0000     199       199 
 04     1.29e+00     8.38e-01      1.0000     197       199 
 05     9.86e-01     7.67e-01      1.0000     213       213 
 06     3.75e-01     3.80e-01      1.0000     215       215 
 07     6.60e-02     1.76e-01      1.0000     219       219 
 08     1.79e-02     2.71e-01      1.0000     215       219 
 09     5.24e-04     2.93e-02      1.0000     223       223 
 10     1.25e-07     2.38e-04      1.0000     168       223 
 11     2.11e-12     1.69e-05      1.0000     269       269 
--------------------------------------------------------------
Newton's method converged in 11 iterations.

If I run the same on development, changing options.damped to options.damping, I get

>> clear classes
Iter.   || du ||   Contraction    stepsize   len(du)   len(u)
--------------------------------------------------------------
 01     2.76e+00          NaN      1.0000     256       256 
 02     1.59e+00     5.77e-01      1.0000     256       256 
 03     1.53e+00     9.65e-01      1.0000     256       256 
 04     1.29e+00     8.38e-01      1.0000     256       256 
 05     9.86e-01     7.67e-01      1.0000     256       256 
 06     3.75e-01     3.80e-01      1.0000     256       256 
 07     6.60e-02     1.76e-01      1.0000     256       256 
 08     1.79e-02     2.71e-01      1.0000     256       256 
 09     5.24e-04     2.93e-02      1.0000     128       256 
 10     4.38e-07     8.35e-04      1.0000      32       256 
 11     4.97e-08     1.14e-01      1.0000      32       256 
 12     2.21e-08     4.44e-01      1.0000      32       256 
 13     1.69e-08     7.67e-01      1.0000      32       256 
 14     8.06e-09     4.77e-01      1.0000      32       256 
 15     3.04e-09     3.77e-01      1.0000      32       256 
 16     9.47e-10     3.12e-01      1.0000      32       256 
 17     2.35e-10     2.48e-01      1.0000      32       256 
 18     3.75e-11     1.60e-01      1.0000      32       256 
--------------------------------------------------------------
Newton's method converged in 18 iterations.

Similarly, solving the van der Pol demo from the GUI:

dom = [0 15];
N = chebop(@(x,u) diff(u,2)-(1-u.^2).*diff(u)+u, dom);
rhs = 0;
N.bc = @(x,u) [u(0)-2; feval(diff(u),0)];
x = chebfun(@(x) x, dom);
N.init =  2.*cos(x);
options = cheboppref();
options.display = 'iter';
options.errTol = 1e-10;
u = solvebvp(N, rhs, options);

gives the following back when I was happy:

Iter.   || du ||   Contraction    stepsize   len(du)   len(u)
--------------------------------------------------------------
 01     6.67e+00     3.61e-01     4.49e-01    151       151 
 02     1.92e+00     4.18e-01      1.0000     251       251 
 03     6.64e-01     3.46e-01      1.0000     318       318 
 04     1.52e-01     2.29e-01      1.0000     324       324 
 05     1.09e-02     7.18e-02      1.0000     328       328 
 06     6.36e-05     5.83e-03      1.0000     339       339 
 07     2.24e-09     3.52e-05      1.0000     434       434 
 08     4.36e-15     1.95e-06      1.0000     546       546 
--------------------------------------------------------------
Newton's method converged in 8 iterations.

But now I get:

Iter.   || du ||   Contraction    stepsize   len(du)   len(u)
--------------------------------------------------------------
 01     6.67e+00     3.61e-01     4.49e-01    256       256 
 02     1.92e+00     4.18e-01      1.0000     256       256 
 03     6.64e-01     3.46e-01      1.0000     256       256 
 04     1.52e-01     2.29e-01      1.0000     256       256 
 05     1.09e-02     7.18e-02      1.0000     256       256 
 06     6.36e-05     5.83e-03      1.0000     256       256 
 07     1.86e-07     2.92e-03      1.0000      32       256 
 08     3.34e-07          NaN      1.0000      32       256 
Warning: Newton iteration failed.

[trefethen]

Yes indeed. I have been slow in achieving convergence on guide10.m because of this.

[nickhale]

Both of these examples work for me in development in I switch to either plateauCheck() or classicCheck().

This suggests to me that V4Check() is too loose. Indeed, looking at the plotcoeffs() of the newton updates, that does seem to be the case.

[tobydriscoll]

I'm not saying we're wrong to revert...but there were complaints of overly long solutions when we were using plateauCheck. I'm concerned about running around in circles.

EDIT: Logically, either I did not really reproduce V4 functionality in the Check, or something else is responsible for variation from V4 experience. And linopV4Check is a pretty simple function. The only intentional difference is for system problems.

[tobydriscoll]

Tonight, when I use the devel branch, the Carrier example above is already fine. If I switch to plateauCheck, both of these are fine.

But the orrsomm test fails. We can fix that with a tolerance increase. However, the reason for the failure is instructive. Because it's more skeptical, plateauCheck demands larger values of N. These cause the lower coefficients to get a lot more polluted (by poor conditioning, I guess).

So aggressive truncation can be more accurate, not less. And if ultraS is as fast as collocation, we might make it the default.

[nickhale]

plateauCheck() at least demands some kind of plateau (hopefully?). My experiments in looking at plotcoeffs images of the Newton updates suggested that linopV4Check() was definitely giving up too soon.

In response to your last comment, I agree that in some cases aggressive truncation can be more accurate. However, I'm not sure that's the case most of the time. Even if I'm wrong, it's certainly true that aggressive truncation can be less accurate, too.

I don't think we're going to converge on this in the next 12 hours (pun intended), but I think the philosophy of plateauCheck() is the right one although it needs improving (particularly in the epslevel estimate).

One final note: I think it's more important that these iterations converge than that we're put off by having powers of two in the solution.

[asgeirbirkis]

This is the output I get from running the undamped Carrier problem in v4, using the same code for controlling the output:

   Iter.     || update ||    Contra.fact.      lambda       length(update)    length(u)
----------------------------------------------------------------------------------------
   1         2.759e+00              NaN       1.000e+00          179          179
   2         1.591e+00        5.767e-01       1.000e+00          207          207
   3         1.535e+00        9.646e-01       1.000e+00          255          255
   4         1.286e+00        8.376e-01       1.000e+00          253          255
   5         9.859e-01        7.669e-01       1.000e+00          271          271
   6         3.747e-01        3.800e-01       1.000e+00          269          271
   7         6.597e-02        1.761e-01       1.000e+00          283          283
   8         1.790e-02        2.713e-01       1.000e+00          269          283
   9         5.243e-04        2.929e-02       1.000e+00          203          283
  10         1.246e-07        2.377e-04       1.000e+00          159          283
  11         1.202e-12        9.640e-06       1.000e+00          107          283

The norms of the updated look identical to the output I posted for the results from e31a3c up until the 10th iteration, and then there is ever so slight change in the 11th iteration, which at least is quite comforting. The lengths of the update evolve in a much better behaved manner than the output I posted for development, where the jump really quickly from 256 -> 128 -> 32. So there appears to be a clear mismatch going on between v4 and the v4-style convergence results currently v5.

[tobydriscoll]

In the near future we should move more of these challenging problems into the test suite. (Preferably with true error checking and not residual checking.)

I think the plateauCheck idea can be tweaked, but you're always making decisions that are likely to be contradicted by another example. We don't (for speed reasons?) have a lot of tests of highly oscillatory problems, where you could get a plateau followed by real convergence. Capturing that while also giving up appropriately for a noise plateau is a tricky balance.

More fundamental thinking is needed, as we know.

[asgeirbirkis]

I'm confused. Even with #969 merged, I still get the following output for the undamped Carrier problem I posted above:

>> u = solvebvp(N, rhs, options);
Iter.   || du ||   Contraction    stepsize   len(du)   len(u)
--------------------------------------------------------------
 01     2.76e+00          NaN      1.0000     256       256 
 02     1.59e+00     5.77e-01      1.0000     256       256 
 03     1.53e+00     9.65e-01      1.0000     256       256 
 04     1.29e+00     8.38e-01      1.0000     256       256 
 05     9.86e-01     7.67e-01      1.0000     256       256 
 06     3.75e-01     3.80e-01      1.0000     256       256 
 07     6.60e-02     1.76e-01      1.0000     256       256 
 08     1.79e-02     2.71e-01      1.0000     256       256 
 09     5.24e-04     2.93e-02      1.0000     128       256 
 10     4.38e-07     8.35e-04      1.0000      32       256 
 11     4.97e-08     1.14e-01      1.0000      32       256 
 12     2.21e-08     4.44e-01      1.0000      32       256 
 13     1.69e-08     7.67e-01      1.0000      32       256 
 14     8.06e-09     4.77e-01      1.0000      32       256 
 15     3.04e-09     3.77e-01      1.0000      32       256 
 16     9.47e-10     3.12e-01      1.0000      32       256 
 17     2.35e-10     2.48e-01      1.0000      32       256 
 18     3.75e-11     1.60e-01      1.0000      32       256 
--------------------------------------------------------------
Newton's method converged in 18 iterations.
Discretization method used: Collocation. 

This behaviour keeps showing up all over the place.

[nickhale]

The problem always seems to occur just as there's a drop in the length of the update.

[nickhale]

This works fine on development now. Just awaiting a test.