matlab按照列kron,Matlab学习笔记kron函数

函数 kron

格式 C&＃61;kron (A,B) %A为m×n矩阵&＃xff0c;B为p×q矩阵&＃xff0c;则C为mp×nq矩阵。

kron即为Kronecker积&＃xff0c;所谓Kronecker积是一种矩阵运算&＃xff0c;其定义可以简单描述成&＃xff1a;

X与Y的Kronecker积的结果是一个矩阵&＃xff1a;

X11*Y X12*Y … X1n*Y

X21*Y X22*Y … X2n*Y

……

Xm1*Y Xm2*Y … Xmn*Y

Matlab中有kron函数用来计算Kronecker积&＃xff0c;我看了代码&＃xff0c;简短而有效&＃xff0c;先列出代码&＃xff0c;然后作简要分析。

function K &＃61; 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] &＃61; size(A);

[mb,nb] &＃61; size(B);

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

% Both inputs full, result is full.

[ia,ib] &＃61; meshgrid(1:ma,1:mb);

[ja,jb] &＃61; meshgrid(1:na,1:nb);

K &＃61; A(ia,ja).*B(ib,jb);

else

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

[ia,ja,sa] &＃61; find(A);

[ib,jb,sb] &＃61; find(B);

ia &＃61; ia(:); ja &＃61; ja(:); sa &＃61; sa(:);

ib &＃61; ib(:); jb &＃61; jb(:); sb &＃61; sb(:);

ka &＃61; ones(size(sa));

kb &＃61; ones(size(sb));

t &＃61; mb*(ia-1)&＃39;;

ik &＃61; t(kb,:)&＃43;ib(:,ka);

t &＃61; nb*(ja-1)&＃39;;

jk &＃61; t(kb,:)&＃43;jb(:,ka);

K &＃61; sparse(ik,jk,sb*sa.&＃39;,ma*mb,na*nb);

end

这个函数的主要部分就是由一个if-else组成的。if用来判断两个输入的矩阵中是否有稀疏矩阵&＃xff0c;只要有一个是稀疏的&＃xff0c;那么就跳转到else部分执行代码。所以else后面的代码都是针对稀疏矩阵的&＃xff1b;首先通过size函数取得了A和B的尺寸信息&＃xff0c;ma、mb分别代表A和B的行数&＃xff0c;na、nb分别代表A和B的列数&＃xff1b;然后使用issparse判断矩阵的类型&＃xff0c;如果输入的两个矩阵A和B都不是稀疏矩阵&＃xff0c;则执行下面的代码&＃xff0c;否则执行else后面的代码。 无论是处理full还是sparse&＃xff0c;四个变量ma,mb,na,nb都被使用着&＃xff0c;它们分别表示矩阵A的行数&＃xff0c;矩阵B的行数&＃xff0c;矩阵A的列数&＃xff0c;矩阵B的列数&＃xff0c;也就是矩阵A和B的尺寸信息&＃xff1b;

(1)如果都不是稀疏矩阵&＃xff1a;

先分析代码的上半部分&＃xff0c;即针对full矩阵的运算。首先通过size函数取得了A和B的尺寸信息&＃xff0c;ma、mb分别代表A和B的行数&＃xff0c;na、nb分别代表A和B的列数&＃xff1b;然后使用issparse判断矩阵的类型&＃xff0c;如果输入的两个矩阵A和B都不是稀疏矩阵&＃xff0c;则执行下面的代码&＃xff0c;否则执行else后面的代码。

两次调用meshgrid构造了4个矩阵ia,ib,ja,jb。其中ia是1:ma按行复制&＃xff0c;ib是1:mb按列复制&＃xff0c;ja是1:na按行复制&＃xff0c;jb是1:nb按列复制。构造两个矩阵A(ia,ja)和B(ib,jb)&＃xff0c;让这两个矩阵对应元素相乘&＃xff0c;得到的新矩阵就是A与B的Kronecker积。

对于full矩阵的kronecker积的代码就是这么少&＃xff0c;非常简洁&＃xff0c;但可能其原理不是那么一目了然。关键就在A(ia,ja)和B(ib,jb)究竟为何物。如果我们将两个数字作为矩阵的下标&＃xff0c;将会得到下标对应的矩阵元素&＃xff0c;例如A&＃61;eye(3);A(2,2)就是1。但是如果以两个向量作为下标对矩阵进行索引&＃xff0c;得到的是什么呢&＃xff1f;做实验如下&＃xff1a;

>> A&＃61;magic(5)

A &＃61;

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&＃61;[1 1 2];ja&＃61;[1 3 5];

>> A(ia,ja)

ans &＃61;

17 1 15

17 1 15

23 7 16

得到了一个新的矩阵。这个新矩阵是这样构建的&＃xff1a;

以ia的第一个元素作为行号&＃xff0c;以ja中的所有元素作为列号&＃xff0c;取出A中的对应元素&＃xff0c;将这些元素摆成一行&＃xff0c;构成新矩阵的第一行&＃xff1b;

以ia的第二个元素作为行号&＃xff0c;以ja中的所有元素作为列号&＃xff0c;取出A中的对应元素&＃xff0c;将这些元素摆成一行&＃xff0c;构成新矩阵的第二行&＃xff1b;

依此重复&＃xff0c;直到ia中所有的元素被取到。

如果ia和ja是矩阵呢&＃xff1f;Matlab会将矩阵形式的ia的列首尾相接&＃xff0c;使其变成一个向量&＃xff0c;ja也是同样处理。

回到程序中&＃xff0c;程序用meshgrid构造了四个矩阵ia,ja,ib,jb。它们的内容分别为&＃xff1a;

ia&＃xff1a; 每一行都是1:ma&＃xff0c;一共有mb行&＃xff1b;

ja&＃xff1a; 每一行都是1:na&＃xff0c;一共有nb行&＃xff1b;

ib&＃xff1a; 每一列都是(1:mb)&＃39;&＃xff0c;一共有ma列&＃xff1b;

jb&＃xff1a; 每一列都是(1:nb)&＃39;&＃xff0c;一共有na列。

之后就有了A(ia,ja)&＃xff0c;Matlab将ia的列首尾相接&＃xff0c;变成向量&＃xff0c;将ja的列首尾相接变成向量&＃xff0c;实际上A(ia,ja)就是A(ia(:),ja(:))&＃xff0c;用ia(:)和ja(:)作为下标对A进行索引&＃xff0c;得到的新矩阵就是&＃xff1a;

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)&＃xff1b;

用ib(:)和jb(:)作为下标对B进行索引&＃xff0c;得到的新矩阵是&＃xff1a;

B B … B

B B … B

……

B B … B每一“行”有na个B&＃xff0c;每一“列”有ma个B&＃xff1b;

A(ia,ja).*B(ib,jb)是两个新矩阵的对应元素乘积&＃xff0c;新矩阵中的“一块”可表示如下&＃xff1a;

Ci,j.*B

&＃61;Ai,j*ones(mb,nb).*B

&＃61;Ai,j*B这就是Kronecker积的定义。

(2)稀疏矩阵的Kronecker代码

else中&＃xff0c;首先用find函数取得了A和B中非零元素的信息&＃xff1a;

ia表示A中非零元素的行号&＃xff0c;ja表示A中非零元素的列号&＃xff0c;sa表示A中非零元素的取值&＃xff1b;ib表示B中非零元素的行号&＃xff0c;jb表示B中非零元素的列号&＃xff0c;sb表示B中非零元素的取值。然后利用ia&＃61;ia(:)这样的方式将这六个变量转化为列向量的形式&＃xff1b;

接着用ones构造了两个尺寸分别与sa和sb相同的全1列向量ka和kb&＃xff1b; 构造行向量t&＃61;mb*(ia-1)&＃39;; 用t和ib构造矩阵ik&＃xff1b; 构造行向量t&＃61;nb*(ja-1)&＃39;; 用t和jb构造矩阵jk&＃xff1b; 然后把这些东西利用sparse函数构造出了Kronecker积的结果。 如果不说这段代码是干什么的&＃xff0c;我肯定会晕头转向。 下面分析一下为什么这就是计算Kronecker积&＃xff1a; K&＃61;kron(A,B)的结果是&＃xff1a; K&＃61; A1,1*BA1,2*B…A1,na*B A2,1*BA2,2*B…A2,na*B …… Ama,1*BAma,2*B…Ama,na*B 假设矩阵B中w行v列有一个非零元素b&＃xff0c;那么在进行Kronecker积的时候&＃xff0c;它会出现在每一个分块Ai,j*B的w行v列&＃xff0c;而这样的分块有(ma*mb)*(na*nb)个&＃xff0c;于是所有b出现的行号可以表示成w&＃43;(i-1)*mb&＃xff0c;b出现的列号可以表示成v&＃43;(j-1)*nb。 所以&＃xff0c;在K中&＃xff0c;所有位于[w&＃43;(i-1)*mb,v&＃43;(j-1)*nb]的元素等于Ai,j*Bw,v&＃xff1b; 由于处理的是稀疏矩阵&＃xff0c;零元素不会保存于结果之中&＃xff0c;因此上式改写成&＃xff1a; Kw&＃43;(i-1)*mb,v&＃43;(j-1)*nb&＃61;Ai,j*Bw,v(Ai,j≠0;Bw,v≠0)(1)&＃xff1b; 这就是程序的思路&＃xff1b; 程序按照这个思路做了&＃xff1a;构造行向量t&＃61;mb*(ia-1)&＃39;&＃xff0c;并将它按行复制&＃xff0c;B中有多少非零元素&＃xff0c;就将t复制成多少行的矩阵&＃xff1a;t(kb,:)。(kb是长度为size(sb)的全1向量)&＃xff1b; 把B中所有非零元素的行号作为一个列向量&＃xff0c;然后按列复制&＃xff0c;A中有多少个非零元素&＃xff0c;就复制多少列&＃xff0c;然后&＃xff1a;ik&＃61;t(kb,:)&＃43;ib(:,ka) 这个操作非常类似meshgrid。回头看一眼式(1)&＃xff0c;kron积的结果K中非零元素是&＃xff1a;Kw&＃43;(i-1)*mb,v&＃43;(j-1)*nb 非零元素的行号w&＃43;(i-1)*mb是B中非零元素的行号w和A中非零元素的行号i的函数&＃xff0c;即K中所非零元素的行号是通过二元函数计算出来的&＃xff0c;上面类似meshgrid的操作ik&＃61;t(kb,:)&＃43;ib(:,ka)是典型的计算二元函数的手法&＃xff1b; 这样计算出来的ik就包含了K中所有非零元素的行号&＃xff0c;非零元素的列号计算与此相同&＃xff0c;就不赘述了&＃xff1b; 得到了K中所有非零元素的行号和列号&＃xff0c;如果再知道它们的值&＃xff0c;K就计算出来了。怎样算呢&＃xff1f;根据式(1)即可。程序中这样实现&＃xff1a; K &＃61; sparse(ik,jk,sb*sa.&＃39;,ma*mb,na*nb); sparse构造矩阵的方式如下&＃xff1a; 将矩阵(sb*sa.&＃39;)的p行q列元素作为K的m行n列的元素&＃xff0c;其中m&＃61;ikp,q,n&＃61;jkp,q&＃xff1b; (sb*sa.&＃39;)p,q为B中第p个非零元素(稀疏矩阵只存放非零元素&＃xff0c;这里“第p个”表示将非零元素排成一排&＃xff0c;取出第p个&＃xff0c;假设该元素行号为w&＃xff0c;列号为v)与A中第q个非零元素(假设行号为i&＃xff0c;列号为j)的积(Bw,v*Ai,j)&＃xff0c;而根据ik和jk的构造方式&＃xff0c;可以算出ikp,q&＃61;w&＃43;mb*(i-1)&＃xff0c;jkp,q&＃61;v&＃43;nb*(j-1)&＃xff0c;这个结果符合式(1)的要求&＃xff0c;即上述方法正确计算了Kronecker积。

转自&＃xff1a;http://blog.sina.com.cn/s/blog_4513dde60100ofoy.html

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