matlab中的kron函数

函数 kron
格式 C=kron (A,B) %A为m×n矩阵，B为p×q矩阵，则C为mp×nq矩阵。

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积，我看了代码，简短而有效，先列出代码，然后作简要分析。

[plain] view plaincopy

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后面的代码。
无论是处理full还是sparse，四个变量ma,mb,na,nb都被使用着，它们分别表示矩阵A的行数，矩阵B的行数，矩阵A的列数，矩阵B的列数，也就是矩阵A和B的尺寸信息；

(1)如果都不是稀疏矩阵：

先分析代码的上半部分，即针对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。但是如果以两个向量作为下标对矩阵进行索引，得到的是什么呢？做实验如下：

[plain] view plaincopy

>> 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进行索引，得到的新矩阵就是：

[plain] view plaincopy

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进行索引，得到的新矩阵是：

[plain] view plaincopy

B B … B
B B … B
……
B B … B

每一“行”有na个B，每一“列”有ma个B；
A(ia,ja).*B(ib,jb)是两个新矩阵的对应元素乘积，新矩阵中的“一块”可表示如下：

[plain] view plaincopy

Ci,j.*B
=Ai,j*ones(mb,nb).*B
=Ai,j*B

这就是Kronecker积的定义。

(2)稀疏矩阵的Kronecker代码

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=
A_1,1*B   A_1,2*B   …   A_1,na*B
A_2,1*B   A_2,2*B   …   A_2,na*B
……
A_ma,1*B  A_ma,2*B  …   A_ma,na*B
      假设矩阵B中w行v列有一个非零元素b，那么在进行Kronecker积的时候，它会出现在每一个分块A_i,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]的元素等于A_i,j*B_w,v；
      由于处理的是稀疏矩阵，零元素不会保存于结果之中，因此上式改写成：
      K_{w+(i-1)*mb,v+(j-1)*nb}=A_i,j*B_w,v(A_i,j≠0;B_w,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中非零元素是：K_{w+(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=ik_p,q,n=jk_p,q；
      (sb*sa.')_p,q为B中第p个非零元素（稀疏矩阵只存放非零元素，这里“第p个”表示将非零元素排成一排，取出第p个，假设该元素行号为w，列号为v）与A中第q个非零元素（假设行号为i，列号为j）的积（B_w,v*A_i,j），而根据ik和jk的构造方式，可以算出ik_p,q=w+mb*(i-1)，jk_p,q=v+nb*(j-1)，这个结果符合式(1)的要求，即上述方法正确计算了Kronecker积。

转自：http://blog.sina.com.cn/s/blog_4513dde60100ofoy.html

http://blog.sina.com.cn/s/blog_4513dde60100ofox.html

matlab中的kron函数相关推荐

matlab计算两向量的乘积,matlab中两个函数相乘
变量名最多不超过63个字符; ? 变量名区分大小写; ? Matlab提供的标准函数名以及命令名必须用小写字母; ? 变量名中不能包含空格.标点.运算符. 1.变量及其...... 中的元素; (2) ...
MATLAB中的常用函数小结
1. MATLAB中的常用函数小结文章目录 1. MATLAB中的常用函数小结 1. MATLAB图像处理工具箱 1.1 图像显示 1.2 图像文件输入/输出 1.3. 图像像素值及其统计 1.4 ...
Matlab中的lsqcurvefit函数的使用
Matlab中的lsqcurvefit函数的使用 lsqcurvefit函数调用示例 lsqcurvefit函数非线性曲线拟合是已知输入向量xdata和输出向量ydata,并且知道输入与输出的函数 ...
Matlab：Matlab中常用的函数、案例详细攻略
Matlab:Matlab中常用的函数.案例详细攻略目录常用函数 1.与文件相关 2.MATLAB GUI不同控件函数间变量传递方法常用函数 Matlab中的bwmorph函数解释 bwmorp ...
matlab作动态函数曲线图,[转载]Matlab中使用Plot函数动态画图方法总结
本帖最后由 sonictl 于 2012-12-31 12:18 编辑请删除我清楚超靠靠靠没办法,一会儿限制这不能发表,那不能发表的.... [转载]Matlab中使用Plot函数动态画图方法总 ...
python实现Matlab中的circshift函数
circshift是Matlab中矩阵循环移位函数,具体使用参照该链接. 但是python中并没有封装好的该函数,因此需要自己实现. 思路:将矩阵分为两部分,然后按照自己的需要堆叠在一起就可以了. n ...
matlab的数学函数,matlab中常见数学函数的使用
matlab中常见数学函数的使用 MATLAB 基本知识 Matlab 的内部常数 pi 圆周率 exp(1) 自然对数的底数 e i 或 j 虚数单位 Inf 或 inf 无穷大 Matlab 的常 ...
matlab里inline定义矩阵,Matlab中的inline函数_matlab中inline函数
Matlab中的inline函数 1．有时为了描述某个数学函数的方便,可以用inline()函数来直接编写该函数,形式相当于M-函数,但无编写一个真正的MATLAB文件,就可以描述出某种数学关系.其调 ...
Matlab中的eig函数和Opecv中eigen（）函数的区别
奇异值分解的理论参见下面的链接 http://www.cnblogs.com/pinard/p/6251584.html https://blog.csdn.net/shenziheng1/artic ...
MATLAB中神经网络train函数使用说明
MATLAB中神经网络train( )函数使用说明函数的语法格式如下: [net, tr]=train(net, P, T, Pi, Ai): train( )函数用于训练创建好的感知器网络,事实上 ...

matlab中的kron函数

matlab中的kron函数相关推荐

最新文章

热门文章