[转载]读Matlab的kron（上）
[p=21, 2, left] 多学习他人的代码有助于提高编程水平。这几天看一篇关于天线阵列的论文，里面使用了Kronecker积的概念。
所谓Kronecker积是一种矩阵运算，其定义可以简单描述成：
X与Y的Kronecker积的结果是一个矩阵：
X11*Y X12*Y … X1n*Y
X21*Y X22*Y … X2n*Y
……
Xm1*Y Xm2*Y … Xmn*Y
Matlab中有kron函数用来计算Kronecker积，我看了代码，简短而有效，先列出代码，然后作简要分析。
[color=#0000FF]function[/color] K = kron(A,B)
[color=#228B22]%KRON Kronecker tensor product.[/color]
[color=#228B22]% KRON(X,Y) is the Kronecker tensor product of X and Y.[/color]
[color=#228B22]% The result is a large matrix formed by taking all possible[/color]
[color=#228B22]% products between the elements of X and those of Y. For[/color]
[color=#228B22]% example, if X is 2 by 3, then KRON(X,Y) is[/color]
[color=#228B22]%[/color]
[color=#228B22]% [ X(1,1)*Y X(1,2)*Y X(1,3)*Y[/color]
[color=#228B22]% X(2,1)*Y X(2,2)*Y X(2,3)*Y ][/color]
[color=#228B22]%[/color]
[color=#228B22]% If either X or Y is sparse, only nonzero elements are multiplied[/color]
[color=#228B22]% in the computation, and the result is sparse.[/color]
[color=#228B22]%[/color]
[color=#228B22]% Class support for inputs X,Y:[/color]
[color=#228B22]% float: double, single [/color]
[color=#228B22]% Previous versions by Paul Fackler, North Carolina State,[/color]
[color=#228B22]% and Jordan Rosenthal, Georgia Tech.[/color]
[color=#228B22]% Copyright 1984-2004 The MathWorks, Inc. [/color]
[color=#228B22]% $Revision: 5.17.4.2 $ $Date: 2004/06/25 18:52:18 $ [/color]
[ma,na] = size(A);
[mb,nb] = size(B);
[color=#0000FF]if[/color] ~issparse(A) && ~issparse(B)
[color=#228B22]% Both inputs full, result is full. [/color]
[ia,ib] = meshgrid(1:ma,1:mb);
[ja,jb] = meshgrid(1:na,1:nb);
K = A(ia,ja).*B(ib,jb);
[color=#0000FF]else[/color]
[color=#228B22]% At least one input is sparse, result is sparse. [/color]
[ia,ja,sa] = find(A);
[ib,jb,sb] = find(B);
ia = ia(:); ja = ja(:); sa = sa(:);
ib = ib(:); jb = jb(:); sb = sb(:);
ka = ones(size(sa));
kb = ones(size(sb));
t = mb*(ia-1)';
ik = t(kb,:)+ib(:,ka);
t = nb*(ja-1)';
jk = t(kb,:)+jb(:,ka);
K = sparse(ik,jk,sb*sa.',ma*mb,na*nb);
[color=#0000FF]end[/color]
[/p][p=21, 2, left] 先分析代码的上半部分，即针对full矩阵的运算。
首先通过size函数取得了A和B的尺寸信息，ma、mb分别代表A和B的行数，na、nb分别代表A和B的列数；
然后使用issparse判断矩阵的类型，如果输入的两个矩阵A和B都不是稀疏矩阵，则执行下面的代码，否则执行else后面的代码。
如果都不是稀疏矩阵：
两次调用meshgrid构造了4个矩阵ia,ib,ja,jb。其中ia是1:ma按行复制，ib是1:mb按列复制，ja是1:na按行复制，jb是1:nb按列复制。[/p][p=21, 2, left] 构造两个矩阵A(ia,ja)和B(ib,jb)，让这两个矩阵对应元素相乘，得到的新矩阵就是A与B的Kronecker积。
对于full矩阵的kronecker积的代码就是这么少，非常简洁，但可能其原理不是那么一目了然。关键就在A(ia,ja)和B(ib,jb)究竟为何物。
如果我们将两个数字作为矩阵的下标，将会得到下标对应的矩阵元素，例如A=eye(3);A(2,2)就是1。但是如果以两个向量作为下标对矩阵进行索引，得到的是什么呢？做实验如下：
>> A=magic(5)[/p][p=21, 2, left]A =[/p][p=21, 2, left] 17 24 1 8 15
23 5 7 14 16
4 6 13 20 22
10 12 19 21 3
11 18 25 2 9[/p][p=21, 2, left]>> ia=[1 1 2];ja=[1 3 5];
>> A(ia,ja)[/p][p=21, 2, left]ans =[/p][p=21, 2, left] 17 1 15
17 1 15
23 7 16[/p][p=21, 2, left] 得到了一个新的矩阵。这个新矩阵是这样构建的：
以ia的第一个元素作为行号，以ja中的所有元素作为列号，取出A中的对应元素，将这些元素摆成一行，构成新矩阵的第一行；
以ia的第二个元素作为行号，以ja中的所有元素作为列号，取出A中的对应元素，将这些元素摆成一行，构成新矩阵的第二行；
依此重复，直到ia中所有的元素被取到。
如果ia和ja是矩阵呢？Matlab会将矩阵形式的ia的列首尾相接，使其变成一个向量，ja也是同样处理。
回到程序中，程序用meshgrid构造了四个矩阵ia,ja,ib,jb。它们的内容分别为：
ia：
每一行都是1:ma，一共有mb行；
ja：
每一行都是1:na，一共有nb行；
ib：
每一列都是(1:mb)'，一共有ma列；
jb：
每一列都是(1:nb)'，一共有na列。
之后就有了A(ia,ja)，Matlab将ia的列首尾相接，变成向量，将ja的列首尾相接变成向量，实际上A(ia,ja)就是A(ia(:),ja(:))，用ia(:)和ja(:)作为下标对A进行索引，得到的新矩阵就是：
C1,1 C1,2 … C1,na
C2,1 C2,2 … C2,na
……
Cma,1 Cma,2… Cma,na
其中Ci,j是矩阵Ai,j*ones(mb,nb)；
用ib(:)和jb(:)作为下标对B进行索引，得到的新矩阵是：
B B … B
B B … B
……
B B … B
每一“行”有na个B，每一“列”有ma个B；
A(ia,ja).*B(ib,jb)是两个新矩阵的对应元素乘积，新矩阵中的“一块”可表示如下：
Ci,j.*B
=Ai,j*ones(mb,nb).*B
=Ai,j*B
这就是Kronecker积的定义。
上面分析的是full矩阵部分，稀疏矩阵的Kronecker代码待到下篇再进行分析。[/p]
页:
[1]