[转载]读Matlab的kron(下)
[p=21, 2, left] 在[url=http://xialulee.spaces.live.com/blog/cns!4EE324C8ACFA82DB!1082.entry]《读Matlab的kron.m(上)》[/url]中,研究了kron.m怎样处理full矩阵的Kronecker积,在本篇中,将要研究kron.m是如何处理稀疏矩阵(sparse)的。还是先看源码(Matlab版本7.6.0):[/p][p=21, 2, left][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] 这个函数的主要部分就是由一个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) [color=#FF0000](1)[/color];
这就是程序的思路;
程序按照这个思路做了:构造行向量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)是典型的计算二元函数的手法,我以前绘制复变函数曲面的时候用的就是这种手法(参见[url=http://xialulee.spaces.live.com/blog/cns!4EE324C8ACFA82DB!614.entry]《绘制一个复变函数》[/url]);
这样计算出来的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积。[/p]
页:
[1]