无穷嵌套矩阵运算在Matlab中的实现.pdf

收稿日期 2010 09 25 基金项目 国家自然科学基金 10901003 黄山学院自然科学研究项目 2007xkj003 作者简介 沈来信 1979 男 安徽怀远人 硕士 讲师 研究方向 数据挖掘 人工智能等 第 32 卷第 1 期淮北师范大学学报 自然科学版 Vol 32No 1 2011 年 3 月 Journal of Huaibei Normal University Natural Science Mar 2011 无穷嵌套矩阵运算在 Matlab 中的实现 沈来信 1 任 勇 2 杨 帆 1 1 黄山学院信息工程学院 安徽 黄山 245021 2 安徽师范大学数学计算机科学学院 安徽 芜湖 241000 摘要 无穷嵌套矩阵的表示和有效运算一直是个难点 前者可以利用 Matlab 提供的单元结构描述 但是单元结构不 支持直接运算 利用逆分块矩阵的思想 实现了嵌套矩阵的加法 乘法等运算 实验表明 基于单元结构和逆分块矩阵 的思想可以有效地解决无穷嵌套矩阵运算问题 关键词 无穷嵌套矩阵 单元结构 逆分块矩阵 中图分类号 TP 393 09文献标识码 A 文章编号 2095 0691 2011 01 0063 03 Matlab 是一种交互式的以矩阵为基础的系统计算平台 它用于科学和工程的计算与可视化 它的优点在 于快速开发计算方法 而不在于计算速度 嵌套矩阵就是矩阵的每个元素又是一个矩阵 QBD Quasi birth and death 过程 1 和 Markov 模型 2 等经 常用到无穷循环的嵌套矩阵 常用的 C Java 等语言都很难准确表示这类矩阵以及它们的运算 本文利用 QBD 过程中典型的 Laplace 变换的求解来说明 Matlab 中矩阵迭代和嵌套矩阵的表示以及其运算的实现 如下面由公式 1 6 确定的变换公式 7 其中 A0 A1 A2 U u v 等都是 n 阶方阵 I 是 n 阶的单位 矩阵 R u v 是带有参数 u v 的 n 阶方阵 R u v A0 U 1 u v 1 G u v U 1 u v A2 2 U u v A1 vB uI R u v A2 3 1矩阵的迭代求法 可以采用如下递归过程实现公式 1 3 的迭代求解 经过一个预定的循环次数后 求出 R u v G u v 和 U u v 算法如下 初始化 U eye n n i 1 循环次数初始化 Iter 200 While i Iter U2 1 inv U R A0 U2 G U2 A2 U A1 v B u eye n R A2 i i 1 结束 64淮北师范大学学报 自然科学版 2011 年 2Matlab 中的单元结构 cell Matlab 中的单元结构 cell 能够实现把不同属性的数据纳入一个变量中 类似于具有 4 5 6 式的矩 阵和多维数组 我们用 cell M M 产生一个 M 阶的方阵 这样每个元素又可以为一个单元 公式 5 可用算 法 IR 实现 算法如下 Function IR IR n M R M 为矩阵 1 的维数 可以取一个很大的值 IR cell M M 定义 IR 为一个 M 阶的方阵 其中每个元素又是个矩阵 For j 1 M 先求出第一行 IR 1 j eye n R j 1 R u v 由相应的函数求出 End For i 2 M 再求出其他行 For j 1 M If i j IR i j eye n 对角线上位 n 阶的单位矩阵 Else If i j IR i j zeros n 下半角为 n 阶的零矩阵 Else IR i j IR i 1 j 从第二行开始 对角线右侧元素等于同列上一行元素 结束 if 和 for U1 1 u v diag U u v U u v U u v 4 I R u v 1 IR u vR u v 2 R u v 3 IR u v 2 R u v 3 IR u v 3 5 I G u v 1 I G u vI G u v 2 G u vI G u v 3 G u v 2 G u vI 6 对于无穷矩阵 我们采用一个很大的维数 M 来模拟 通过不断增加其值 得到最后无穷矩阵的近似值 通过观察矩阵的收敛特征 可知 n 增大到一个程度时 无穷嵌套矩阵中元素趋于 0 这样就把无穷矩阵变为 有穷矩阵进行计算 公式 4 和 6 可以采用类似的函数实现 这里不再一一叙述 3基于逆分块思想的 cell 结构运算 分块矩阵是用穿过矩阵的横线和竖线将矩阵 A 分割成若干个子块 以这些子块为元素的矩阵 A 在运 算时 先把子块当作元素运算 然后子块再运算 当计算 C A B 时 A 和 B 的分块矩阵要求满足 前一个 矩阵列的分法与后一个矩阵的行的分法完全一样 计算 C AB 时 A 和 B 的分块矩阵要求满足 A 的列块 数 B 的行块数 即 Aik的列数 Bkj的行数 由于 cell 结构不支持直接运算 用 cell 结构表示的嵌套矩阵就不能完成相应的加法和乘法 将采用逆分 块矩阵的思想 可以先把 cell 结构表示的嵌套矩阵用 Matlab 提供的 cell2mat 函数把 cell 结构转换为普通矩 阵的形式 然后参与运算 这就相当于把两个分块矩阵的加法和乘法转换成了相应的普通矩阵的加法和乘 法 有分块矩阵的性质可以知道 这是一个等价的转换 定理由单元结构表示的 每个元素是同维方阵的任意多个嵌套矩阵的运算等于它们所表示的原始矩 阵的运算 证明因为如果任意两个 M 阶的嵌套矩阵 A 和 B 的所有对应的元素都是同维的 n 阶方阵 那么他们 正好可以看成是两个分块矩阵 并且他们对应的分块都是同维的方阵 那么这种分法就满足分块矩阵的加 法和乘法等运算的条件 前一个矩阵列的分法与后一个矩阵的行的分法完全一样 都是对一个 n M 阶的 方阵每隔 n 行画一个横线和每隔 n 列画一个竖线分割而成 且 A 的列块数 B 的行块数 M 因此这两个 矩阵是可以进行加法和乘法运算的 根据分块矩阵的性质可知 分块矩阵的运算等于它们所表示的原始矩 阵的运算 故定理得证 第 1 期沈来信等 无穷嵌套矩阵运算在 Matlab 中的实现65 根据上述定理 转换后的两个普通矩阵的加法和乘法就是原分块矩阵的加法和乘法的逆分块 把分块矩 阵中的分块取消 得到的普通矩阵 如公式 7 h 是 1 M 阶矩阵 其中每个元素都是 I I G u v 和公 式 5 类似 是个 M 阶的下三角阵 用 IG n M G 函数求出 U 1 u v 是个 M 阶的对角阵 对角线上每个 元素是个 n 阶的矩阵 用 U1 n M U 求出 f u v h I G u v 1U 1 u v I R u v 7 f cell2mat h n M cell2mat IG n M G cell2mat U1 n M U cell2mat IR n M R 用 来实现 7 其中 h n M 函数实现了一个 1 M 阶的行向量 每个元素都是个单元矩阵 cell2mat 是把相应 的单元结构即嵌套矩阵转变成普通的矩阵 注意此时的矩阵的维数变成了 n M 维 其中 n 为原单元结构 每一个元素的维数 M 为整个单元结构的维数 把结果变成单元结构类型即 f1 mat2cell f n n M 表示把 n M 阶方阵转换成相应的每个元素为 n 阶的方阵的 M 阶嵌套矩阵 4实验结果与分析 我们在程序分别采用 10 阶 10 阶 20 阶的嵌套矩阵 其中每个元素又分别是 4 阶 10 阶 10 阶的方阵 这样公式 7 就分别对 3 个 40 阶 100 阶 200 阶的普通矩阵做乘法运算 并用 clock 结合 etime 函数来测试 运行时间 最后系统运行时间为 4 875 s 5 657 s 24 859 s 当 n 大于 20 阶时 无穷嵌套矩阵 IR 和 IG 中 的元素就趋于 0 了 公式 7 具有明显的收敛特征 所以如果无穷嵌套矩阵具有收敛性质 那么无穷嵌套矩阵 的运算就可以在有限步内得到求解 这也充分体现了 Matlab 超强的大矩阵和嵌套矩阵的运算能力 如表 1 所示 表 1 不同阶数的嵌套矩阵运行用时及误差 嵌套外矩阵内矩阵普通矩阵参与运算阶数用时最大误差 10 阶 4 阶 40 阶 10 40 40 404 87510 6 10 阶 10 阶 100 阶 10 100 100 1005 65710 6 20 阶 10 阶 200 阶 20 200 200 20024 85910 7 实验结果表明 运用 Matlab 中单元结构和逆分块矩阵的思想可以有效地解决无穷大嵌套矩阵的表示和 运算的问题 参考文献 1 LI Q ZHAO Y Block structured fluid queues driven by QBD processes J Stochastic Analysis and Applications 2005 23 1 087 1 112 2 任志刚 关于 Z 矩阵与 M 矩阵广义 Perron 补矩阵的研究 D 成都 电子科技大学 2007 3 黄永丽 李晔 三维重构中层次性断层数据的处理技术 J 计算机工程与设计 2008 24 6 407 6 409 The Implementation of Infinity Nested Matrices Operation in Matlab SHEN Lai xin1 REN Yong2 YANG Fan1 1 College of Information Engineering Huangshan University 245021 Huangshan Anhui China 2 College of Mathematics and Computer Science Anhui Normal University 241000 Wuhu Anhui China Abstract The expression and effective operations of infinity nested matrix is always difficult the former can use the function cell to implement in Matlab but cell based structure does not support its direct operations we use the idea of inverse block matrix to implement the addition and multiplication of nes