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

3 次查看（过去 30 天）
Patrick Mboma
Patrick Mboma 2022-8-8
评论： Patrick Mboma 2022-8-8
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
Matt J 2022-8-8
编辑：Matt J 2022-8-8
Ran in:
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
  17 个评论
显示 15更早的评论
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by