Programming
Why is MATLAB so fast in matrix multiplication
When delving into the realm of numerical computation, especially in scientific and engineering fields, one often encounters MATLAB. A common observation that frequently sparks curiosity is: why is MATLAB so fast in matrix multiplication? It’s a critical question for anyone working with large datasets or complex simulations, as efficient matrix operations are the backbone of many computational tasks, from signal processing to machine learning. MATLAB’s reputation for speed in linear algebra isn’t just anecdotal; it’s a meticulously engineered feature rooted in several core design principles and underlying technologies. This remarkable performance isn’t a simple trick, but rather a sophisticated synergy of highly optimized libraries, intelligent compilation, and efficient resource management. Understanding these mechanisms not only demystifies MATLAB’s speed but also helps users write more performant code themselves.
The Cornerstone: Optimized Linear Algebra Libraries (BLAS & LAPACK)
At the heart of MATLAB’s exceptional speed for matrix operations lies its strategic reliance on highly optimized, external libraries, primarily the Basic Linear Algebra Subprograms (BLAS) and the Linear Algebra PACKage (LAPACK). These are not proprietary MATLAB inventions; rather, they are industry-standard collections of routines for performing common vector and matrix operations. MATLAB acts as a high-level interface, calling these low-level, highly tuned functions behind the scenes whenever you perform operations like matrix multiplication, inversion, or decomposition.
The true power of BLAS and LAPACK stems from their diverse implementations. While generic versions exist, major hardware vendors like Intel (with Intel Math Kernel Library or MKL), AMD (with AMD AOCL), and even OpenBLAS provide their own highly optimized versions. These vendor-specific libraries are often hand-tuned in assembly language to exploit the unique architectural features of their processors, including instruction sets like AVX, FMA, and efficient cache utilization. This means that when you install MATLAB on a system with an Intel processor, it typically links to Intel MKL, which is optimized precisely for that hardware, giving a significant performance boost that off-the-shelf code cannot easily replicate.
For instance, a simple matrix multiplication (A B) in MATLAB doesn’t involve MATLAB interpreting each element-wise operation. Instead, it translates this into a call to a highly optimized BLAS routine (specifically, a DGEMM operation for double-precision general matrix multiplication). This delegation to highly specialized, pre-compiled code is a fundamental reason why MATLAB is so fast in matrix multiplication, allowing it to achieve near bare-metal performance for critical numerical tasks. According to Intel, using Intel MKL can provide up to 10x speedup for linear algebra routines compared to unoptimized code. You can learn more about these powerful libraries on the Netlib BLAS official page.
Just-In-Time Compilation and Vectorization
Beyond optimized external libraries, MATLAB employs sophisticated internal mechanisms to enhance performance, notably its Just-In-Time (JIT) compilation and strong emphasis on vectorization. The JIT compiler, often referred to as the “accelerator,” dynamically translates portions of MATLAB code into machine code during execution. This means that frequently executed loops or array operations can bypass the typical interpreter overhead, running much closer to the speed of compiled C or Fortran code. This is particularly effective for array-based operations, including matrix multiplication, where patterns of computation are often repetitive and predictable.
For users seeking to understand why MATLAB is so fast in matrix multiplication, it’s crucial to grasp the concept of vectorization. Instead of writing explicit loops to process array elements one by one, MATLAB encourages performing operations on entire arrays or matrices at once. For example, instead of a nested loop for matrix multiplication, you simply use the `` operator. This vectorized approach allows the JIT compiler and the underlying BLAS libraries to efficiently map these high-level operations onto optimized CPU instructions that can process multiple data elements simultaneously (SIMD - Single Instruction, Multiple Data). This significantly reduces the overhead associated with loop control and individual memory accesses.
By combining JIT compilation with vectorization, MATLAB dramatically improves cache efficiency. When data is accessed sequentially in large blocks, it’s more likely to reside in faster CPU caches (L1, L2, L3), reducing the need to fetch data from slower main memory. This data locality is a cornerstone of high-performance computing. For instance, in a vectorized matrix multiplication, the JIT compiler can generate machine code that pre-fetches data, ensuring the CPU always has data ready for processing, which directly translates to faster execution times.
Leveraging Parallel Computing and Hardware Acceleration
Modern CPUs are multi-core, and graphics processing units (GPUs) offer immense parallel processing capabilities. MATLAB is designed to automatically take advantage of these hardware Question & Answer :
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.
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’t know what to think about the C# results. The algorithm is just the same as C++ and Java, but there’s a giant jump 2048 from 1024.
How is MATLAB performing matrix multiplication so fast?
C++ Code:
float temp = 0; timer.start(); for(int j = 0; j < rozmer; j++) { for (int k = 0; k < rozmer; k++) { temp = 0; for (int m = 0; m < rozmer; m++) { temp = temp + matice1[j][m] * matice2[m][k]; } matice3[j][k] = temp; } } timer.stop();
This kind of question is recurring and should be answered more clearly than “MATLAB uses highly optimized libraries” or “MATLAB uses the MKL” for once on Stack Overflow.
History:
Matrix multiplication (together with Matrix-vector, vector-vector multiplication and many of the matrix decompositions) is (are) the most important problems in linear algebra. Engineers have been solving these problems with computers since the early days.
I’m not an expert on the history, but apparently back then, everybody just rewrote his FORTRAN version with simple loops. Some standardization then came along, with the identification of “kernels” (basic routines) that most linear algebra problems needed in order to be solved. These basic operations were then standardized in a specification called: Basic Linear Algebra Subprograms (BLAS). Engineers could then call these standard, well-tested BLAS routines in their code, making their work much easier.
BLAS:
BLAS evolved from level 1 (the first version which defined scalar-vector and vector-vector operations) to level 2 (vector-matrix operations) to level 3 (matrix-matrix operations), and provided more and more “kernels” so standardized more and more of the fundamental linear algebra operations. The original FORTRAN 77 implementations are still available on Netlib’s website.
Towards better performance:
So over the years (notably between the BLAS level 1 and level 2 releases: early 80s), hardware changed, with the advent of vector operations and cache hierarchies. These evolutions made it possible to increase the performance of the BLAS subroutines substantially. Different vendors then came along with their implementation of BLAS routines which were more and more efficient.
I don’t know all the historical implementations (I was not born or a kid back then), but two of the most notable ones came out in the early 2000s: the Intel MKL and GotoBLAS. Your Matlab uses the Intel MKL, which is a very good, optimized BLAS, and that explains the great performance you see.
Technical details on Matrix multiplication:
So why is Matlab (the MKL) so fast at dgemm (double-precision general matrix-matrix multiplication)? In simple terms: because it uses vectorization and good caching of data. In more complex terms: see the article provided by Jonathan Moore.
Basically, when you perform your multiplication in the C++ code you provided, you are not at all cache-friendly. Since I suspect you created an array of pointers to row arrays, your accesses in your inner loop to the k-th column of “matice2”: matice2[m][k] are very slow. Indeed, when you access matice2[0][k], you must get the k-th element of the array 0 of your matrix. Then in the next iteration, you must access matice2[1][k], which is the k-th element of another array (the array 1). Then in the next iteration you access yet another array, and so on… Since the entire matrix matice2 can’t fit in the highest caches (it’s 8*1024*1024 bytes large), the program must fetch the desired element from main memory, losing a lot of time.
If you just transposed the matrix, so that accesses would be in contiguous memory addresses, your code would already run much faster because now the compiler can load entire rows in the cache at the same time. Just try this modified version:
timer.start(); float temp = 0; //transpose matice2 for (int p = 0; p < rozmer; p++) { for (int q = 0; q < rozmer; q++) { tempmat[p][q] = matice2[q][p]; } } for(int j = 0; j < rozmer; j++) { for (int k = 0; k < rozmer; k++) { temp = 0; for (int m = 0; m < rozmer; m++) { temp = temp + matice1[j][m] * tempmat[k][m]; } matice3[j][k] = temp; } } timer.stop();
So you can see how just cache locality increased your code’s performance quite substantially. Now real dgemm implementations exploit that to a very extensive level: They perform the multiplication on blocks of the matrix defined by the size of the TLB (Translation lookaside buffer, long story short: what can effectively be cached), so that they stream to the processor exactly the amount of data it can process. The other aspect is vectorization, they use the processor’s vectorized instructions for optimal instruction throughput, which you can’t really do from your cross-platform C++ code.
Finally, people claiming that it’s because of Strassen’s or Coppersmith–Winograd algorithm are wrong, both these algorithms are not implementable in practice, because of hardware considerations mentioned above.
