Strange matrix multiplication slowdown

edit: It looks like it's a problem with floating point underflow. Many values in the W matrix are very close to zero and are represented as denormal numbers. According to Wikipedia:

"But even when denormal values are entirely computed in hardware, the speed of computation is significantly reduced on most modern processors; in extreme cases, instructions involving denormal operands may run as much as 100 times slower."

So that explains things. Learn something new about floating point every day!

While training a neural net, I came across some odd behavior where a matrix multiplication sees a significant slowdown as training progresses. Here is a copy/paste from the command window illustrating the strange behavior:

load svhn_3_13_12_10_21_53_47
whos

Name      Size               Bytes  Class     Attributes
nn        1x1             33075406  struct              

A = rand(2500, size(nn.W{1},2), 'single');
whos

Name         Size                 Bytes  Class     Attributes

A         2500x1025            10250000  single              
nn           1x1               33075406  struct              

tic; 
A*nn.W{1}';
toc

Elapsed time is 5.353390 seconds.

tic;
A*rand(size(nn.W{1}), 'single')';
toc

 Elapsed time is 0.130106 seconds.

nn.W{1} is a single-precision matrix. Multiplying A by nn.W{1}' takes 5.3 seconds, but multiplying A by the same sized matrix takes 0.13 seconds. I can think of no reasonable explanation why any two matrix-matrix multiplication of the same size can differ so much in execution time.

Curiously, if I recast both A and nn.W{1} as a double, I can reclaim significant performance. If I only recast one of them, I am back to the same lousy performance:

tic; 
double(A) * double(nn.W{1}'); 
toc
Elapsed time is 0.175379 seconds.

tic; double(A) * W'; toc
Elapsed time is 0.794721 seconds.

tic;
A * double(nn.W{1}'); 
toc
Elapsed time is 5.670088 seconds.

Ideas?

[Cedric]

I would be interested to get an answer to this question as well, possibly from Mathworks(?) I think that your investigations hit the point..

 >> nnz(W)
 ans =
     1537499

 >> V = abs(W) ;
 >> v = V(V>0) ;
 >> nnz(v < realmin('single'))     % Most n-z elements are denormal singles.
 ans =
     1138506

 >> nnz(v < realmin('double'))
 ans =
     0

and obviously

 >> nnz(rand(size(W)) < realmin('double'))
 ans =
     0

which is also the case with singles.

Now it's difficult to say more when we don't know how MATLAB treats denormal floating-points internally.

[Matt J]

Makes one wonder whether simple scaling, to undo the denormalization, would be an alternative to casting

 s=2^30;

 tic
 (A*(nn.W{1}*s))/s
 toc