I think, as seberg suggested, this is an issue with the BLAS library used. If you look at how numpy.dot is implemented here and here you'll find a call to cblas_dgemm() for the double-precision matrix-times-matrix case.
This C program, which reproduces some of your examples, gives the same output when using "plain" BLAS, and the right answer when using ATLAS.
#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;
}
Output with 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]
Output with 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 seems to have expected behaviour when the first operand has a NaN, and the wrong when the first operand is zero and the second has a NaN.
Anyway, I don't think this bug is in the Numpy layer; it's in BLAS. It appears to be possible to workaround by using ATLAS instead.
Above generated on Ubuntu 14.04, using Ubuntu-provided gcc, BLAS, and ATLAS.