bugGNU Octave - Bugs: bug #61762, qp should check positive...

 
 

bug #61762: qp should check positive definiteness

Submitter:  Olaf Till <i7tiol>
Submitted:  Tue 04 Jan 2022 07:54:43 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Missed Error or Warning
Status:  Fixed Assigned to:  None
Originator Name:  Olaf Till Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 27 Jun 2023 09:17:29 AM UTC, comment #9: 

After some additional tests, it appears that the behaviour of _qp_ (as it is) is not reliable when H is not positive definite.

As long as a better QP solver is not used, checking for positive definiteness before calling _qp_ is therefore a safe choice.

Sorry for the noise.

Julien Bect <jbect>
Tue 27 Jun 2023 07:33:52 AM UTC, comment #8: 

Re-opening this thread in relation with https://savannah.gnu.org/bugs/?64346.

I have looked at the source code of _qp_ and I am not convinced that the algorithm is restricted to the case where H is positiv definite.  It seems to handle positive semi-definite matrices as well.

I have tried to reproduce the problem in the first post (after commenting out the test for positiv definiteness in qp.m of course):

% Objective: 1/2 x'*H*x + q'*x
H = [21, 30, 39; 30, 45, 60; 39, 60, 81];
q = [-40; -65; -90];

[x, fx, info, lambda] = qp ([], H, q)


This is a convex quadratic problem, where H is positive semi-definite, with rank 2.  As a consequence, the objective function has a minimum, but the minimizer is not unique.

Anyway, qp is able to find a solution, and no error is raised:

x =

  -1.4324
  -0.1352
   1.9010

fx = -52.500
info =

  scalar structure containing the fields:

    solveiter = 17
    info = 0

lambda = [](0x1)


Checking optimality:

>> max (abs (H*x + q))
ans = 3.7414e-10


I cannot reproduce the original error.

Julien Bect <jbect>
Sun 09 Jan 2022 12:00:10 PM UTC, comment #7: 

Ok. Closing as fixed.

Markus Mützel <mmuetzel>
Group administrator
Sat 08 Jan 2022 07:55:59 PM UTC, comment #6: 

OK I see what you mean by code duplication now. I think it's easier in that case to not move the code up and not return early.

I was using "not quadratic" in the sense of "not suitable for a quadratic minimization problem" but yes it was an ambiguous comment.

Arun Giridhar <arungiridhar>
Group Member
Sat 08 Jan 2022 06:14:37 PM UTC, comment #5: 

The order in which fields appear in a structure is important. (Or it can cause issues when trying to concatenate these structures.)
That's the kind of error that I was trying to hint at with my comment. (Though admittedly "a lot" was probably an exaggeration.)

Also in the code further down in qp.m, there are some checks for nargout. Would we need to replicate those here?

I'm still not sure if it is a good idea to replicate that code (from a maintainability point of view).
I haven't looked into it in much detail. But might that cause the function to return a value when it otherwise would error out at some of the tests that are skipped?

(Also, I haven't come across "quadratic" as a synonym for "positive definite". Do those mean the same?)

Markus Mützel <mmuetzel>
Group administrator
Wed 05 Jan 2022 07:39:49 PM UTC, comment #4: 

I was thinking of something like this, moving the check up where the other checks are done on the matrix, but if this will cause problems we don't need to do it.


diff -r 00ab5e929111 scripts/optimization/qp.m
--- a/scripts/optimization/qp.m Wed Jan 05 19:30:43 2022 +0100
+++ b/scripts/optimization/qp.m Wed Jan 05 14:36:42 2022 -0500
@@ -182,6 +182,13 @@
   ## Validate the quadratic penalty.
   if (! issquare (H))
     error ("qp: quadratic penalty matrix must be square");
+  elseif (isdefinite (H) != 1)  # not quadratic, return early
+    x = x0;
+    obj = nan;  # or something indicating problem instance is unbounded below
+    INFO.info = 2;
+    INFO.solveiter = 0;
+    lambda = [];
+    return
   elseif (! ishermitian (H))
     ## warning ("qp: quadratic penalty matrix not hermitian");
     H = (H + H')/2;
@@ -348,10 +355,6 @@

   info = 0;

-  if (isdefinite (H) != 1)
-    info = 2;
-  endif
-
   if (info == 0 && (eq_infeasible || in_infeasible))
     ## The initial guess is not feasible.
     ## First, define an xbar that is feasible with respect to the


Arun Giridhar <arungiridhar>
Group Member
Wed 05 Jan 2022 10:50:00 AM UTC, comment #3: 

@Arun: I don't see how that could be done without duplicating a lot of code. But maybe I'm missing your point...

Markus Mützel <mmuetzel>
Group administrator
Tue 04 Jan 2022 02:44:23 PM UTC, comment #2: 

I notice that there is input validation of the matrix near line 185. Would it be better to check for definiteness there and return early setting info=2 etc? That way it can skip the rest of the validation code.

Arun Giridhar <arungiridhar>
Group Member
Tue 04 Jan 2022 01:15:18 PM UTC, comment #1: 

Thanks for the patch. Looks good to me.
I added your motivating example as a BIST and pushed here:
https://hg.savannah.gnu.org/hgweb/octave/rev/2f1fae9dd79d

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Tue 04 Jan 2022 07:54:43 AM UTC, original submission:  


qp ([], [21, 30, 39; 30, 45, 60; 39, 60, 81], [-40; -65; -90])
error: __qp__: operator *: nonconformant arguments (op1 is 2x2, op2 is 3x1)
error: called from
    qp at line 409 column 25


Running under gdb shows the exception is thrown by the last of the following lines of _qp_.cc (lines 200--213):

              // Inverting the Hessian.  Using the Cholesky
              // factorization since the Hessian is positive
              // definite.

              math::chol<Matrix> cholH (H);

              R = cholH.chol_matrix ();

              Matrix Hinv = math::chol2inv (R);

              // Computing the unconstrained step.
              // p = -Hinv * g;

              p = -Hinv * g;


So H, given as an argument to _qp_, is required to be positive
definite. qp.m should check H to be positive definite before passing it to _qp_. If the check fails, a suitable error flag should be returned. A patch will be submitted as soon as a number is assigned to this bug report.

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

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -email is unavailable- added by jbect (Posted a comment)
  • -email is unavailable- added by siko1056
  • -email is unavailable- added by arungiridhar (Posted a comment)
  • -email is unavailable- added by mmuetzel (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
    2023-11-16 arungiridhar Dependencies- bugs #64901 is dependent
    2022-01-09 mmuetzel StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2022-01-05 siko1056 Carbon-Copy- Added siko1056
    2022-01-04 mmuetzel StatusNone Ready For Test
    2022-01-04 i7tiol Attached File- Added qp-check-positive-definite.patch, #52610

    Back to the top

    Powered by Savane 3.13-bb6a.
    Corresponding source code