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

9 ビュー (過去 30 日間)
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?
  2 件のコメント
Matt J
Matt J 2022 年 8 月 8 日
You are missing a closing ')' somewhere in your equation. I assume it is supposed to be this?
o=f*( kron(__)+kron(___)+kron(___)+kron(___) )*B

サインインしてコメントする。

採用された回答

Matt J
Matt J 2022 年 8 月 8 日
編集済み: Matt J 2022 年 8 月 8 日
With this FEX package
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
  17 件のコメント
Matt J
Matt J 2022 年 8 月 8 日
編集済み: Matt J 2022 年 8 月 8 日
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
Patrick Mboma 2022 年 8 月 8 日
Sorry, my mistake! I have rechecked and you are right, the results are numerically the same.
Thank you so much !

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeLogical についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by