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积,我看了代码,简短而有效,先列出代码,然后作简要分析。
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。但是如果以两个向量作为下标对矩阵进行索引,得到的是什么呢?做实验如下:
>> A=magic(5)
A =17 24 1 8 1523 5 7 14 164 6 13 20 2210 12 19 21 311 18 25 2 9
>> ia=[1 1 2];ja=[1 3 5];
>> A(ia,ja)
ans =17 1 1517 1 1523 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进行索引,得到的新矩阵就是:
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代码
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积。
转自:http://blog.sina.com.cn/s/blog_4513dde60100ofoy.html
http://blog.sina.com.cn/s/blog_4513dde60100ofox.html
Matlab学习笔记 kron函数相关推荐
- Matlab学习笔记 figure函数
Matlab学习笔记 figure函数 matlab中的 figure 命令,能够创建一个用来显示图形输出的一个窗口对象.每一个这样的窗口都有一些属性,例如窗口的尺寸.位置,等等.下面一一介绍它们. ...
- matlab 调用子函数返回值,matlab学习笔记13_1 函数返回值
一起来学matlab-matlab学习笔记13函数 13_1 函数返回值 觉得有用的话,欢迎一起讨论相互学习~Follow Me 函数返回一个值 返回值不必使用return语句,而是直接将需要返回的变 ...
- [MATLAB学习笔记]peaks函数1013(2)
>> Z = peaksZ =1 至 10 列0.0001 0.0001 0.0002 0.0004 0.0007 0.0011 0.0017 0.0025 0.0034 0.00430. ...
- Matlab学习笔记——find()函数
写在这里的初衷,一是备忘,二是希望得到高人指点,三是希望能遇到志同道合的朋友. 目录 find 1.功能 2.格式 3.说明 拓展 1.稀疏矩阵 2.魔方矩阵 find 1.功能 查找非零元素的值和下 ...
- MATLAB学习笔记 :函数文件的定义和使用
数学建模比赛MATLAB从入门到精通教程_哔哩哔哩_bilibili function语法 1.编写函数文件,求半径为r的圆的面积和周长 (1)新建->函数 (2)编辑代码,保存 (3)回命令行 ...
- matlab学习笔记 bsxfun函数
最近总是遇到 bsxfun这个函数,前几次因为无关紧要只是大概看了一下函数体去对比结果,今天再一次遇见了这个函数,想想还是有必要掌握的,遂查了些资料总结如下. 函数bsxfun [功能描述]两个数组间 ...
- [MATLAB学习笔记]Rng函数
'twister':梅森旋转 'simdTwister':面向 SIMD 的快速梅森旋转算法 'combRecursive':组合多递归 'philox':执行 10 轮的 Philox 4×32 生 ...
- Matlab学习笔记 figure函数
matlab中的 figure 命令,能够创建一个用来显示图形输出的一个窗口对象.每一个这样的窗口都有一些属性,例如窗口的尺寸.位置,等等.下面一一介绍它们. 一.概述 总的来说,figure 的使 ...
- linspace函数matlab_从零开始的matlab学习笔记——(29)泰勒逼近函数
matlab应用--求极限,求导,求积分,解方程,概率统计,函数绘图,三维图像,拟合函数,动态图....更多内容尽在个人专栏:matlab学习 上一节我们成功制作了能自己转圈的三维螺旋线,这里我们再来 ...
- matlab数组平方的计算自定义函数_从零开始的matlab学习笔记——(38)简单数论计算函数:取整,gcd,lcm,质数,全排列...
matlab应用--求极限,求导,求积分,解方程,概率统计,函数绘图,三维图像,拟合函数,动态图,傅里叶变换,随机数,优化问题....更多内容尽在个人专栏:matlab学习 翻了翻优化工具箱,发现内容 ...
最新文章
- html判断对错,Html翻转校园试题
- FAT16文件系统结构扇区数据分析
- 精通python网络爬虫-精通Python网络爬虫:核心技术、框架与项目实战
- C# rdlc 报表学习总结
- 顺序表的所有基本操作
- 过了一会的gduuu
- redis最基础的入门教程
- html/css静态网页制作
- android 可执行程序 root权限,非Root权限的Android上运行可执行文件
- linux系统的总父目录,Linux虚拟文件系统-资料路径名的解析(2)-回退父目录
- 学好英语对IT软件工程师的影响
- 计算机一级中替换,08年计算机一级辅导:实战WPS转义符在查找替换中的应用
- 数据结构与算法-进阶(十二)最短路径Dijkstra 算法
- 通信原理包络是什么意思_为什么齿轮不能少于17个齿,少于17个齿,齿轮传动会如何?...
- 女儿拿着小天才电话手表问我App启动流程
- 计算机无法打开命令,电脑点击运行cmd打不开怎么办
- 比较motif和一条长序列的相似性
- 论成长型思维的重要性
- Cloud Paks地理数据研究成果|IBM
- 计算机毕业设计(附源码)python校园社团管理系统