Hadamard product?

[personalityson]

Elementwise vector product (like dot, but without the reduction to one number)
Does it exist?

In MKL its vdMul, in cuBLAS its DHAD (so, not part of the standard BLAS function list)...

Hadamard product is completely central in machine learning, and is the only reason which prevents me from using OpenBLAS...

[martin-frbg]

Does not exist currently but sounds interesting. Would you happen to have some rough number for the speedup from using MKL compared to just letting a modern compiler optimize a naive implementation ?

[brada4]

It used to be in BLAS.... Still in lapack 3.1.1
it will be 3/2 execution time of DOT (3 memory accesses in place of 2 for dot per multiplication), and _dot is only C
Probably C macro will perform much better than a function call...

[martin-frbg]

@brada4 the only trace of this function appears to be a spurious entry in the intro_blas1 manpage distributed with lapack 3.1.1

[personalityson]

Just learned there is a walkaround by using sbmv
http://stackoverflow.com/questions/7621520/element-wise-vector-vector-multiplication-in-blas

Elementwise matrix*matrix product would still be nice...

[martin-frbg]

Not sure abusing sbmv for this would be more efficient than a simple loop over all elements ? BTW the only references for *HAD functions in cuBLAS appear to be people asking if there are any...

[brada4]

Now imagine your sbmv multiplying 2 million-element vectors ...
There is less wasteful way with _gemv or _gemm (i.e free dimension(s) ==1) (by magnitude slower than a loop)
dot is simple loop. I say macro because for small samples library call would be enormous time waste.

[personalityson]

less wasteful way with _gemv or _gemm (i.e free dimension(s) ==1)

Please explain... Do you extract the diagonal from the resulting matrix somehow? How is it less wasteful?

[brada4]

Diagonal? It is not square.

[martin-frbg]

Not sure I get the gemm/gemv idea either, but it could be instructive to see a quick benchmark comparison of that sbmv suggestion from stackoverflow against a simple loop that the compiler can unroll as needed.
As far as I could find out, the *HAD functions were a brainchild of Cray/SGI in the mid-1990s (proposed for future inclusion in LAPACK at some netlib developer meeting in 1995) and part of their SCSF library but never universally supported.

[personalityson]

The only way I see with gemv/gemm is to multiply two vectors v1T * v2 (of the same size), which returns a square matrix, and the diagonal is actually the hadamard product.
I spent the whole weekend looking if there were some tricks with lda/incx, but I'm too dumb to find anything. Maybe there is? :)

sbmv allows the input matrix to be stored in banded matrix form. With number of super-diagonals = 0 you only need to supply the diagonal, which in our case is one of the vectors, so it's not so bad...

I should explain the context to my problem:
I'm on a work computer where I do not have admin rights, and cannot install or compile anything.
I use Excel VBA to declare function calls to libopenblas.dll in order to do some matrix crunching.
(When I write that hadamard is "the only reason which prevents me from using OpenBLAS" -- I'm bluffing. mkl_rt.dll does not allow stdcall declarations from VBA, or maybe there is another reason, in any case I can only use OpenBLAS.)
So, my two options are sbmv or a C Macro VBA macro for hadamard, and ideally I would want to have it outside of VBA...

[brada4]

swap dimensions of v2 and get HAD
swap dimensions of v1 and get DOT

[grisuthedragon]

The Hadamard product is obviously a memory bounded operation, so the only thing one had to take care of is that the matrices are not traversed orthogonal to their storage scheme. If the matrices are accesses in the right way one has 3 loads, 1 multiplication and one store. The only point I see where one can optimize something is if we assume $ C := \alpha A\circle B + \beta C $ with $\beta=0$ and using non-temporal stores for C. In this case we only have 2 loads and one store per multiplication but requires aligned access to C and C must be larger than the L3 cache of current Intel-CPUs to give good results.

[brada4]

in common case one can treat matrices as 1:(M*N) vectors and apply marginal case of gemm / gemv

[personalityson]

I still don't see it. Please give us a concrete example
Assume size(v1) = size(v2)