BUG: bvls can fail catastrophically to converge

[divenex]

The function bvls, which is called by optimize.lsq_linear(... method='bvls') can fail catastrophically to converge due to an incorrect implementation of the published BVLS algorithm.

It is supposed to implement the BVLS algorithm by Stak & Parker 1995, which is a small generalization of the classic NNLS algorithm 23.10 by Lawson & Hanson (1995).

However, the current Scipy implementation fails to implement the so-called loop-B in the notation of the above two papers.

The proposed changes fix the deficiency and allow the algorithm to converge in the cases where it currently fails.

Reference issue

Here is a minimal example where the current implementation fails

seed = 159
np.random.seed(seed)

n, m = np.sort(np.random.randint(2, 4, size=2))
A = 1.*np.random.randint(-99, 99, size=[m, n])
b = 1.*np.random.randint(-99, 99, size=[m])
bounds = 1.*np.sort(np.random.randint(-99, 99, size=(2, n)), 0)

print("seed", seed, "dims:", m, n, "rank: ", np.linalg.matrix_rank(A))
print(A)
print(b)
print(bounds)
print("")

x2 = optimize.lsq_linear(A, b, bounds=bounds, method='bvls', verbose=2).x
x3 = optimize.lsq_linear(A, b, bounds=bounds, method='trf', verbose=2).x

print(f"   x_bvls = {x2}")
print(f"    x_trf = {x3}\n")

chi2_2 = chi2(A @ x2 - b)
chi2_3 = chi2(A @ x3 - b)
chi2min = min(chi2_2, chi2_3)

print(f"rel_err   bvls: {chi2_2/chi2min - 1:.2e}")
print(f"rel_err    trf: {chi2_3/chi2min - 1:.2e}")

And here is the corresponding output of Scipy 1.5

seed 159 dims: 3 3 rank:  3

[[ 49.  41. -32.]
 [-19. -32.  -8.]
 [-13.  10.  69.]]

[-41. -90.  47.]

[[ 31. -44.  26.]
 [ 54. -32.  28.]]

Iteration        Cost      Cost reduction    Step norm     Optimality   
       0         7.4506e+05                                    4.80e+04    
       1         3.5282e+05      3.92e+05       1.64e+01       1.68e+04    
       2         2.4232e+05      1.11e+05       1.03e+01       1.14e+04    
The relative change of the cost function is less than `tol`.
Number of iterations 3, initial cost 7.4506e+05, final cost 2.4232e+05, first-order optimality 1.14e+04.

   Iteration        Cost      Cost reduction    Step norm     Optimality   
       0         8.0576e+05                                    1.05e+06    
       1         3.9417e+05      4.12e+05       8.77e+00       2.59e+05    
       2         2.8801e+05      1.06e+05       5.67e+00       6.47e+04    
       3         2.3206e+05      5.60e+04       5.84e+00       1.08e+04    
       4         2.1981e+05      1.22e+04       1.60e+00       6.26e+02    
       5         2.1920e+05      6.11e+02       2.46e-01       1.73e+02    
       6         2.1915e+05      4.84e+01       1.21e-01       4.13e+01    
       7         2.1914e+05      7.86e+00       5.52e-02       8.58e+00    
       8         2.1914e+05      7.29e-01       1.95e-02       1.15e+00    
       9         2.1914e+05      1.83e-02       3.49e-03       4.08e-02    
      10         2.1914e+05      2.47e-05       1.33e-04       4.93e-05    
      11         2.1914e+05      2.43e-10       1.61e-07       1.35e-10    
The relative change of the cost function is less than `tol`.
Number of iterations 12, initial cost 8.0576e+05, final cost 2.1914e+05, first-order optimality 1.35e-10.

   x_bvls = [ 54.         -39.82549256  26.        ]
    x_trf = [ 54.         -43.89055258  26.        ]

rel_err   bvls: 1.06e-01
rel_err    trf: 0.00e+00

Additional information

The current implementation fails in various other examples, with more realistic larger arrays, and especially when the matrix A is rank deficient. However, I selected a very simple example to make it easy to understand what the code is doing differently in the proposed changes.

[tylerjereddy]

Would it be practical to add a regression test based on the minimal example provided?

[nmayorov]

Is rcond parameter significant for this algorithm? If so we should pass it through lsq_linear.

Can we get away with default rcond=-1 which uses machine epsilon to cutoff small singular values? I wouldn't introduce another parameter if there are no cases for it.

[divenex]

The usefulness of rcond in bvls is the same as for linalg.lstsq. In most cases the default is enough, but the extra generality can be useful in some instances, as the default test for rank deficiency is too restrictive e.g. when the matrix contain measurement errors in addition to just numerical errors.

[nmayorov]

I think this algorithm was implemented based only on the BVLS paper. They don't go into details on implementation and only mention "Loop B" briefly. As far as I understand the problem arises because of finite precision computations, but of course it is important for numerical methods.

Would it be practical to add a regression test based on the minimal example provided?

I think at least one test to trigger this problem is a must. And, if possible, I would suggest extending the tests for this algorithm significantly as it proved to be fragile. Adding a test with "with more realistic larger arrays" would be very suitable.

[divenex]

Although "loop B" is briefly mentioned, the algorithm in Sec. 2 of Stark & Parker (1995), which closely follows Lawson & Hanson, states that one should go to step 2 (exit loop B) if the lstsq solution in step 6 is feasible, but go back to step 6 (remain in loop B) if it is unfeasible and one had to compute alpha to interpolate the solution onto the feasible region.

Indeed dealing with numerical errors are critical here. The typical problem is that one may select the wrong variable to be added to the free set and, when using a test on the relative cost decrease, the algorithm may incorrectly conclude that no decrease is possible and terminate.

The Lawson & Hanson and Stark & Parker algorithms avoid having a simple test of convergence based on the relative cost decrease between iterations, to avoid premature convergence. Instead, they rely on the KKT test only. However, this requires using slightly more complex techniques to avoid cycling of the active set (discussed in the papers), which can cause endless loops. In my experience, the current fixed implementation seems robust enough.

NOTE: The need for loop B is because, just after the linear interpolation of the solution, in general, the solution will not be a minimum of the quadratic function for the given active set. For this reason, the test of convergence in loop A will be invalid: the procedure should never stop just after interpolating the solution. Iterating within Loop B makes sure that the solution is a minimum of the quadratic function for the given active set before exiting the loop and testing for convergence in loop A.

[divenex]

Here is a more realistic example with large rank-deficient matrices

seed = 0
np.random.seed(seed)
n, m = np.sort(np.random.randint(2, 1000, size=2))
m *= 2
A = 1.*np.random.randint(-99, 99, size=[m, n])
b = 1.*np.random.randint(-99, 99, size=[m])
bounds = 1.*np.sort(np.random.randint(-99, 99, size=(2, n)), 0)
bounds[1, :] += 1  # make up > lb

# Make the A matrix strongly rank deficient by replicating some columns
w = np.random.choice(n, n)  # Select random columns with duplicates
A = A[:, w]

which gives before the fix

seed 0 dims: 1372 561 rank:  357

rel_err   bvls: 8.67e-01
rel_err    trf: 0.00e+00

and after the fix

rel_err   bvls: 0.00e+00
rel_err    trf: 3.68e-13

[nmayorov]

@divenex it all sounds very good. Will you add 2 tests on which current implementation fails?

As for rcond parameter -- add it to lsq_linear signature and pass to all lstsq calls.

[divenex]

@nmayorov where/how would I add the tests?

[nmayorov]

This is the file with tests https://github.com/scipy/scipy/blob/master/scipy/optimize/tests/test_lsq_linear.py

Add tests into BaseMixin. If you have any issues with it, I'll be happy to hep.

[nmayorov]

Any problems with adding tests from #12816 as another commit here?

[rgommers]

Any problems with adding tests from #12816 as another commit here?

I blindly cherry-picked those tests, they had some syntax issues. Fixed those and made the linter happy too, but the tests are failing - this needs work it looks like.

[nmayorov]

@rgommers

It was a trivial syntax error. Now tests are active and passing.

I think we should merge it finally.

[rgommers]

Great, glad to see this merged. Thanks @divenex and @nmayorov