Numpy.dot Fehler? Inkonsistentes NaN-Verhalten

Question

Ich bemerkte ein inkonsistentes Verhalten in numpy.dot, wenn nan s und Nullen beteiligt sind.

Kann jemand einen Sinn daraus machen? Ist das ein Fehler? Ist das spezifisch für die dot Funktion?

Ich benutze numpy v1.6.1, 64bit, läuft auf Linux (auch auf v1.6.2 getestet). Ich habe auch auf v1.8.0 auf Windows 32bit getestet (so kann ich nicht sagen, ob die Unterschiede auf die Version oder Betriebssystem oder Arch zurückzuführen sind).

from numpy import * 
0*nan, nan*0 
=> (nan, nan) # makes sense 

#1 
a = array([[0]]) 
b = array([[nan]]) 
dot(a, b) 
=> array([[ nan]]) # OK 

#2 -- adding a value to b. the first value in the result is 
#  not expected to be affected. 
a = array([[0]]) 
b = array([[nan, 1]]) 
dot(a, b) 
=> array([[ 0., 0.]]) # EXPECTED : array([[ nan, 0.]]) 
# (also happens in 1.6.2 and 1.8.0) 
# Also, as @Bill noted, a*b works as expected, but not dot(a,b) 

#3 -- changing a from 0 to 1, the first value in the result is 
#  not expected to be affected. 
a = array([[1]]) 
b = array([[nan, 1]]) 
dot(a, b) 
=> array([[ nan, 1.]]) # OK 

#4 -- changing shape of a, changes nan in result 
a = array([[0],[0]]) 
b = array([[ nan, 1.]]) 
dot(a, b) 
=> array([[ 0., 0.], [ 0., 0.]]) # EXPECTED : array([[ nan, 0.], [ nan, 0.]]) 
# (works as expected in 1.6.2 and 1.8.0) 

Fall # 4 scheint richtig in v1.6.2 und v1.8.0 zu arbeiten, aber nicht Fall # 2 ...

---

EDIT: @seberg wies darauf hin, das ist ein Problem blas , hier so ist die Information über die blas Installation I, indem sie from numpy.distutils.system_info import get_info; get_info('blas_opt') gefunden:

1.6.1 linux 64bit 
/usr/lib/python2.7/dist-packages/numpy/distutils/system_info.py:1423: UserWarning: 
    Atlas (http://math-atlas.sourceforge.net/) libraries not found. 
    Directories to search for the libraries can be specified in the 
    numpy/distutils/site.cfg file (section [atlas]) or by setting 
    the ATLAS environment variable. 
    warnings.warn(AtlasNotFoundError.__doc__) 
{'libraries': ['blas'], 'library_dirs': ['/usr/lib'], 'language': 'f77', 'define_macros': [('NO_ATLAS_INFO', 1)]} 

1.8.0 windows 32bit (anaconda) 
c:\Anaconda\Lib\site-packages\numpy\distutils\system_info.py:1534: UserWarning: 
    Blas (http://www.netlib.org/blas/) sources not found. 
    Directories to search for the sources can be specified in the 
    numpy/distutils/site.cfg file (section [blas_src]) or by setting 
    the BLAS_SRC environment variable. 
warnings.warn(BlasSrcNotFoundError.__doc__) 
{} 

(ich persönlich weiß nicht, was daraus zu machen)

Answer 1

Ich denke, wie Seberg vorgeschlagen hat, ist dies ein Problem mit der verwendeten BLAS-Bibliothek. Wenn Sie sich ansehen, wie numpy.dot implementiert ist here und here, finden Sie einen Aufruf von cblas_dgem() für den doppelt-genauen Matrix-Zeit-Matrix-Fall.

Dieses C-Programm, das einige Ihrer Beispiele reproduziert, gibt die gleiche Ausgabe, wenn Sie "einfaches" BLAS verwenden, und die richtige Antwort, wenn Sie ATLAS verwenden.

#include <stdio.h> 
#include <math.h> 

#include "cblas.h" 

void onebyone(double a11, double b11, double expectc11) 
{ 
    enum CBLAS_ORDER order=CblasRowMajor; 
    enum CBLAS_TRANSPOSE transA=CblasNoTrans; 
    enum CBLAS_TRANSPOSE transB=CblasNoTrans; 
    int M=1; 
    int N=1; 
    int K=1; 
    double alpha=1.0; 
    double A[1]={a11}; 
    int lda=1; 
    double B[1]={b11}; 
    int ldb=1; 
    double beta=0.0; 
    double C[1]; 
    int ldc=1; 

    cblas_dgemm(order, transA, transB, 
       M, N, K, 
       alpha,A,lda, 
       B, ldb, 
       beta, C, ldc); 

    printf("dot([ %.18g],[%.18g]) -> [%.18g]; expected [%.18g]\n",a11,b11,C[0],expectc11); 
} 

void onebytwo(double a11, double b11, double b12, 
       double expectc11, double expectc12) 
{ 
    enum CBLAS_ORDER order=CblasRowMajor; 
    enum CBLAS_TRANSPOSE transA=CblasNoTrans; 
    enum CBLAS_TRANSPOSE transB=CblasNoTrans; 
    int M=1; 
    int N=2; 
    int K=1; 
    double alpha=1.0; 
    double A[]={a11}; 
    int lda=1; 
    double B[2]={b11,b12}; 
    int ldb=2; 
    double beta=0.0; 
    double C[2]; 
    int ldc=2; 

    cblas_dgemm(order, transA, transB, 
       M, N, K, 
       alpha,A,lda, 
       B, ldb, 
       beta, C, ldc); 

    printf("dot([ %.18g],[%.18g, %.18g]) -> [%.18g, %.18g]; expected [%.18g, %.18g]\n", 
     a11,b11,b12,C[0],C[1],expectc11,expectc12); 
} 

int 
main() 
{ 
    onebyone(0, 0, 0); 
    onebyone(2, 3, 6); 
    onebyone(NAN, 0, NAN); 
    onebyone(0, NAN, NAN); 
    onebytwo(0, 0,0, 0,0); 
    onebytwo(2, 3,5, 6,10); 
    onebytwo(0, NAN,0, NAN,0); 
    onebytwo(NAN, 0,0, NAN,NAN); 
    return 0; 
} 

Ausgang mit BLAS:

dot([ 0],[0]) -> [0]; expected [0] 
dot([ 2],[3]) -> [6]; expected [6] 
dot([ nan],[0]) -> [nan]; expected [nan] 
dot([ 0],[nan]) -> [0]; expected [nan] 
dot([ 0],[0, 0]) -> [0, 0]; expected [0, 0] 
dot([ 2],[3, 5]) -> [6, 10]; expected [6, 10] 
dot([ 0],[nan, 0]) -> [0, 0]; expected [nan, 0] 
dot([ nan],[0, 0]) -> [nan, nan]; expected [nan, nan] 

Ausgabe mit ATLAS:

dot([ 0],[0]) -> [0]; expected [0] 
dot([ 2],[3]) -> [6]; expected [6] 
dot([ nan],[0]) -> [nan]; expected [nan] 
dot([ 0],[nan]) -> [nan]; expected [nan] 
dot([ 0],[0, 0]) -> [0, 0]; expected [0, 0] 
dot([ 2],[3, 5]) -> [6, 10]; expected [6, 10] 
dot([ 0],[nan, 0]) -> [nan, 0]; expected [nan, 0] 
dot([ nan],[0, 0]) -> [nan, nan]; expected [nan, nan] 

BLAS scheint Verhalten zu erwarten hat, wenn der erste Operand eine NaN hat, und das Unrecht, wenn der erste Operand ist Null und das zweite hat ein NaN.

Jedenfalls glaube ich nicht, dass dieser Fehler in der Numpy-Ebene ist; Es ist in BLAS. Es scheint möglich zu sein, Problemumgehung zu verwenden, indem Sie stattdessen ATLAS verwenden.

Oben auf Ubuntu 14.04 generiert, mit Ubuntu bereitgestellten gcc, BLAS und ATLAS.