函数 kron
格式 C=kron (A,B)
kron即为Kronecker积,所谓Kronecker积是一种矩阵运算,其定义可以简单描述成:
X与Y的Kronecker积的结果是一个矩阵:
X11*Y
X21*Y
……
Xm1*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
这个函数的主要部分就是由一个if-else组成的。if用来判断两个输入的矩阵中是否有稀疏矩阵,只要有一个是稀疏的,那么就跳转到else部分执行代码。所以else后面的代码都是针对稀疏矩阵的;首先通过size函数取得了A和B的尺寸信息,ma、mb分别代表A和B的行数,na、nb分别代表A和B的列数;然后使用issparse判断矩阵的类型,如果输入的两个矩阵A和B都不是稀疏矩阵,则执行下面的代码,否则执行else后面的代码。 (1)如果都不是稀疏矩阵:
先分析代码的上半部分,即针对full矩阵的运算。首先通过size函数取得了A和B的尺寸信息,ma、mb分别代表A和B的行数,na、nb分别代表A和B的列数;然后使用issparse判断矩阵的类型,如果输入的两个矩阵A和B都不是稀疏矩阵,则执行下面的代码,否则执行else后面的代码。
>> 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
得到了一个新的矩阵。这个新矩阵是这样构建的:
回到程序中,程序用meshgrid构造了四个矩阵ia,ja,ib,jb。它们的内容分别为:
ia:
ja:
ib:
jb:
C1,1 C1,2 … C1,na
C2,1 C2,2 … C2,na
……
Cma,1 Cma,2… Cma,na
其中C
i,j是矩阵A
i,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积的定义。
(2)稀疏矩阵的Kronecker代码
接着用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
A2,1*B
……
Ama,1*B
K = sparse(ik,jk,sb*sa.',ma*mb,na*nb);
转自:http://blog.sina.com.cn/s/blog_4513dde60100ofoy.html
http://blog.sina.com.cn/s/blog_4513dde60100ofox.html