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

Sat Apr 5 14:26:24 CDT 2008

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 Ma. 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)?

Ben, here you are

octave-3.0.0:7> [X, OBJ, INFO, LAMBDA] = qp (x0, H, Q, A, B, LB, UB, A_LB, A_IN, 
A_UB);
octave-3.0.0:8> tol = 10*eps
tol =  2.2204e-15
octave-3.0.0:9> all ( X > LB-tol & X < UB+tol)
ans = 0
octave-3.0.0:10> all ( A_LB-tol < A_IN*X & A_UB+tol > A_IN*X)
ans = 0
octave-3.0.0:11> result = X'*H*X + X'*Q
result =  1.1752e+59

Further ideas?

   Thomas