[转载]读Matlab的kron（上）

   多学习他人的代码有助于提高编程水平。
      这几天看一篇关于天线阵列的论文，里面使用了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积，我看了代码，简短而有效，先列出代码，然后作简要分析。
function K = kron(A,B)
%KRON   Kronecker tensor product.
%   KRON(X,Y) is the Kronecker tensor product of X and Y.
%   The result is a large matrix formed by taking all possible
%   products between the elements of X and those of Y.   For
%   example, if X is 2 by 3, then KRON(X,Y) is
%
%      [ X(1,1)*Y  X(1,2)*Y  X(1,3)*Y
%        X(2,1)*Y  X(2,2)*Y  X(2,3)*Y ]
%
%   If either X or Y is sparse, only nonzero elements are multiplied
%   in the computation, and the result is sparse.
%
%   Class support for inputs X,Y:
%      float: double, single

%   Previous versions by Paul Fackler, North Carolina State,
%   and Jordan Rosenthal, Georgia Tech.
%   Copyright 1984-2004 The MathWorks, Inc.
%   $Revision: 5.17.4.2 $ $Date: 2004/06/25 18:52:18 $

[ma,na] = size(A);
[mb,nb] = size(B);

if ~issparse(A) && ~issparse(B)

   % Both inputs full, result is full.

   [ia,ib] = meshgrid(1:ma,1:mb);
   [ja,jb] = meshgrid(1:na,1:nb);
   K = A(ia,ja).*B(ib,jb);

else

   % At least one input is sparse, result is sparse.

   [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);

end

      先分析代码的上半部分，即针对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按列复制。

      构造两个矩阵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)

A =

    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

>> ia=[1 1 2];ja=[1 3 5];
>> A(ia,ja)

ans =

    17     1    15
    17     1    15
    23     7    16

      得到了一个新的矩阵。这个新矩阵是这样构建的：
      以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代码待到下篇再进行分析。