Fastest way to compute a multiplication of matrices times a sequence of kronecker products

I would like to compute the following operation

o=f*(kron(a0,a1,a1,a1,a1)+kron(a1,a0,a1,a1,a1)+kron(a1,a1,a0,a1,a1)+kron(a1,a1,a1,a0,a1)+kron(a1,a1,a1,a1,a0))*B

All matrices involved in the operation are very sparse. In particular,

a0 is of size 38 x 6 with 141 nonzero elements

a1 of size 38 x 6 with 32 nonzero elements

f is of size 25 x 79,235,168 with 4698 nonzero elements

B is of size 7776 x 7776 with 18 nonzero elements

In the notation above kron(a0,a1,a1,a1,a1) is a shorthand for kron(kron(kron(kron(a0,a1),a1),a1),a1)

How fast can this problem be solved and what is the fastest way to solve it?

[Matt J]

You are missing a closing ')' somewhere in your equation. I assume it is supposed to be this?

o=f*(  kron(__)+kron(___)+kron(___)+kron(___)    )*B

[Patrick Mboma]

Yes Matt,

Thank you

Answer 1

[Matt J]

With this FEX package

https://www.mathworks.com/matlabcentral/fileexchange/25969-efficient-object-oriented-kronecker-product-manipulation

it takes about 5 sec. on my machine: (EDIT: now it's much faster if we aggregate some of the operands into explicit Kronecker products, about 0.02 sec.)

a0=sprand(38,6,141/38/6);
a1=sprand(38,6,38/38/6);
ft=sprand(79235168,25,4698/25/79235168); %f tranposed - less memory
B=sprand(7776,7776,18/7776^2);

tic
K{1}=KronProd({kron(a0,kron(a1,a1)),kron(a1,a1)}).';
K{2}=KronProd({kron(a1,kron(a0,a1)),kron(a1,a1)}).';
K{3}=KronProd({kron(a1,kron(a1,a0)),kron(a1,a1)}).';
K{4}=KronProd({kron(a1,kron(a1,a1)),kron(a0,a1)}).';
K{5}=KronProd({kron(a1,kron(a1,a1)),kron(a1,a0)}).';
toc%Elapsed time is 0.007476 seconds.

tic;    
o=0;
for i=1:5
    o=o+K{i}*ft;
end
o=o.'*B;
toc%Elapsed time is 0.015946 seconds.

My code works with the transpose of f, since this consumes a lot less memory. Compare:

f=ft.';
whos f ft

  Name             Size                      Bytes  Class     Attributes

  f               25x79235168            633956520  double    sparse    
  ft        79235168x25                      75376  double    sparse    

[Patrick Mboma]

Essentially you would recommend computing and storing the kronecker products before multiplying to f

I will try and let you know how well it works.

Thanks

[Matt J]

Essentially you would recommend computing and storing the kronecker products before multiplying to f

No, not at all! The KronProd object avoids explicit computation of the product matrix.

[Patrick Mboma]

Hi Matt,

I got the following error

Error using procShape

Cannot determine reshaping rule

Error in * (line 96)

[wasrowshaped,layers]=procShape(M.rangesizes,size(X),1);

[Matt J]

My bad. I didn't notice you had 5 kron operands, rather than 4.

[Patrick Mboma]

OK. I think I can fix that

[Matt J]

I've adjusted my posted demo. It takes about 5 sec. on my machine.

[Patrick Mboma]

I was hoping for less than one second as I will be running this thousands of times and it is just one operation out of many others but this is already an improvement.

Thanks a lot Matt

[Matt J]

Using a parfor loop will cut the time down somewhat - for me to about 2.5 sec.

[Patrick Mboma]

One last question. Since K1, K2, K3, K4 and K5 are all permutations of the same tensor, would it be possible to generate K2, K3, K4, K5 out of K1?

It that were possible, there wouldn't be a need to generate K2, K3, K4 and K5 and in that case the calculations would simplify to something like

 o=0;
    for i=1:5
        [Li,Ri]=generate_rotation(orderi)
      o=o+Bt*(Li*K1*Ri*ft);
    end

[Matt J]

The creation of the K{i} are a negligible part of the computation. It won't save you anything.

But a simpler way to permute the order would be,

    for i=1:5
       K1.opinds=orderi;
      o=o+Bt*(K1*ft);
    end

[Matt J]

I find it to be quite a bit faster if you re-express the KronProd objects in terms of smaller, explicit Kronecker products, even without parfor:

a0=sprand(38,6,141/38/6);
a1=sprand(38,6,38/38/6);
ft=sprand(79235168,25,4698/25/79235168);
B=sprand(7776,7776,18/7776^2);

tic
    K{1}=KronProd({kron(a0,kron(a1,a1)),kron(a1,a1)}).';
    K{2}=KronProd({kron(a1,kron(a0,a1)),kron(a1,a1)}).';
    K{3}=KronProd({kron(a1,kron(a1,a0)),kron(a1,a1)}).';
    K{4}=KronProd({kron(a1,kron(a1,a1)),kron(a0,a1)}).';
    K{5}=KronProd({kron(a1,kron(a1,a1)),kron(a1,a0)}).';
toc%Elapsed time is 0.007476 seconds.


tic;    

    o=0;
    for i=1:5
      o=o+K{i}*ft;
    end
    o=o.'*B;

toc%Elapsed time is 0.015946 seconds.

[Patrick Mboma]

This behavior is strange. I just replaced the first part of the block with the earlier syntax. The second part becomes awfully slow ... and the results are different too

tic
    K{1}=KronProd({a0,a1},[1,2,2,2,2]).';
    K{2}=KronProd({a0,a1},[2,1,2,2,2]).';
    K{3}=KronProd({a0,a1},[2,2,1,2,2]).';
    K{4}=KronProd({a0,a1},[2,2,2,1,2]).';
    K{5}=KronProd({a0,a1},[2,2,2,2,1]).';
toc%Elapsed time is 0.005760 seconds.
tic;    
    o=0;
    for i=1:5
      o=o+K{i}*ft;
    end
    o=o.'*B;
toc%Elapsed time is 5.049683 

[Matt J]

Did you see my last post? I have a faster recommendation.

[Patrick Mboma]

Hi Matt,

yes I saw your last recommendation, thanks a lot.

I do not understand, however, why the results are different when I use the old setup. That is not just in the computation time, but also in the actual numbers produced.

[Matt J]

Which two versions are being compared? I've lost track.

Note however, that KronProd({A,B,C}) is the same as kron(C,kron(B,A)).

[Matt J]

I did the comparison below. The computation time differs substantially (as expected), but I do not see a significant numerical difference in the result:

a0=sprand(38,6,141/38/6);
a1=sprand(38,6,38/38/6);
ft=sprand(79235168,25,4698/25/79235168);
B=sprand(7776,7776,18/7776^2);

tic
K{1}=KronProd({kron(a0,kron(a1,a1)),kron(a1,a1)}).';
K{2}=KronProd({kron(a1,kron(a0,a1)),kron(a1,a1)}).';
K{3}=KronProd({kron(a1,kron(a1,a0)),kron(a1,a1)}).';
K{4}=KronProd({kron(a1,kron(a1,a1)),kron(a0,a1)}).';
K{5}=KronProd({kron(a1,kron(a1,a1)),kron(a1,a0)}).';
toc


    o=0;
    for i=1:5
      o=o+K{i}*ft;
    end
    o=o.'*B;

o1=o;


    K{1}=KronProd({a0,a1},[1,2,2,2,2]).';
    K{2}=KronProd({a0,a1},[2,1,2,2,2]).';
    K{3}=KronProd({a0,a1},[2,2,1,2,2]).';
    K{4}=KronProd({a0,a1},[2,2,2,1,2]).';
    K{5}=KronProd({a0,a1},[2,2,2,2,1]).';

  
    o=0;
    for i=1:5
      o=o+K{i}*ft;
    end
    o=o.'*B;

o2=o;

PercentDiscrepancy = max(abs(o2(:)-o1(:)))/max(abs(o2(:)))*100   %2.2200e-14 

[Patrick Mboma]

Sorry, my mistake! I have rechecked and you are right, the results are numerically the same.

Thank you so much !