I know `@` is for decorators, but what is `@=` for in Python? Is it just reservation for some future idea?

This is just one of my many questions while reading `tokenizer.py`.

What is the &#39;@=&#39; symbol for in Python?

I am making some benchmarks with CUDA, C++, C#, Java, and using MATLAB for verification and matrix generation.  When I perform matrix multiplication with MATLAB, `2048x2048` and even bigger matrices are almost instantly multiplied.

&lt;!-- language: none --&gt;

                 1024x1024   2048x2048   4096x4096
                 ---------   ---------   ---------
    CUDA C (ms)      43.11      391.05     3407.99
    C++ (ms)       6137.10    64369.29   551390.93
    C# (ms)       10509.00   300684.00  2527250.00
    Java (ms)      9149.90    92562.28   838357.94
    MATLAB (ms)      75.01      423.10     3133.90

Only CUDA is competitive, but I thought that at least C++ will be somewhat close and not 60 times slower. I also don&#39;t know what to think about the C# results. The algorithm is just the same as C++ and Java, but there&#39;s a giant jump `2048` from `1024`.

How is MATLAB performing matrix multiplication so fast? 

C++ Code:

```cpp
float temp = 0;
timer.start();
for(int j = 0; j &lt; rozmer; j++)
{
	for (int k = 0; k &lt; rozmer; k++)
    {
        temp = 0;
        for (int m = 0; m &lt; rozmer; m++)
        {
			temp = temp + matice1[j][m] * matice2[m][k];
        }
        matice3[j][k] = temp;
    }
}
timer.stop();
```

Why is MATLAB so fast in matrix multiplication?

I recently moved to Python 3.5 and noticed the [new matrix multiplication operator (@)][1] sometimes behaves differently from the [numpy dot][2] operator. In example, for 3d arrays:

    import numpy as np
    
    a = np.random.rand(8,13,13)
    b = np.random.rand(8,13,13)
    c = a @ b  # Python 3.5+
    d = np.dot(a, b)

The `@` operator returns an array of shape:

    c.shape
    (8, 13, 13)

while the `np.dot()` function returns:

    d.shape
    (8, 13, 8, 13)

How can I reproduce the same result with numpy dot? Are there any other significant differences?


  [1]: https://docs.python.org/3/whatsnew/3.5.html#whatsnew-pep-465
  [2]: http://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html

Difference between numpy dot() and Python 3.5+ matrix multiplication @

I have two matrices 

    a = np.matrix([[1,2], [3,4]])
    b = np.matrix([[5,6], [7,8]])

and I want to get the element-wise product, `[[1*5,2*6], [3*7,4*8]]`, equaling 

`[[5,12], [21,32]]`

I have tried

    print(np.dot(a,b)) 

and

    print(a*b)

but both give the result

`[[19 22], [43 50]]`

which is the matrix product, not the element-wise product. How can I get the the element-wise product (aka Hadamard product) using built-in functions?

How to get element-wise matrix multiplication (Hadamard product) in numpy?

The numpy docs recommend using array instead of matrix for working with matrices. However, unlike octave (which I was using till recently), * doesn&#39;t perform matrix multiplication, you need to use the function matrixmultipy(). I feel this makes the code very unreadable.

Does anybody share my views, and has found a solution?

how does multiplication differ for NumPy Matrix vs Array classes?

I am making some matrix multiplication benchmarking, as previously mentioned in 
https://stackoverflow.com/questions/6058139/why-is-matlab-so-fast-in-matrix-multiplication

Now I&#39;ve got another issue, when multiplying two 2048x2048 matrices, there is a big difference between C# and others. When I try multiply only 2047x2047 matrices, it seems normal. Added some others for comparsion too.

1024x1024 - 10 seconds.

1027x1027 - 10 seconds.

2047x2047 - 90 seconds.

2048x2048 - 300 seconds.

2049x2049 - 91 seconds. (update)

2500x2500 - 166 seconds

That is three and a half minute difference for the 2k by 2k case.

using 2dim arrays


    //Array init like this
    int rozmer = 2048;
    float[,] matice = new float[rozmer, rozmer];

    //Main multiply code
    for(int j = 0; j &lt; rozmer; j++)
    {
       for (int k = 0; k &lt; rozmer; k++)
       {
         float temp = 0;
         for (int m = 0; m &lt; rozmer; m++)
         {
           temp = temp + matice1[j,m] * matice2[m,k];
         }
         matice3[j, k] = temp;
       }
     }

Why is there huge performance hit in 2048x2048 versus 2047x2047 array multiplication?

I have a matrix multiply code that looks like this:

    for(i = 0; i &lt; dimension; i++)
        for(j = 0; j &lt; dimension; j++)
            for(k = 0; k &lt; dimension; k++)
                C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];

Here, the size of the  matrix is represented by `dimension`.
Now, if the size of the matrices is 2000, it takes 147 seconds to run this piece of code, whereas if the size of the matrices is 2048, it takes 447 seconds. So while the difference in no. of multiplications is (2048*2048*2048)/(2000*2000*2000) = 1.073, the difference in the timings is 447/147 = 3. Can someone explain why this happens? I expected it to scale linearly, which does not happen. I am not trying to make the fastest matrix multiply code, simply trying to understand why it happens. 

Specs: AMD Opteron dual core node (2.2GHz), 2G RAM, gcc v 4.5.0

Program compiled as `gcc -O3 simple.c`

I have run this on Intel&#39;s icc compiler as well, and seen similar results.

EDIT:

As suggested in the comments/answers, I ran the code with dimension=2060 and it takes 145 seconds.

Heres the complete program:

    #include &lt;stdlib.h&gt;
    #include &lt;stdio.h&gt;
    #include &lt;sys/time.h&gt;
    
    /* change dimension size as needed */
    const int dimension = 2048;
    struct timeval tv; 
    
    double timestamp()
    {
            double t;
            gettimeofday(&amp;tv, NULL);
            t = tv.tv_sec + (tv.tv_usec/1000000.0);
            return t;
    }
    
    int main(int argc, char *argv[])
    {
            int i, j, k;
            double *A, *B, *C, start, end;
    
            A = (double*)malloc(dimension*dimension*sizeof(double));
            B = (double*)malloc(dimension*dimension*sizeof(double));
            C = (double*)malloc(dimension*dimension*sizeof(double));
    
            srand(292);
    
            for(i = 0; i &lt; dimension; i++)
                    for(j = 0; j &lt; dimension; j++)
                    {   
                            A[dimension*i+j] = (rand()/(RAND_MAX + 1.0));
                            B[dimension*i+j] = (rand()/(RAND_MAX + 1.0));
                            C[dimension*i+j] = 0.0;
                    }   
    
            start = timestamp();
            for(i = 0; i &lt; dimension; i++)
                    for(j = 0; j &lt; dimension; j++)
                            for(k = 0; k &lt; dimension; k++)
                                    C[dimension*i+j] += A[dimension*i+k] *
                                            B[dimension*k+j];
    
            end = timestamp();
            printf(&quot;\nsecs:%f\n&quot;, end-start);
    
            free(A);
            free(B);
            free(C);
    
            return 0;
    }

Matrix multiplication: Small difference in matrix size, large difference in timings

I&#39;m new to the CUDA paradigm. My question is in determining the number of threads per block, and blocks per grid. Does a bit of art and trial play into this? What I&#39;ve found is that many examples have seemingly arbitrary number chosen for these things.

I&#39;m considering a problem where I would be able to pass matrices - of any size - to a method for multiplication. So that, each element of C (as in C = A * B) would be calculated by a single thread. How would you determine the threads/block, blocks/grid in this case?

CUDA determining threads per block, blocks per grid

I&#39;ve stumbled upon the weird way (in my view) that Matlab is dealing with [empty matrices][1]. For example, if two empty matrices are multiplied the result is:

    zeros(3,0)*zeros(0,3)
    ans =

     0     0     0
     0     0     0
     0     0     0

Now, this already took me by surprise, however, a quick search got me to the link above, and I got an explanation of the somewhat twisted logic of why this is happening. 

**However**, nothing prepared me for the following observation. I asked myself, how efficient is this type of multiplication vs just using `zeros(n)` function, say for the purpose of initialization?  I&#39;ve used [`timeit`][2] to answer this:

    f=@() zeros(1000)
    timeit(f)
    ans =
        0.0033

vs:

    g=@() zeros(1000,0)*zeros(0,1000)
    timeit(g)
    ans =
        9.2048e-06

Both have the same outcome of 1000x1000 matrix of zeros of class `double`, but the empty matrix multiplication one is ~350 times faster! (a similar result happens using `tic` and `toc` and a loop)

How can this be? are `timeit` or `tic,toc` bluffing or have I found a faster way to initialize matrices? 
(this was done with matlab 2012a, on a win7-64 machine,  intel-i5 650 3.2Ghz...)

**EDIT:**

After reading your feedback, I have looked more carefully into this peculiarity, and tested on 2 different computers (same matlab ver though 2012a) a code that examine the run time vs the size of matrix n. This is what I get:

![enter image description here][3]

The code to generate this used `timeit` as before, but a loop with `tic` and `toc` will look the same. So, for small sizes, `zeros(n)` is comparable. However, around `n=400` there is a jump in performance for the empty matrix multiplication. The code I&#39;ve used to generate that plot was:

    n=unique(round(logspace(0,4,200)));
    for k=1:length(n)
        f=@() zeros(n(k));
        t1(k)=timeit(f);

        g=@() zeros(n(k),0)*zeros(0,n(k));
        t2(k)=timeit(g);
    end

    loglog(n,t1,&#39;b&#39;,n,t2,&#39;r&#39;);
    legend(&#39;zeros(n)&#39;,&#39;zeros(n,0)*zeros(0,n)&#39;,2);
    xlabel(&#39;matrix size (n)&#39;); ylabel(&#39;time [sec]&#39;);

Are any of you experience this too?

**EDIT #2:**

Incidentally, empty matrix multiplication is not needed to get this effect. One can simply do:

    z(n,n)=0;

where n&gt; some threshold matrix size seen in the previous graph,  and get the **exact** efficiency profile as with empty matrix multiplication (again using timeit). 

![enter image description here][4]

Here&#39;s an example where it improves efficiency of a code:

    n = 1e4;
    clear z1
    tic
    z1 = zeros( n ); 
    for cc = 1 : n
        z1(:,cc)=cc;
    end
    toc % Elapsed time is 0.445780 seconds.

    %%
    clear z0
    tic
    z0 = zeros(n,0)*zeros(0,n);
    for cc = 1 : n
        z0(:,cc)=cc;
    end
    toc % Elapsed time is 0.297953 seconds.

However, using `z(n,n)=0;` instead yields similar results to the `zeros(n)` case.

  [1]: http://www.mathworks.com/help/matlab/math/empty-matrices-scalars-and-vectors.html
  [2]: http://www.mathworks.com/matlabcentral/fileexchange/18798
  [3]: http://i.stack.imgur.com/DSZn7.png
  [4]: http://i.stack.imgur.com/EZVOI.png

Faster way to initialize arrays via empty matrix multiplication? (Matlab)

I want to calculate the row-wise dot product of two matrices of the same dimension as fast as possible. This is the way I am doing it:

&lt;!--language: python--&gt;

    import numpy as np
    a = np.array([[1,2,3], [3,4,5]])
    b = np.array([[1,2,3], [1,2,3]])
    result = np.array([])
    for row1, row2 in a, b:
        result = np.append(result, np.dot(row1, row2))
    print result

and of course the output is:

    [ 26.  14.]

Vectorized way of calculating row-wise dot product two matrices with Scipy

I was trying to figure out the fastest way to do matrix multiplication and tried 3 different ways:

- Pure python implementation: no surprises here.
- Numpy implementation using `numpy.dot(a, b)`
- Interfacing with C using `ctypes` module in Python.

This is the C code that is transformed into a shared library:

    #include &lt;stdio.h&gt;
    #include &lt;stdlib.h&gt;

    void matmult(float* a, float* b, float* c, int n) {
    	int i = 0;
    	int j = 0;
    	int k = 0;

    	/*float* c = malloc(nay * sizeof(float));*/
    
    	for (i = 0; i &lt; n; i++) {
    		for (j = 0; j &lt; n; j++) {
    			int sub = 0;
    			for (k = 0; k &lt; n; k++) {
    				sub = sub + a[i * n + k] * b[k * n + j];
    			}
    			c[i * n + j] = sub;
    		}
    	}
    	return ;
    }

And the Python code that calls it:

    def C_mat_mult(a, b):
        libmatmult = ctypes.CDLL(&quot;./matmult.so&quot;)

        dima = len(a) * len(a)
        dimb = len(b) * len(b)

        array_a = ctypes.c_float * dima
        array_b = ctypes.c_float * dimb
        array_c = ctypes.c_float * dima

        suma = array_a()
        sumb = array_b()
        sumc = array_c()

        inda = 0
        for i in range(0, len(a)):
            for j in range(0, len(a[i])):
                suma[inda] = a[i][j]
                inda = inda + 1
            indb = 0
        for i in range(0, len(b)):
            for j in range(0, len(b[i])):
                sumb[indb] = b[i][j]
                indb = indb + 1

        libmatmult.matmult(ctypes.byref(suma), ctypes.byref(sumb), ctypes.byref(sumc), 2);

        res = numpy.zeros([len(a), len(a)])
        indc = 0
        for i in range(0, len(sumc)):
            res[indc][i % len(a)] = sumc[i]
            if i % len(a) == len(a) - 1:
                indc = indc + 1
        
        return res

I would have bet that the version using C would have been faster ... and I&#39;d have lost ! Below is my benchmark which seems to show that I either did it incorrectly, or that `numpy` is stupidly fast:

![benchmark][1]

I&#39;d like to understand why the `numpy` version is faster than the `ctypes` version, I&#39;m not even talking about the pure Python implementation since it is kind of obvious.

  [1]: http://i.stack.imgur.com/TLhDi.png

Why is matrix multiplication faster with numpy than with ctypes in Python?

I know that, in the 1D case, the convolution between two vectors, `a` and `b`, can be computed as `conv(a, b)`, but also as the product between the `T_a` and `b`, where `T_a` is the corresponding Toeplitz matrix for `a`.

Is it possible to extend this idea to 2D?

Given `a = [5 1 3; 1 1 2; 2 1 3]` and `b=[4 3; 1 2]`, is it possible to convert `a` in a Toeplitz matrix and compute the matrix-matrix product between `T_a` and `b` as in the 1-D case?

2-D convolution as a matrix-matrix multiplication

When working with matrices in R, one can put them side-by-side or stack them top of each other using cbind and rbind, respectively. What is the equivalent function for stacking matrices or arrays in other dimensions?

For example, the following creates a pair of 2x2 matrices, each having 4 elements:

    x = cbind(1:2,3:4)
    y = cbind(5:6,7:8)

What is the code to combine them into a 2x2x2 array with 8 elements?


All Matrix Multiplication Solutions on Gang of Coders

Total of 13 Matrix Multiplication Solutions

What is the '@=' symbol for in Python?