[Fix] decrease tolerance on forecast

[skoudoro]

The aim of this PR is to fix #1654 by decreasing the solver tolerance. At the moment, the current tolerance is 0.001, which I think is too high for getting any quality solutions.

• Advantage: The solution is more reliable and stable, and hopefully it will fix the negative variables that violate our positivity constraints
• Drawback: We have a tiny loss of speed.

As you can see here, cvxpy (and osqp) improved its solver. In fact, some unnecessary extra variables were introduced during the problem resolution. Below, the previous version verbose result:

with cvxpy v1.0.9
-----------------------------------------------------------------
           OSQP v0.4.1  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2018
-----------------------------------------------------------------
problem:  variables n = 249, constraints m = 272
          nnz(P) + nnz(A) = 7247
settings: linear system solver = qdldl,
          eps_abs = 1.0e-03, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: on

iter   objective    pri res    dua res    rho        time
   1   0.0000e+00   1.00e+00   1.11e+04   1.00e-01   1.61e-03s
 175   1.4997e-01   2.58e-04   9.77e-04   5.50e-01   7.11e-03s

status:               solved
solution polish:      unsuccessful
number of iterations: 175
optimal objective:    0.1500
run time:             7.85e-03s
optimal rho estimate: 1.53e-01

here, you can see that variables n = 249 is really strange because our input matrix, in this test, has a shape of (193, 28). It will make more sense to have 221 variables. This is what we get with the new cvxpy release as you can see below.

with cvxpy v1.0.10
-----------------------------------------------------------------
           OSQP v0.4.1  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2018
-----------------------------------------------------------------
problem:  variables n = 221, constraints m = 244
          nnz(P) + nnz(A) = 7191
settings: linear system solver = qdldl,
          eps_abs = 1.0e-03, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: on

iter   objective    pri res    dua res    rho        time
   1   0.0000e+00   1.00e+00   1.11e+04   1.00e-01   1.27e-03s
 200   1.3808e-01   4.50e-03   2.42e-05   1.00e-01   6.75e-03s
 250   1.4711e-01   1.49e-03   9.48e-04   5.71e-01   8.38e-03s

status:               solved
solution polish:      unsuccessful
number of iterations: 250
optimal objective:    0.1471
run time:             9.35e-03s
optimal rho estimate: 3.20e-01

To me, this indicates that the 2 versions of cvxpy are solving a different problem. Fixing the tolerance issue may be a good starting point for our test to stop failing. Does it make sense @maurozucchelli @arokem @Garyfallidis ?

[arokem]

Why not use the following?

Suggested change

	for _ in range(-l, l + 1):
	for l in range(0, sh_order + 1, 2):
	diag_lb[l//2] = (l * (l + 1)) ** 2

And get rid of the counter variable? Wouldn't that have the same effect?

[skoudoro]

I do not think so, because in your case, we will not fill diag_lb, it misses some iteration or another trick

[arokem]

Yes. This looks like a good solution to me, and fixes the test failure on my machine.

[codecov-io]

Codecov Report

❗ No coverage uploaded for pull request base (master@c2a3150). Click here to learn what that means.
The diff coverage is 100%.

@@            Coverage Diff            @@
##             master    #1655   +/-   ##
=========================================
  Coverage          ?   87.31%           
=========================================
  Files             ?      246           
  Lines             ?    32613           
  Branches          ?     3552           
=========================================
  Hits              ?    28477           
  Misses            ?     3275           
  Partials          ?      861

Impacted Files	Coverage Δ
dipy/reconst/forecast.py	92.22% <100%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update c2a3150...f887b77. Read the comment docs.

[arokem]

Nice!

[arokem]

I am +1 for merging this. Might be good for others to also take a look.

[Garyfallidis]

Thanks @skoudoro. All PRs should rebase.