bugGNU Octave - Bugs: bug #56037, qp exits with an indexing error...

 
 

bug #56037: qp exits with an indexing error for a specific input with OpenBLAS

Submitter:  None
Submitted:  Sat 30 Mar 2019 04:34:33 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Unexpected Error or Warning
Status:  Fixed Assigned to:  None
Originator Name:  Rohit Beriwal Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 5.1.0
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 04 Apr 2019 09:20:37 PM UTC, comment #21: 

You're right that it is probably down to the different eigenvectors.  With most iterative optimization algorithms the eventual local minimum found depends strongly on the starting position.  I accept that.  But when you start at the same place, and use the same algorithm, I find it odd to arrive at different destinations.

 Rik <rik5>
Group administrator
Thu 04 Apr 2019 07:57:59 PM UTC, comment #20: 

The last part should have read

With OpenBLAS the eigenvectors are


[v, ~] = eig (p)
v =

  -0.22089  -0.02616  -0.00000   0.16667   0.16667   0.83333  -0.27481  -0.35355
   0.22089   0.02616   0.00000  -0.16667   0.83333   0.16667   0.27481   0.35355
   0.22089   0.02616   0.00000   0.83333  -0.16667   0.16667   0.27481   0.35355
   0.22089   0.02616   0.00000  -0.50000  -0.50000   0.50000   0.27481   0.35355
  -0.11281   0.44394   0.70711  -0.00000  -0.00000   0.00000  -0.40643   0.35355
  -0.11281   0.44394  -0.70711   0.00000   0.00000  -0.00000  -0.40643   0.35355
  -0.85950  -0.26116   0.00000  -0.00000  -0.00000   0.00000   0.26086   0.35355
   0.20156  -0.73136   0.00000   0.00000   0.00000   0.00000  -0.54725   0.35355


And with netlib BLAS, the eigenvectors are returned as


[v, ~] = eig (p)
v =

  -0.07138   0.16667   0.16667   0.83333   0.22361   0.01465  -0.26399  -0.35355
   0.07138   0.83333  -0.16667   0.16667  -0.22361  -0.01465   0.26399   0.35355
   0.07138  -0.16667   0.83333   0.16667  -0.22361  -0.01465   0.26399   0.35355
   0.07138  -0.50000  -0.50000   0.50000  -0.22361  -0.01465   0.26399   0.35355
   0.07138  -0.00000  -0.00000   0.00000   0.89443  -0.01465   0.26399   0.35355
  -0.45433  -0.00000  -0.00000   0.00000  -0.00000  -0.71838  -0.39053   0.35355
  -0.56908   0.00000   0.00000  -0.00000   0.00000   0.68695  -0.28152   0.35355
   0.66653   0.00000   0.00000   0.00000   0.00000   0.10468  -0.64791   0.35355


The first 7 columns correspond to eigenvalues of zero or near to zero, so is it reasonable that the eigenvectors are not identical? Seems reasonble to me.

I have no idea if there is a bug in how the qp algorithm is implemented when the system has one or more zero eigenvalues.

 Mike Miller <mtmiller>
Group Member
Thu 04 Apr 2019 07:37:13 PM UTC, comment #19: 

I get the following three solutions with netlib, OpenBLAS, and ATLAS, now that this patch is applied.

With OpenBLAS


x =

   0.12500
   0.00000
   0.00000
   0.66721
   0.00000
   0.00000
   0.10390
   0.10390

obj =  0.56250
info =

  scalar structure containing the fields:

    solveiter =  5
    info = 0


With netlib BLAS


x =

   0.12500
   0.23988
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.63512

obj =  0.56250
info =

  scalar structure containing the fields:

    solveiter =  6
    info = 0


With ATLAS


x =

   0.12500
   0.23988
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.63512

obj =  0.56250
info =

  scalar structure containing the fields:

    solveiter =  6
    info = 0


The difference between BLAS and OpenBLAS can also be seen by simply getting the eigenvectors of the system. The different qp results are likely due entirely to the different eigenvectors returned by the selected BLAS library.

With OpenBLAS


[v, ~] = eig (p)
v =

  -0.22089  -0.02616  -0.00000   0.16667   0.16667   0.83333  -0.27481  -0.35355
   0.22089   0.02616   0.00000  -0.16667   0.83333   0.16667   0.27481   0.35355
   0.22089   0.02616   0.00000   0.83333  -0.16667   0.16667   0.27481   0.35355
   0.22089   0.02616   0.00000  -0.50000  -0.50000   0.50000   0.27481   0.35355
  -0.11281   0.44394   0.70711  -0.00000  -0.00000   0.00000  -0.40643   0.35355
  -0.11281   0.44394  -0.70711   0.00000   0.00000  -0.00000  -0.40643   0.35355
  -0.85950  -0.26116   0.00000  -0.00000  -0.00000   0.00000   0.26086   0.35355
   0.20156  -0.73136   0.00000   0.00000   0.00000   0.00000  -0.54725   0.35355



[v, ~] = eig (p)
v =

  -0.07138   0.16667   0.16667   0.83333   0.22361   0.01465  -0.26399  -0.35355
   0.07138   0.83333  -0.16667   0.16667  -0.22361  -0.01465   0.26399   0.35355
   0.07138  -0.16667   0.83333   0.16667  -0.22361  -0.01465   0.26399   0.35355
   0.07138  -0.50000  -0.50000   0.50000  -0.22361  -0.01465   0.26399   0.35355
   0.07138  -0.00000  -0.00000   0.00000   0.89443  -0.01465   0.26399   0.35355
  -0.45433  -0.00000  -0.00000   0.00000  -0.00000  -0.71838  -0.39053   0.35355
  -0.56908   0.00000   0.00000  -0.00000   0.00000   0.68695  -0.28152   0.35355
   0.66653   0.00000   0.00000   0.00000   0.00000   0.10468  -0.64791   0.35355


 Mike Miller <mtmiller>
Group Member
Thu 04 Apr 2019 07:20:11 PM UTC, comment #18: 

@Mike: Is there another bug report here in that, for the same problem, qp reaches a different solution based on the particular BLAS backend used?

 Rik <rik5>
Group administrator
Thu 04 Apr 2019 06:35:15 PM UTC, comment #17: 

I made this additional change on the stable branch to clean this up for safety and readability

https://hg.savannah.gnu.org/hgweb/octave/rev/4a2cb3392014

Closing as fixed.

 Mike Miller <mtmiller>
Group Member
Wed 03 Apr 2019 08:00:13 PM UTC, comment #16: 

It's a small point, but in the test script try1.m the order of output arguments is incorrect.  "info" and "obj" (value of objective function at final solution) should be reversed.

It should read


[x,obj,info,lambda] = qp([],p,q,a,b,lb,ub)


Apparently I am using reference BLAS because I also get a final value for the objective function of 0.5625.  This is the case even when I start very close to the other solution.


[x,obj,info,lambda] = qp([0.5 0 0.51 0 0 0 0 0],p,q,a,b,lb,ub)


I would suggest inspecting the variables before and after line 413


    [x, lambda, info, iter] = __qp__ (x0, H, q, A, b, Ain, bin, maxit);


when using a working and a non-working BLAS library.


Rik <rik5>
Group administrator
Wed 03 Apr 2019 07:22:32 PM UTC, comment #15: 

I do not get error code 3 any more, but the answers are different depending on blas version. Perhaps the problem is ill-posed...


(default, with openblas)
octave -q -f
octave:1> try1
x =

   0.50000
   0.00000
   0.50000
   0.00000
   0.00000
  -0.00000
   0.00000
   0.00000

info =    1.4050e-31
obj =

  scalar structure containing the fields:

    solveiter =  8
    info = 0

lambda =

  -1.7790e-47
   0.0000e+00
   0.0000e+00
   4.9304e-31
   0.0000e+00
   0.0000e+00
   0.0000e+00
   4.9304e-31
   0.0000e+00
   4.9304e-31
   0.0000e+00
   4.9304e-31
   0.0000e+00
   4.9304e-31
   0.0000e+00
   4.9304e-31
   0.0000e+00


With atlas:


LD_PRELOAD=/usr/lib64/atlas/libtatlas.so octave -q -f
octave:1> try1
x =

   0.50000
  -0.00000
   0.00000
   0.00000
   0.35884
  -0.00000
   0.00000
   0.14116

info =    3.0830e-33
obj =

  scalar structure containing the fields:

    solveiter =  14
    info = 0

lambda =

   4.5213e-17
   0.0000e+00
   0.0000e+00
   4.5213e-17
   0.0000e+00
   4.5213e-17
   0.0000e+00
   0.0000e+00
   0.0000e+00
   0.0000e+00
   0.0000e+00
   4.5213e-17
   0.0000e+00
   4.5213e-17
   0.0000e+00
   0.0000e+00
   0.0000e+00


with reference blas:


LD_PRELOAD=/usr/lib64/libblas.so.3 octave -q -f
octave:1> try1
x =

   0.12500
   0.23988
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.63512

info =  0.56250
obj =

  scalar structure containing the fields:

    solveiter =  6
    info = 0

lambda =

   0.50000
   0.00000
   0.00000
   0.00000
   0.00000
   1.00000
   0.00000
   1.00000
   0.00000
   1.00000
   0.00000
   1.00000
   0.00000
   1.00000
   0.00000
   0.00000
   0.00000


Dmitri.

Dmitri A. Sergatskov <dasergatskov>
Wed 03 Apr 2019 07:02:48 PM UTC, comment #14: 

That's probably why it just got coded as a for loop then.  The only optimization I see would be to use xelem instead of () for indexing the columnvector


if (mRrH == eigenvalH(i))
=>
if (mRrH == eigenvalH.xelem (i))


Since the range of i is already known to be within the bounds of the columnvector there isn't a requirement to check for that each time a value is extracted.

I can't test the patch because the code doesn't error on my system.

Rik <rik5>
Group administrator
Wed 03 Apr 2019 06:38:36 PM UTC, comment #13: 

Ok, I pushed this fix on the stable branch

https://hg.savannah.gnu.org/hgweb/octave/rev/2f3e8ea7e387

I don't think the min approach is as simple as that, there doesn't seem to be a 'min' method in the Array<T> base class. Only the Matrix<T> and NDArray<T> classes have a 'min' method that takes idx and dim arguments. The ColumnVector and RowVector classes only have a 'min(void)' method.

 Mike Miller <mtmiller>
Group Member
Wed 03 Apr 2019 05:49:13 PM UTC, comment #12: 

Regarding getting rid of the for loop...

One thing would be to call the Octave min() function (available as Fmin) from C++ since the second return argument is the index that you seek.  But, this does involve wrapping the inputs in octave_value_list and then unpacking the outputs that are also in an octave_value_list.  It probably isn't worth that complication.

Another possibility is to use the C++ code that min() eventually calls.


Array<octave_idx_type> idx;
min = array.min (idx, dim);


If I understand it correctly, eigenvalrH is a column vector so the declaration could be just


octave_idx_type idx;
min = eigenvalrH.min (idx, dim);



 Rik <rik5>
Group administrator
Wed 03 Apr 2019 04:19:54 PM UTC, comment #11: 

Nice that the bug fix is so simple in the end.

For the commit message, my understanding of the GNU standards which Octave has chosen to follow is

1) First sentence is a concise summary of the overall changeset
 1a) Octave-specific modification is that if the change occurred in an m-file the name of the m-file followed by a colon starts the line.
 1b) Octave occasionally use the prefixes "doc:", "maint:", "test:" if the changeset is for one of these purposes.
2) Additional sentences describe more detail about the overall change, perhaps including reasoning for a particular approach over another.
3) Changes on a per-file basis document what was done, not why it was done (that is the global reasoning in parts 1 & 2)

So for the commit message I might try something like


Fix obscure corner case when calculating qp caused by a typo (bug #56037).

* __qp__.cc (qp): Fix typo where "eigenvalH" was used instead of "eigenvalrH".



 Rik <rik5>
Group administrator
Wed 03 Apr 2019 02:56:47 AM UTC, comment #10: 


> it will probably take some gdb and/or printf-debugging to see where it comes from.


I did this now.

Attached is vars.txt, which is a set of local variables extracted from the 'qp.m' function file at the point of the bad call to the '__qp__' function. The error is reproducible for me with


./run-octave -cli --eval "load vars.txt; __qp__ (x0, H, q, A, b, Ain, bin, maxit);"


This error message comes from the octave::out_of_range exception in liboctave, so I simply set a breakpoint on the exception constructor, and the error is raised in '__qp__.cc' line 313 (on the stable branch):

https://hg.savannah.gnu.org/hgweb/octave/file/c3716220d5b9/libinterp/corefcn/__qp__.cc#l313

Which led me to discover a typo in that section of code. The attached changeset fixes this on the stable branch. I could use some help wording the commit message, since this is not my area of expertise. It would also be better if we had a more well-tested and reliable function to get both the minimum value and index from a vector instead of using these ad-hoc loops.

(file #46703, file #46704)

 Mike Miller <mtmiller>
Group Member
Tue 02 Apr 2019 11:47:00 PM UTC, comment #9: 

Yeah, I did that earlier, but I guess I didn't post the details.

The error comes from inside '__qp__', so it will probably take some gdb and/or printf-debugging to see where it comes from.

 Mike Miller <mtmiller>
Group Member
Tue 02 Apr 2019 11:37:40 PM UTC, comment #8: 

@Mike: If I look at line 413 of qp.m in 5.1.0 I see


    [x, lambda, info, iter] = __qp__ (x0, H, q, A, b, Ain, bin, maxit);


But your error message refers to a variable index.  In fact, I can't find "index" anywhere in qp.m.  Could you set


debug_on_error (1)


and then run your failing example from comment #3 and find out what is going on?

Rik <rik5>
Group administrator
Tue 02 Apr 2019 08:12:00 PM UTC, comment #7: 

I got z=0.5 solution when i use atlas:


LD_PRELOAD=/usr/lib64/atlas/libtatlas.so octave -f -q
octave:1> try1
x =

   0.50000
   0.50000
   0.00000
  -0.00000
   0.00000
   0.00000
   0.00000
   0.00000

info =    3.7309e-34
obj =

  scalar structure containing the fields:

    solveiter =  9
    info = 0

lambda =

   0
   0
   0
   0
   0
   0
   0
   0
   0
   0
   0
   0
   0
   0
   0
   0
   0


Dmitri.

Dmitri A. Sergatskov <dasergatskov>
Tue 02 Apr 2019 05:06:07 PM UTC, comment #6: 

@Rohit Notice that obj.info is 3, which means the maximum number of iterations was reached. So the solution that Dmitri gets, and the different solution that I get with BLAS, explicitly indicates that it has not really found a solution. It's up to the caller to check obj.info to see that a global solution was not reached. But that is not what this bug is about.

 Mike Miller <mtmiller>
Group Member
Tue 02 Apr 2019 04:56:06 PM UTC, comment #5: 

Just to help a bit more, by solving it by hand, one of the possible answer is (x1=0.5, x2=0.5, x3=x4=…=x8=0). The objective value of the function is ZERO (the global minimum) and not 0.62229 as calculated by @Dmitir. Since the function is not reaching the minimum, I suppose something is really wrong with it.

(Also, I used whatever default implementation of BLAS that comes with Octave 4.4.1, so mostly this bug should be reproducible in that scenario for sure, cant say about other BLAS or OS implementation.)

Thank you

 Rohit Beriwal <rohitberiwal>
Tue 02 Apr 2019 02:43:49 PM UTC, comment #4: 

The results in my previous post were done
with openblas.
I get different answers when using ref blas and
different again when using atlas, but it would not
crash for me.

Dmitri.

Dmitri A. Sergatskov <dasergatskov>
Tue 02 Apr 2019 02:37:29 PM UTC, comment #3: 

I can confirm the error reported here with both Octave 4.4 and 5 using OpenBLAS


p = 2 * [-1; ones(7,1)] * [-1, ones(1,7)];
q = zeros (8,1);
lb = zeros (1,8);
ub = ones (1,8);
a = ones (1,8);
b = 1;
[x, info, obj, lambda] = qp ([], p, q, a, b, lb, ub)


error: index (_,4): but object has size 3x3
error: called from
    qp at line 413 column 29:


If I use BLAS instead, I get a solution, but a different result than the one shown by Dmitri in comment #2.

 Mike Miller <mtmiller>
Group Member
Tue 02 Apr 2019 05:39:18 AM UTC, comment #2: 

Seems to work for me with octave 5.1.1:


 try1
x =

   0.10557
   0.00000
   0.00000
   0.22361
   0.22361
   0.22361
   0.00000
   0.22361

info =  0.62229
obj =

  scalar structure containing the fields:

    solveiter =  200
    info =  3

lambda =

   0.94663
   0.00000
   0.00000
   0.63108
   0.00000
   0.63108
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.63108
   0.00000
   0.00000
   0.00000

octave:2> x
x =

   0.10557
   0.00000
   0.00000
   0.22361
   0.22361
   0.22361
   0.00000
   0.22361


Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Sat 30 Mar 2019 05:00:45 PM UTC, comment #1: 

I don't see the file

 Doug Stewart <dastew>
Sat 30 Mar 2019 04:34:33 PM UTC, original submission:  

Current OS: Windows 10, Build no: 1809,
System: 16GB RAM, i7-7700HQ,
Octave version: 4.4.1(G.U.I)

I was trying to solve the a minimization problem (file linked below). And I Have rechecked each and every matrix size and input syntax for QP() function, but for my particular problem, it gives the following error:

error: index (_,4): but object has size 3x3
error: called from
    qp at line 409 column 25
    try1 at line 14 column 20

 Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #46703:  vars.txt added by mtmiller (1KiB - text/plain)
file #46704:  qp.diff added by mtmiller (918B - text/x-patch)
file #46693:  try1.m added by None (419B - 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 rohitberiwal (Posted a comment)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by dastew (Posted a comment)
  • -email is unavailable- added by None (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 13 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-04-04 mtmiller StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2019-04-03 mtmiller StatusPatch Reviewed Ready For Test
    2019-04-03 rik5 StatusPatch Submitted Patch Reviewed
    2019-04-03 mtmiller Attached File- Added vars.txt, #46703
        Attached File- Added qp.diff, #46704
        StatusConfirmed Patch Submitted
    2019-04-02 rik5 Release4.4.1 5.1.0
        Operating SystemMicrosoft Windows Any
    2019-04-02 mtmiller StatusNeed Info Confirmed
        SummaryQp function crashes for a specific input qp exits with an indexing error for a specific input with OpenBLAS
    2019-04-02 None Attached File- Added try1.m, #46693
    2019-03-30 mtmiller StatusNone Need Info

    Back to the top

    Copyright © 2024 Free Software Foundation, Inc.
    Verbatim copying and distribution of this entire article is permitted in any medium, provided this notice is preserved.
    The Levitating, Meditating, Flute-playing Gnu logo is a GNU GPL'ed image provided by the Nevrax Design Team.
    Page source code

    Powered by Savane 3.14-f13d.
    Corresponding source code