## Do we have an efficient way to do this matrix multiplication?

### Sinwoo Jeong (view profile)

on 22 Jan 2019
Latest activity Commented on by Steven Lord

on 22 Jan 2019

### Jan (view profile)

What I wanna do is just simple matrix multiplication of two matrix.
I hope I can do this in a more efficient way.
Here we have the MATLAB code for doing this calculations.
a = 10;
b = 8;
n = 500;
A_MAT = zeros(a*n, b*n);
B_MAT = zeros(b*n, a*n);
for j = 1:100
for i = 1:n
A_MAT(1+(i-1)*a:i*a, 1+(i-1)*b:i*b) = rand(a,b);
B_MAT(1+(i-1)*b:i*b, 1+(i-1)*a:i*a) = rand(b,a);
end
A_MAT = sparse(A_MAT);
B_MAT = sparse(B_MAT);
% 1st way (naive MATLAB * operation)
tic;
C_MAT = A_MAT*B_MAT;
time_1(j) = toc;
% 2nd way (for-loop operation)
tic;
for i = 1:n
C_MAT_2(1+(i-1)*a:i*a, 1+(i-1)*a:i*a) = A_MAT(1+(i-1)*a:i*a, 1+(i-1)*b:i*b)*B_MAT(1+(i-1)*b:i*b, 1+(i-1)*a:i*a);
end
time_2(j) = toc;
time_ratio(j) = time_2(j)/time_1(j);
end
plot(time_ratio)
Here, a, b, and n can have different values (the values depend on problems I want to solve).
You'll notice that the 1st way (naive MATLAB * operation) is much faster than for-loop after running this code. Therefore, I concluded that for-loop is not an answer for my question.
I want to find more efficient way than naive MATLAB * operation (I know MATLAB * operation is already quite optimized).
The reason why I think there may be more efficient way is that matrix "A_MAT" and "B_MAT" are not only quite sparse but also very regularly organized. (You can see how well they are organized by using "spy(A_MAT)" or "spy(B_MAT)" in MATLAB).
Therefore, I thought there might be a special algoritm suitable for doing this matrix multiplication.

on 22 Jan 2019
Edited by Jan

### Jan (view profile)

on 22 Jan 2019

Matlab calls optimized functions for sparse matrix multiplications already. Sometime we had an acceleration with hard-coded sparse matrix multiplications written in Fortran77, but I assume, that the time to implement this will exceed the time you win during running. Maybe you can call one of the established libraries through a Mex function, see Armadillo, Blaze, Boost.uBlas, CSparse, Eigen, SuiteSparse .
Even if the sparsity pattern looks obvious for the human eye, it is not trivial to implement the pattern in a hard coded algorithm. Therefore the general purpose codes in the sparse libs called by Matlab are efficient already and hard to bet.

Steven Lord

### Steven Lord (view profile)

on 22 Jan 2019
I agree with what Jan said, but going a little further you may be able to gain depending on what you're trying to do.
Some of the sparse matrix functions can accept either a sparse matrix or a function that performs a computation on the sparse matrix (that you provide) and a vector (that the function provides.) One example of this is the iterative solver gmres.
See the "Using gmres with a Matrix Input" and "Using gmres with a Function Handle" examples on that page. The latter example takes advantage of the structure of the Wilkinson matrix (created using the gallery function) to compute A*x without actually needing to build A. Depending on how big A is, avoiding the need to create A explicitly (and the associated allocation of memory) may save you time or even allow you to solve problems that you otherwise could not solve.
But this is going to depend upon you manually creating the function afun and taking advantage of the specific structure of your problem. If you're looking for a general matrix multiplication that beats the one already in MATLAB -- if we were aware of one that we felt was better from a accuracy, robustness, and performance standpoint, why wouldn't we be using that one?