qp() in Octave 3.0.0 returns result egregiously violating input constraints

Sat Apr 5 15:04:17 CDT 2008

On Apr 5, 2008, at 3:26 PM, Thomas Treichl wrote:
> Ben Abbott schrieb:
>> On Apr 5, 2008, at 2:20 PM, Joshua Redstone wrote:
>>> On Sat, Apr 5, 2008 at 10:52 AM, Ben Abbott <bpabbott at mac.com>  
>>> wrote:
>>>
>>> On Apr 5, 2008, at 1:16 PM, Joshua Redstone wrote:
>>> Another data point: for kicks, I also tried running the Octave   
>>> 2.9.18 release from sourceforge, running it directly from  
>>> the .dmg  (i.e., without installing),
>>> and I got the same kind of numerical error.
>>> Josh
>>>
>>> On Sat, Apr 5, 2008 at 10:06 AM, Joshua Redstone  <redstone at gmail.com 
>>> > wrote:
>>> It's encouraging you're getting reasonable numbers on the example  
>>> I  sent - hopefully I can too!
>>> On the question you asked, I get:
>>> octave-3.0.0:10> zzz=randn(1000); max(max(zzz*inv(zzz)))
>>> ans =  1.0000
>>>
>>> Re libraries - any suggestions how to find which one may be  
>>> broken?   I installed Octave with the stock
>>> osx binary package from http://octave.sourceforge.net/.  If the   
>>> cause is broken libraries, would they have had to come from
>>> that release?
>>> I'm trying to compile octave from source code so I can put a  
>>> bunch  of printfs in __qp__(), but am getting some make errors  
>>> about shell  options.
>>> Thanks,
>>> Josh
>>>
>>>
>>> On Fri, Apr 4, 2008 at 11:06 PM, Dmitri A. Sergatskov <dasergatskov at gmail.com
>>>> wrote:
>>> On Fri, Apr 4, 2008 at 8:24 PM, Joshua Redstone  
>>> <redstone at gmail.com>  wrote:
>>>> Hi John,
>>>> I've attached a data file created with (on OSX):
>>>>    save -binary -z "./example.dat" x0 H Q A B LB UB A_LB A_IN A_UB
>>>>
>>>> The following should reproduce the problem for you:
>>>>   load "example.dat"
>>>>   [X, OBJ, INFO, LAMBDA] = qp(x0, H, Q, A, B, LB, UB, A_LB, A_IN,
>>> A_UB)
>>>> With my qp.m that I modified so the last arguement is maxiter, I
>>> get the
>>>> follow first few lines of output:
>>>> octave-3.0.0:5> [X, OBJ, INFO, LAMBDA] = qp(x0, H, Q, A, B, LB,
>>> UB, A_LB,
>>>> A_IN, A_UB, 187)
>>>> X =
>>>>
>>>>   0.17954
>>>>    0.00000
>>>>   0.01742
>>>>    0.00304
>>>>   0.16220
>>>>    [... output truncated.... ]
>>> I cannot reproduce the problem you were describing:
>>>
>>> octave:11> X
>>> X =
>>>
>>> 1.8652e-01
>>> -1.6807e-18
>>> 1.1024e-02
>>> 2.4571e-03
>>> 1.5873e-01
>>> 4.5925e-18
>>> 3.0011e-17
>>> ....
>>>
>>> octave:12> max(X)
>>> ans =  0.18652
>>>
>>> 2.2584e-18
>>> octave:13> sum(X)
>>> ans =  1.00000
>>>
>>> I suspect that some of your library is broken.
>>> (BLAS/LAPACK is a popular candidate.)
>>> Just for kicks what do you get for:
>>>
>>> zzz=randn(1000); max(max(zzz*inv(zzz)))
>>>
>>> ?
>>>
>>> Sincerely,
>>>
>>> Dmitri.
>>>
>>> I'm running Octave 3.0.0+ on Mac OSX. My binary was built using   
>>> gfortran and Fink's package manager, and Apple's lapack/blas.
>>>
>>> octave:14> load "example.dat"
>>> octave:15>[X, OBJ, INFO, LAMBDA] = qp (x0, H, Q, A, B, LB, UB,  
>>> A_LB,  A_IN, A_UB);
>>> octave:16> X
>>>
>>> X =
>>>
>>>  2.0000e-01
>>>  2.7044e-17
>>>  7.7140e-17
>>>  9.8057e-18
>>>  5.5853e-18
>>>  4.0810e-02
>>>  2.7094e-17
>>> -1.3130e-16
>>>  8.6621e-02
>>> -6.4476e-17
>>>  9.4064e-17
>>>  1.3185e-01
>>> ...
>>>
>>> octave:17> max(X)
>>> ans =  0.20000
>>> octave:18> sum(X)
>>> ans =  1.00000
>>> octave:19> tol = 10*eps;
>>> octave:20> all ( X > LB-tol & X < UB+tol)
>>> ans =  1
>>> octave:21> all ( A_LB-tol < A_IN*X & A_UB+tol > A_IN*X)
>>> ans =  1
>>>
>>> I am unfamiliar with qp, but I this looks ok, yes?
>>>
>>> I've copied Thomas Treichl who I understand built the binary Josh  
>>> is  running.
>>>
>>> Ben
>>>
>>> Hi Ben and Dmitri,
>>> So I tried the problematic example on a version of Octave 3.0.0   
>>> compiled from source on ubuntu, and it looks like it worked fine.
>>> I got the same answer as Dmitri reported (though different from   
>>> Ben's).
>>> So this suggests the problem is particular to the osx sourceforge   
>>> package and/or my mac.
>>> I'm downloading xcode so I can try fink.
>>> Josh
>>>
>> hmmm ... I assume the different answers are related to differences  
>> in  lapack/blas.
>> Might someone check their solution against mine to determine so  
>> that I  might have an idea of how significant the difference is?
>> octave:1> load "example.dat"
>> octave:2> [X, OBJ, INFO, LAMBDA] = qp (x0, H, Q, A, B, LB, UB,  
>> A_LB,  A_IN, A_UB);
>> octave:3> tol = 10*eps;
>> octave:4> all ( X > LB-tol & X < UB+tol)
>> ans =  1
>> octave:5> all ( A_LB-tol < A_IN*X & A_UB+tol > A_IN*X)
>> ans =  1
>> octave:6> result = X'*H*X + X'*Q
>> result =  5.5352e-06
>> Ben
>
> Thanks Ben for forwarding your last email. Just read it...
>
> I just had a quick look at the problem. It seems to me that at least  
> the solution of that qr() call looks the same on my Mac as on  
> Joshua's Mac. However, I really don't know very much about this qr()  
> and I think debugging this *huge* problem is quite difficult, isn't  
> it?
>
> So how to continue from here? BTW Octave.app uses Apple's -framework  
> vecLib which is installed locally on every Mac. As it seems not to  
> happen with a self made blas/lapack like on Ben's Mac - do we need  
> another test inside acx_blas_with_f77_func for this in the default  
> branch (where I can see this result too)?

I've built octave using Apple's lapack/blas (i.e. vecLib), so the  
problem isn't directly related to the lapack/blas, correct? I

However, when I tried the release-3-0-x branch which I build yesterday  
(and passed all tests), it fails (also built with vecLib).

octave:8> result = 0.5*X'*H*X + X'*Q
result =  5.4445e+58
octave:9> tol = 10*eps;
octave:10> all ( X > LB-tol & X < UB+tol)
ans = 0
octave:11> all ( A_LB-tol < A_IN*X & A_UB+tol > A_IN*X)
ans = 0

Thomas, when did you build last? The sources for one I'm running were  
pulled from John's archive archive 2-3 days ago.

ls /sw/bin/octave
lrwxr-xr-x 1 root admin 13 2008-04-03 19:22 /sw/bin/octave ->  
octave-3.0.0+

I'll pull the most recent sources, build the default branch, and then  
test the result.

Ben