[转载]读Matlab的kron（下）

    在《读Matlab的kron.m（上）》中，研究了kron.m怎样处理full矩阵的Kronecker积，在本篇中，将要研究kron.m是如何处理稀疏矩阵（sparse）的。
      还是先看源码（Matlab版本7.6.0）：

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

      这个函数的主要部分就是由一个if-else组成的。if用来判断两个输入的矩阵中是否有稀疏矩阵，只要有一个是稀疏的，那么就跳转到else部分执行代码。所以else后面的代码都是针对稀疏矩阵的；
      无论是处理full还是sparse，四个变量ma,mb,na,nb都被使用着，它们分别表示矩阵A的行数，矩阵B的行数，矩阵A的列数，矩阵B的列数，也就是矩阵A和B的尺寸信息；
      else中，首先用find函数取得了A和B中非零元素的信息：ia表示A中非零元素的行号，ja表示A中非零元素的列号，sa表示A中非零元素的取值；ib表示B中非零元素的行号，jb表示B中非零元素的列号，sb表示B中非零元素的取值。然后利用ia=ia(:)这样的方式将这六个变量转化为列向量的形式；
      接着用ones构造了两个尺寸分别与sa和sb相同的全1列向量ka和kb；
      构造行向量t=mb*(ia-1)';
      用t和ib构造矩阵ik；
      构造行向量t=nb*(ja-1)';
      用t和jb构造矩阵jk；
      然后把这些东西利用sparse函数构造出了Kronecker积的结果。
      如果不说这段代码是干什么的，我肯定会晕头转向。
      下面分析一下为什么这就是计算Kronecker积：
      K=kron(A,B)的结果是：
K=
A1,1*B   A1,2*B   …   A1,na*B
A2,1*B   A2,2*B   …   A2,na*B
……
Ama,1*B  Ama,2*B  …   Ama,na*B
      假设矩阵B中w行v列有一个非零元素b，那么在进行Kronecker积的时候，它会出现在每一个分块Ai,j*B的w行v列，而这样的分块有(ma*mb)*(na*nb)个，于是所有b出现的行号可以表示成w+(i-1)*mb，b出现的列号可以表示成v+(j-1)*nb。
      所以，在K中，所有位于[w+(i-1)*mb,v+(j-1)*nb]的元素等于Ai,j*Bw,v；
      由于处理的是稀疏矩阵，零元素不会保存于结果之中，因此上式改写成：
      Kw+(i-1)*mb,v+(j-1)*nb=Ai,j*Bw,v(Ai,j≠0;Bw,v≠0)      (1)；
      这就是程序的思路；
      程序按照这个思路做了：构造行向量t=mb*(ia-1)'，并将它按行复制，B中有多少非零元素，就将t复制成多少行的矩阵：t(kb,:)。（kb是长度为size(sb)的全1向量）；
      把B中所有非零元素的行号作为一个列向量，然后按列复制，A中有多少个非零元素，就复制多少列，然后：
ik=t(kb,:)+ib(:,ka)
      这个操作非常类似meshgrid。回头看一眼式(1)，kron积的结果K中非零元素是：
Kw+(i-1)*mb,v+(j-1)*nb
      非零元素的行号w+(i-1)*mb是B中非零元素的行号w和A中非零元素的行号i的函数，即K中所非零元素的行号是通过二元函数计算出来的，上面类似meshgrid的操作ik=t(kb,:)+ib(:,ka)是典型的计算二元函数的手法，我以前绘制复变函数曲面的时候用的就是这种手法（参见《绘制一个复变函数》）；
      这样计算出来的ik就包含了K中所有非零元素的行号，非零元素的列号计算与此相同，就不赘述了；
      得到了K中所有非零元素的行号和列号，如果再知道它们的值，K就计算出来了。怎样算呢？根据式(1)即可。程序中这样实现：
K = sparse(ik,jk,sb*sa.',ma*mb,na*nb);
      sparse构造矩阵的方式如下：
      将矩阵(sb*sa.')的p行q列元素作为K的m行n列的元素，其中m=ikp,q,n=jkp,q；
      (sb*sa.')p,q为B中第p个非零元素（稀疏矩阵只存放非零元素，这里“第p个”表示将非零元素排成一排，取出第p个，假设该元素行号为w，列号为v）与A中第q个非零元素（假设行号为i，列号为j）的积（Bw,v*Ai,j），而根据ik和jk的构造方式，可以算出ikp,q=w+mb*(i-1)，jkp,q=v+nb*(j-1)，这个结果符合式(1)的要求，即上述方法正确计算了Kronecker积。