bugGNU Octave - Bugs: bug #31742, sqp: incorrect order of removal of...

 
 

bug #31742: sqp: incorrect order of removal of Inf bounds

Submitter:  Olaf Till <i7tiol>
Submitted:  Thu 02 Dec 2010 11:55:12 AM UTC
   
 
Category:  None Severity:  3 - Normal
Priority:  5 - Normal Item Group:  None
Status:  Fixed Assigned to:  None
Originator Name:  Olaf Till Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 29 Dec 2010 07:23:22 PM UTC, comment #10: 

This issue has been resolved.  A second bug report is being filed to report new issues discovered with sqp while implementing this change.

Rik <rik5>
Group administrator
Thu 16 Dec 2010 03:45:27 PM UTC, comment #9: 

I credited you for the patch in the more important ChangeLog file kept within Mercurial.  I reverted the ",1" indexing change I made in this changeset (http://hg.savannah.gnu.org/hgweb/octave/rev/c767bb1afa03).

I think you need to keep the tmp_bnds variable to handle the case of different numbers of bounds.


lb = [1; 2; 3];
ub = [10; Inf; 20];

There is no upper constraint on the second variable, which leads to __sqp_lb__ being of size 3x1 and __sqp_ub__ beof of size 2x1 and the comparison "any (__sqp_lb__ > __sqp_ub__)" would fail.


This could be handled by calculating the index and not applying it until after the comparison above.  However, I think the code would be less clear.  It is easier to understand calculating the index and removing the inapplicable constraints as one block of code.

Rik <rik5>
Group administrator
Wed 15 Dec 2010 07:45:59 AM UTC, comment #8: 

(I did not find the applied changeset first because it has the wrong author logged.)

You have edited it and made a change that should be revoked: the lines

__sqp_lb__ = __sqp_lb__(tmp_idx);
...
__sqp_ub__ = __sqp_ub__(tmp_idx);

must indeed be

__sqp_lb__ = __sqp_lb__(tmp_idx, 1);
...
__sqp_ub__ = __sqp_ub__(tmp_idx, 1);

as I had it. Otherwise, if lb (or ub, respectively) is scalar and -Inf (or +Inf, respectively), _sqp_lb_ (or _sqp_ub_) get dimensions 0x0, which are the wrong dimensions for the later calculations.

I also saw now that I did not see that actually tmp_lb and tmp_ub are not needed. They can be removed everywhere and the test

if (any (tmp_lb > tmp_ub))

can be made with _sqp_lb_ and _sqp_ub_ instead. (Since the removed bounds still present in tmp_lb can not be greater than the respective bounds in tmp_ub anyway, and vice versa.)

Olaf Till <i7tiol>
Group Member
Mon 13 Dec 2010 07:28:29 PM UTC, comment #7: 

Excellent testing.  I always found the global variables here a bit of a hack, but had never gone to the trouble of actually testing their performance.  I applied your patch here http://hg.savannah.gnu.org/hgweb/octave/rev/2726132f77f6.

While testing your patch with the example in the documentation I found a few more problems with sqp (not your patch).  First, it doesn't always respect bounds if only one bound is specified, such as -1,[].  It seems that in this case the algorithm converges on a value for x(1) near -1.7 and the algorithm stops because the stepsize becomes too small.  The return value for this case is still 101, success, which it clearly isn't.  Perhaps we should use a different return value to indicate this sort of stopping condition?  Secondly, the problem seems to be with the magnitude of the upper bound in this case.  If I use a smaller upper bound such as 10, rather than +realmax, then sqp works.  In the same vein, if I use a vector of all Infs, then sqp also works because the constraint is removed with your new code.  Perhaps, if the bound is unspecified, the code should reduce to using +/-Inf rather than +/-realmax so that the bound is never tested?

The second problem is a little more serious.  If I run the demo code and set a smallish upper bound then the function crashes because lambda is empty.

[x, obj, info, iter, nf, lambda] = sqp (x0, @phi, @g, [], -10, 1)
error: sqp: operator *: nonconformant arguments (op1 is 5x13, op2 is 0x0)
error: called from:
error:   /home/rik/wip/Projects_Mine/octave-doc/scripts/optimization/sqp.m at line 481, column 9

Since you seem to understand the algorithm pretty well it would help if you took a look at that.  Be sure to pull a new copy of sqp.m from the repository since I also changed a bit of the documentation which had repercussions for some variable names in the code.

Rik <rik5>
Group administrator
Mon 13 Dec 2010 09:35:56 AM UTC, comment #6: 

Thanks for the thorough review.

The first issue (scalar bounds) is resolved in the now attached patch (without expanding the scalar bounds).

The third suggestion (isinf() ->  ==) is also in the new patch. I did not know that == is overloaded for +-Inf in current Octave, it was not so before.

I did not follow your second suggestion (global variables) for the following reasons:

1. I don't think global variables give a speed advantage, see this test:


octave:1> function f_global (x)
> global y;
> for id = 1:1000
> x == y;
> endfor
> endfunction
octave:2> function f_anonymous (x, y)
> for id = 1:1000
> x == y;
> endfor
> endfunction
octave:3> global y;
octave:4> y = 1;
octave:5> z = 1;
octave:6> af = @ (x) f_anonymous (x, z);
octave:7> tic; f_global (1), toc
Elapsed time is 0.00171609 seconds.
octave:8> tic; af (1), toc
Elapsed time is 0.00165111 seconds.
octave:9> tic; f_global (1), toc
Elapsed time is 0.00172499 seconds.
octave:10> tic; af (1), toc
Elapsed time is 0.00164196 seconds.
octave:11> tic; f_global (1), toc
Elapsed time is 0.00170795 seconds.
octave:12> tic; af (1), toc
Elapsed time is 0.00164512 seconds.
octave:13>


2. I think the original author used global variables simply because anonymouns function were not available in Octave at this time.

3. The most important point: Global variables break recursive calls of sqp. But recursive optimization is sometimes necessary. I think all global variables in sqp have to be replaced by using anonymous functions in the future. So I would not add a further global variable now.


(file #22192)

Olaf Till <i7tiol>
Group Member
Sun 12 Dec 2010 03:57:50 PM UTC, comment #5: 

I finally had time this weekend to review your patch.  I think your second approach is the the correct one.  I have a few comments that I hope you can incorporate into a final patch.

First, using scalars for upper and lower bounds no longer works.  That functionality could be restored by always extending any scalar bound to a vector the size of x.  But perhaps you want to take a look at it and see if there is a better approach?

Second, the original author seems to have used global variables to pass information to subfunctions.  Presumably this was done for efficiency since each variable in an ordinary function call would have to be pushed onto the stack and popped back off.  I would recommend adopting this strategy for your function declarations:

     ci_fun = @ (x) cf_ub_lb (x, lower_bounds_idx, upper_bounds_idx);
     ci_grd = @ (x) cigrad_ub_lb (x, bounds_grad);
                to
     ci_fun = @cf_ub_lb;
     ci_grd = @cigrad_ub_lb;

The original author used _sqp_lb_ for the lower bounds so something like _sqp_lbidx_ would make sense.

Third, I know the original author used isinf() but it is less efficient and less clear than using a direct test for Inf.  I would change the following code.

(isinf (lb) & lb < 0)
        to
(lb == -Inf)

The original construction requires iterating over the length of lb 3 times (isinf, lb < 0, '&') whereas the second construction makes only one pass and clearly shows what you are testing for.


Rik <rik5>
Group administrator
Mon 06 Dec 2010 02:14:37 PM UTC, comment #4: 

A new patch is attached.

1. Never violated Inf bounds are now already excluded from computation in the function for inequality constrained. This is better than in the first patch. The matter of checking for inequality constraints being +Inf in each iteration is not touched by this change (but see 2.). This fixes the original bug.

2. Removing inequality constraints getting +Inf in each iteration is dropped for the following reasons. Although such a removal might make some sense (see previous comments), it can not have worked and I see no fix. Tests for this are difficult to construct. According to the code, such removal would trigger the same problem as the original bug, since it produces a dimension mismatch (line 478 in the original version) between lambda (which is shortened by the removal) and (A_new - A) (which is not shortened by the removal). It would be no solution to also shorten A_new (which later is assigned to A), since A_new and A would then have different dimensions if the number of inequality constraints being +Inf changes between iterations, so (A_new - A) is not defined then. (The dimensions might even be equal while the rows refer to different constraints!)

(I noticed that my editor introduced some tabs for indentation. Please tell me if this is a problem. I faintly remember that there has been some un-tabbifying in Octave.)

(file #22132)

Olaf Till <i7tiol>
Group Member
Sun 05 Dec 2010 08:07:55 PM UTC, comment #3: 

I just noticed the mistake in my thinking. Sorry. The previous code checked for h(x) being +Inf, so only for non-violated constraints which are Inf. So it would indeed make some sense to weed them out each iteration, even if I still doubt that this was the intention since the original comment mentioned bounds.

I'll think about that and maybe supply a new patch, but not today.

Olaf Till <i7tiol>
Group Member
Sun 05 Dec 2010 07:09:31 PM UTC, comment #2: 

To avoid misunderstandings: Of corse the value of the inequality
constraints is still checked in each iteration. Only isinf() is
checked only once, and the yielded index is applied unchanged in each iteration.

Elements of h(x0) which are +Inf appear because some elements of x were specified by the user to have no bounds by setting these elements bounds to -Inf or +Inf. This means that these elements of h(x) are actually no constraints, so they can be removed forever (the original comment at these lines in sqp.m already mentions bounds).

The algorithm requires a gradient of the constraint function for
sensible action. Thats why h(x) is required to be continuous.

In particular, if elements of h(x) are Inf or -Inf, their gradient can not be computed at all. One could weed out such values if they appear, by checking isinf(h(x)) in each iteration. Under certain circumstances, this weeding out could lead to finally achieving a solution nevertheless, but only by chance (if h(x) gets non-Inf again by chance), and not by the intentional course of the algorithm. I would not suggest to do that, I'd say such a h(x) is just erroneous. Furthermore, this was clearly also not the intention before my patch, since the check was not just for isinf(h(x)), but for (isinf(h(x)) & h(x) > 0), which is true in the no-bounds specification (see 1st paragraph); also, see the original comment:

## ... never be active

(also true for the no-bounds specification, but not in our just
constructed case of h(x) sometimes returning Inf or -Inf).

Olaf Till <i7tiol>
Group Member
Sun 05 Dec 2010 03:16:24 AM UTC, comment #1: 

I have a question about your patch.  You remove inequality constraints that are -Inf at the initial evaluation point X0.  This is done just once for the entire algorithm.  The old code updated which constraints were valid at every iteration of the loop so that as the search point, X, moved some inequality constraints might become valid or invalid.  Why are you certain that the inequality constraints won't change?  I see you reference continuity of h(x), but I'd like your thoughts on how this applies.

I know sqp uses a linear model for the constraints and a linear model is obviously a continuous model.  However, everything I see is also about applying the model near the solution point.  It seems that far from the solution point one might have discontinuous functions.  For example, what if my objective function was just sin(x).  The minimum (-1) is at ...-5pi/2, -pi/2, 3*pi/2, ...  It would be easy to choose the minimum around 3*pi/2 by using the lower bound/upper bound syntax but what if I instead wrote a constraint inequality function that was -1 (not met) outside of the region [pi, 2*pi] and Inf (completely met) within that region.  It seems like this constraint would be valid and would need to be switched on or off as X moved from outside the region to within the region.

Rik <rik5>
Group administrator
Thu 02 Dec 2010 11:55:12 AM UTC, original submission:  

This is a repost after data removal due to an intrusion into savannah at 2010-11-22/24. The original report was from 2010-11-23 and had the number 31747.

This bug was reported to the former bug-list, together with a patch, in february 2010, but the patch was never applied. I would have brought the issue to the tracker earlier, but I thought I had tested that some later changes by Rik fixed it --- but there was obviously some mistake in my testing, since I see now that the bug and the respective parts of the code are the same as before. So here is it again:

There is a bug in the order in which sqp removes never-active Inf bounds, triggered e.g. by (I had moved the file to a non-standard directory for testing):


octave:3> [x, objf] = sqp ([1; 1], @ (x) sumsq ([1, 2, 4, 8, 20] - x(1) * exp
(x(2) * (1:5))), [], [], [-1; -inf], [2; 2])
error: operator *: nonconformant arguments (op1 is 2x4, op2 is 3x1)
error: called from:
error:   /home/olaf/devel/octave/optim/sqp.m at line 478, column 9


A patch is attached (also includes a minor correction:

if (__sqp_lb__ > __sqp_ub__) ->
if (any (__sqp_lb__ > __sqp_ub__))

).

Olaf Till <i7tiol>
Group Member

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #22192:  sqp-3rd-inf-bounds-removal.changeset added by i7tiol (4KiB - application/octet-stream)
file #22132:  sqp-2nd-inf-bounds-removal.changeset added by i7tiol (4KiB - application/octet-stream)
file #22087:  sqp-inf-bounds-removal.changeset added by i7tiol (3KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by i7tiol (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2010-12-29 rik5 StatusNeed Info Fixed
        Open/ClosedOpen Closed
    2010-12-13 i7tiol Attached File- Added sqp-3rd-inf-bounds-removal.changeset, #22192
    2010-12-06 i7tiol Attached File- Added sqp-2nd-inf-bounds-removal.changeset, #22132
    2010-12-05 rik5 StatusNone Need Info
    2010-12-02 i7tiol Attached File- Added sqp-inf-bounds-removal.changeset, #22087

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code