@John Fredsted
This thread is 3 years old, I don't wish to upset anyone by "reviving" it, forgive me.
I came to this thread as I was searching for information on how to write efficient procedures.
I learned a great deal by looking at how others write a proc.
Now the LinearAlgebra package implements a KroneckerProduct so the need for user-written procedures to compute the Kronecker products must be slim. But still I got curious about how different approaches would fair with Maple14/Standard/Windows with the latest i7processor and plenty of RAM.
The built-in LinearAlgebra KroneckerProduct came top. Nice! Acer's and DJ Keenan's third proc do very well too.
As pointed out, the results vary greatly according to the size and shape of the Matrices.
I didn't do a scientific study at all, I was just fooling around. I'll copy the run times below, they can be read below the line time() - t;
--------------------------------------------------------------------------------
# Comparison of various procedures to compute the Kronecker Product of 2 Matrice
# References:
# http://www.mapleprimes.com/posts/41702-Maple-12--Wish-List
# http://www.mapleprimes.com/questions/42797-KroneckerTensor-Products
restart:
A := Matrix([[1,2,3],[4,5,6]]):
B := Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]]):
LinearAlgebra:-KroneckerProduct(A,B);
Matrix(%id = 186680652)
A := Matrix(5,5,(i,j) -> a||i||j):
B := Matrix(5,5,(i,j) -> b||i||j):
t := time():
for i from 1 to 1000 do
LinearAlgebra:-KroneckerProduct(A,B)
end do:
time() - t;
2.090
restart:
myKroneckerProduct := proc(a::Matrix,b::Matrix) # by Acer
local i,j,aRow,aCol,bRow,bCol,p,k,l;
# description "Returns the Kronecker product of A and B";
aRow,aCol := LinearAlgebra:-Dimension(a);
bRow,bCol := LinearAlgebra:-Dimension(b);
p := Matrix(aRow * bRow,aCol * bCol);
for i from 1 to aRow do
for j from 1 to aCol do
for k from 1 to bRow do
for l from 1 to bCol do
p[(i-1)*bRow+k,(j-1)*bCol+l]:=a[i,j]*b[k,l];
end do;
end do;
end do;
end do;
p
end proc:
A := Matrix([[1,2,3],[4,5,6]]):
B := Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]]):
myKroneckerProduct(A,B);
Matrix(%id = 186680460)
A := Matrix(5,5,(i,j) -> a||i||j):
B := Matrix(5,5,(i,j) -> b||i||j):
t := time():
for i from 1 to 1000 do
myKroneckerProduct(A,B)
end do:
time() - t;
2.684
restart:
myKroneckerProduct := proc(a::Matrix,b::Matrix) # by Jon Fredsted
local i,j,aRow,aCol,bRow,bCol,p;
aRow,aCol := LinearAlgebra:-Dimension(a);
bRow,bCol := LinearAlgebra:-Dimension(b);
p := Matrix(aRow * bRow,aCol * bCol);
for i from 1 to aRow do
for j from 1 to aCol do
p[ (i-1)*bRow+1..i*bRow, (j-1)*bCol+1..j*bCol ] := a[i,j]*b
end do;
end do;
p
end proc:
A := Matrix([[1,2,3],[4,5,6]]):
B := Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]]):
myKroneckerProduct(A,B);
Matrix(%id = 204764300)
A := Matrix(5,5,(i,j) -> a||i||j):
B := Matrix(5,5,(i,j) -> b||i||j):
t := time():
for i from 1 to 1000 do
myKroneckerProduct(A,B)
end do:
time() - t;
14.321
restart:
myKroneckerProduct := proc(A::Matrix, B::Matrix) # by DJ Keenan v.1
local a, b, C, Catenate, i;
uses LinearAlgebra;
# description "Returns the Kronecker product of A and B";
Catenate:= (d, A, B)-> `if`(LinearAlgebra:-Dimension(A)[d]=0, B, ArrayTools:-Concatenate(d, A, B));
a := Dimension(A); b:= Dimension(B);
C := map(curry(ScalarMultiply,B), A);
foldl(curry(Catenate,1), Matrix(0,a[2]*b[2]), seq(foldl(curry(Catenate,2), Matrix(b[1],0), op(convert(C[i,1..-1], list))), i=1..a[1]))
end proc:
A := Matrix([[1,2,3],[4,5,6]]):
B := Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]]):
myKroneckerProduct(A,B);
Matrix(%id = 188661516)
A := Matrix(5,5,(i,j) -> a||i||j):
B := Matrix(5,5,(i,j) -> b||i||j):
t := time():
for i from 1 to 1000 do
myKroneckerProduct(A,B)
end do:
time() - t;
20.592
restart:
myKroneckerProduct := proc(A::Matrix, B::Matrix) # by DJ Keenan v.1.b
local C, i;
uses ArrayTools,LinearAlgebra;
# description "Returns the Kronecker product of A and B";
C:= map(curry(ScalarMultiply, B), A);
Concatenate(1, seq(Concatenate(2, op(convert(C[i,1..-1], list))), i=1..Dimension(A)[1]))
end proc:
A := Matrix([[1,2,3],[4,5,6]]):
B := Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]]):
myKroneckerProduct(A,B);
Matrix(%id = 194741132)
A := Matrix(5,5,(i,j) -> a||i||j):
B := Matrix(5,5,(i,j) -> b||i||j):
t := time():
for i from 1 to 1000 do
myKroneckerProduct(A,B)
end do:
time() - t;
14.383
restart:
myKroneckerProduct := proc(A::Matrix, B::Matrix) # by DJ Keenan, v.2
local a, b, i, j;
# description "Returns the Kronecker product of A and B";
uses ArrayTools, LinearAlgebra;
a:= Dimension(A); b:= Dimension(B);
if member(0,[a,b]) then Matrix(a[1]*b[1],a[2]*b[2])
else Concatenate(1, seq(Concatenate(2, seq(A[i,j]*B, j=1..a[2])), i=1..a[1]))
end if;
end proc:
A := Matrix([[1,2,3],[4,5,6]]):
B := Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]]):
myKroneckerProduct(A,B);
Matrix(%id = 188660236)
A := Matrix(5,5,(i,j) -> a||i||j):
B := Matrix(5,5,(i,j) -> b||i||j):
t := time():
for i from 1 to 1000 do
myKroneckerProduct(A,B)
end do:
time() - t;
18.080
restart:
myKroneckerProduct := proc(A::Matrix, B::Matrix) # by DJ Keenan, v.3
local a, b;
# description "Returns the Kronecker product of A and B";
uses LinearAlgebra;
a:= Dimension(A); b:= Dimension(B);
Matrix(a[1]*b[1], a[2]*b[2], (i,j)-> A[1+iquo(i-1,b[1]), 1+iquo(j-1,b[2])]*B[1+irem(i-1,b[1]), 1+irem(j-1,b[2])])
end proc:
A := Matrix([[1,2,3],[4,5,6]]):
B := Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]]):
myKroneckerProduct(A,B);
Matrix(%id = 186680972)
A := Matrix(5,5,(i,j) -> a||i||j):
B := Matrix(5,5,(i,j) -> b||i||j):
t := time():
for i from 1 to 1000 do
myKroneckerProduct(A,B)
end do:
time() - t;
3.838
restart:
myKroneckerProduct := proc(A, B) # by Robert Israel
# options `Maple Advisor Database 1.01 for Maple 6`,
# `Copyright (c) 1998 by Robert B. Israel. All rights reserved`;
local Ap, Bp, i,j;
if nargs > 2 then RETURN(myKroneckerProduct(myKroneckerProduct(A,B),args[3..nargs])) fi;
if type(A,{vector,list(algebraic)}) and type(B,{vector,list(algebraic)}) then
vector([seq(seq(A[i]*B[j], j=1..linalg[vectdim](B)), i=1..linalg[vectdim](A))])
else # otherwise result is matrix
if type(A,matrix) then Ap:= A
elif type(A,listlist) then Ap:= convert(A,matrix)
elif type(A,list) then Ap:= matrix(map(t->[t],A))
elif type(A,specfunc(list,transpose)) then Ap:= matrix([op(A)])
else Ap:= convert(A,matrix)
fi;
if type(B,matrix) then Bp:= B
elif type(B,listlist) then Bp:= convert(B,matrix)
elif type(B,list) then Bp:= matrix(map(t->[t],B))
elif type(B,specfunc(list,transpose)) then Bp:= matrix([op(B)])
else Bp:= convert(B,matrix)
fi;
linalg[stackmatrix](seq(linalg[augment](seq(linalg[scalarmul](Bp,Ap[i,j]),
j = 1 .. linalg[coldim](Ap))),
i = 1 .. linalg[rowdim](Ap)));
fi
end proc:
A := Matrix([[1,2,3],[4,5,6]]):
B := Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]]):
myKroneckerProduct(A,B);
A := Matrix(5,5,(i,j) -> a||i||j):
B := Matrix(5,5,(i,j) -> b||i||j):
t := time():
for i from 1 to 1000 do
myKroneckerProduct(A,B)
end do:
time() - t;
41.652
restart:
myKroneckerProduct := proc(A::Matrix,B::Matrix) # by Preben Alsholm
local m,n;
m,n:=LinearAlgebra:-Dimensions(A);
Matrix([seq([seq(A[j,i]*B,i=1..n)],j=1..m)])
end proc:
A := Matrix([[1,2,3],[4,5,6]]):
B := Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]]):
myKroneckerProduct(A,B);
Matrix(%id = 204765580)
A := Matrix(5,5,(i,j) -> a||i||j):
B := Matrix(5,5,(i,j) -> b||i||j):
t := time():
for i from 1 to 1000 do
myKroneckerProduct(A,B)
end do:
time() - t;
34.851
restart:
myKroneckerProduct := module() # by Alec Mihailovs
export `*`;
option package;
`*`:=overload([proc(A::rtable,B::rtable) option overload;
rtable(<seq(seq(`*`(A[i,[rtable_dims(A)][2..-1][]],
B[j,[rtable_dims(B)][2..-1][]]),
j=[rtable_dims(B)][1]),i=[rtable_dims(A)][1])>,
select(has,{rtable_options(A)} intersect {rtable_options(B)},
{datatype,subtype,order,rectangular,sparse})[]) end,
:-`*`])
end:
with(myKroneckerProduct):
A := Matrix([[1,2,3],[4,5,6]]):
B := Matrix([[1,2,3],[4,5,6],[7,8,9],[10,11,12]]):
A*B;
Matrix(%id = 168930508)
A := Matrix(5,5,(i,j) -> a||i||j):
B := Matrix(5,5,(i,j) -> b||i||j):
t := time():
for i from 1 to 1000 do
A*B
end do:
time() - t;
26.676
