Stopping criteria confusing/Python code doesn't match Matlab

[jhtu]

I am trying to use SPGL1 to solve the Basis Pursuit Denoise problem:
minimize ||x||_1 subject to ||Ax - b||_2 <= sigma.

When I run spgl1 with the sigma keyword argument, it sometimes stops with exit status 1 (root found), but the sigma criterion is not met. I am confused why this is. It seems to depend on the opt_tol keyword argument, but the effect of this parameter is not clear from the BPDN equation nor from the spgl1 documentation.

I have attached a dataset and sample script to show the behavior I am talking about.

Furthermore, I noticed while preparing this script that the output of the Python and Matlab libraries is not identical. Is this to be expected?

spgl1_bug.zip

[mrava87]

Hi,
Thanks for reporting this. I will investigate and let you know what I find.

When you mean the output of Matlab and python is not identical, you mean the actual inverted x is different or the output on screen?

[jhtu]

Hi,

I am not at my work computer where I can run my scripts, but I recall that the residual rnorm and the number of iterations did not match. Since A, b, sigma, opt_tol, and the maximum number of iterations are all set the same, I would expect the algorithm to generate the same output.

[mrava87]

Hi,
I suggest trying to run the following modified code :

# Fix some parameters
sigma = 1e-4
iter_lim = 2000

# Run for different cases of opt_tol
for opt_tol in [1e-1,1e-2, 1e-4, 1e-7, 1e-10]:
    # Use spgl1 to solve BPDN problem
    x, r, g, info = spgl1.spgl1(
        A, b, sigma=sigma, opt_tol=opt_tol, iter_lim=iter_lim, verbosity=1)

    # Run some checks
    print('opt_tol = {:.2e}'.format(opt_tol))
    print('    rnorm = {:.2e}'.format(info['rnorm']))
    print('    iters = {:d}'.format(info['niters']))
    print('    rnorm == ||A x - b||_2: {}'.format(
        info['rnorm'] == np.linalg.norm(A.dot(x) - b)))
    print('    rnorm <= sigma: {}'.format(info['rnorm'] <= sigma))
    print('    rnorm - sigma <= opttol : {}'.format(info['rnorm'] - sigma <= opt_tol))

Note the last print.

Basically despite theory says that BPDN finds x such that the residual is smaller of equal than sigma this is not exactly the case in the code (in the original MATLAB code and so we kept the same for the python code) as you check that rnorm - sigma <= opt_tol (or rnorm <= sigma + opt_tol), this is why you sometimes see that your final condition is False: the algorithm has exited as the condition is satisfied but this is not the case if you take away opt_tol. The new condition I have put is always satisfied

Here is the condition in python code:

rerror1 = abs(aerror1) / max(1., rnorm)
...
if rerror1 <= opt_tol:
                    stat = EXIT_ROOT_FOUND # Found approx root.

and MATLAB:

rError1 = abs(aError1) / max(1,rNorm);
...
test3 = rError1     <=  optTol;
...
if test3, stat=EXIT_ROOT_FOUND;   end  % Found approx root.

[mrava87]

Hi,

I am not at my work computer where I can run my scripts, but I recall that the residual rnorm and the number of iterations did not match. Since A, b, sigma, opt_tol, and the maximum number of iterations are all set the same, I would expect the algorithm to generate the same output.

Ok no problem. Now that i know what is different I can try and monitor various parameters through iterations, will let you know.

[mrava87]

Finally you say

When I run spgl1 with the sigma keyword argument, it sometimes stops with exit status 1 (root found), but the sigma criterion is not met. I am confused why this is.

I totally agree and had the same problem in the past before digging into the code further and understand why (opt_tol). I think we should not change the python code as otherwise we are going to really start creating inconsistencies with the matlab one but I will update the doctring to make this less confusing :)

What do you think about this description of opt_tol:

Optimality tolerance. More specifically, when using basis pursuit denoise, 
the optimility condition is met when the absolute difference between the 
L2 norm of the residual and the ``sigma`` is smaller than``opt_tol``

[jhtu]

The docstring modification you suggested makes sense to me. I had dug into the source code and saw opt_tol was involved in some of the exit condition checks, but I hadn't figured out how. Your explanation makes a lot of sense to me.

When I was looking at the source code, I thought I saw opt_tol was involved in some other checks. I am asking because from what you said, it seems that I could set opt_tol = 0, in which case sigma for the BPDN problem would act as intended. However, could it be that setting opt_tol = 0 would mess up something else in the code?

By the way, I started a parallel inquiry into this with the original authors of the method and the Matlab code. I will let you know if I hear anything back from them. I understand the desire to keep the Python and Matlab libraries consistent.

[mrava87]

Hi,
great to know you got in touch with the original authors, would be great to hear what they say :)

And yes, opt_tol is used in several places as a tolerance, if the condition is below this tolerance the condition is met (I think reason why being that you will never get exactly zero when doing things numerically). I think that explaining all places where it is used may be too heavy but I see that the one I wrote can help as it really affects a condition people would expect to be met in BPDN (and check thelseves as you did). Let's see if we get any feedback and change things accordingly.

[mrava87]

Just done some comparison with matlab.

I noticed that the number of Products with A and Line search its was clearly different and wrong from our side (due to counting one less per _spg_line_curvy and _spg_line), I fixed it.

I also checked iterations with verbose=2 for both BPDN and LASSO and I see that early iterations (below 30) are exact to the number of displayed decimals. Things start to slowly become different in later iterations and of course that leads to tau updates happening at different iterations (as criteria are sometimes met at different iterations). The end results are however very similar and for real life problems I always found to get good results with both codes.

Provided the problem is not just due to different ways of rounding off in Matlab and python, we would need to check every small function one-by-one to nail down what makes this difference appear. I would probably do that but it may take some time before I have the time to commit to this. Happy if you want to go ahead and try to discover where discrepancies occur, otherwise we can just keep this Issue open until we solve this.

[jhtu]

Regarding opt_tol, I understand about having this tolerance for computations where we are driving a residual down to zero. I didn't expect this for the BPDN problem, but maybe it was easier for the authors to not have a special case for this. In any case, I think you are right that it makes sense to update the documentation, and hopefully that will help most people. I suppose in the future one could set opt_tol = 0 for the BPDN problem always, or something like that.

As far as matching Matlab, I'm not really that familiar with the SPGL1 algorithm, so I won't be of much help in hunting down the discrepancy. It is interesting that it is a slow divergence. I wonder what could be the source. But I agree that for practical purposes it probably doesn't matter that much at this point, as the sparse solutions are generally "good" sparse approximations of x either way.