做矩阵分块解析运算（证明也解决了，软件是不存在的，更新在3页）
本帖最后由 orionsnow 于 2011-1-12 19:55 编辑有什么软件可以做矩阵分块运算的？
mathematica 不行
maple？
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Dear Sean,
Thank you very much for the information.
the infor you gaved me is matrix decomposition.
but the one i want to use is block wise matrix decomposition.
maybe i am not make my question not clear
in detail.
here i have a big matrix E, which have 4 blocks
[A,B
C,D]
latex code
E = $\displaystyle \left[\begin{array}{cc} A & B \\ C & D
\end{array}\right]^{-1}$
A, B, C, D are full rank square matrix.
so the block wise inverse E^-1is
[ S_d, S_d^-1 .B .D^-1,
..., ..............]
latex code:
$\displaystyle \left[\begin{array}{cc} S^{-1}_{D} & -S^{-1}_{D}BD^{-1}
\\ -D^{-1}CS^{-1}_{D} & D^{-1} + D^{-1}CS^{-1}_{D}BD^{-1}
\end{array}\right]$
in short,
i what repeat the precedure that listed here by mathematica7. that i
give ABCD, then mathematica inverse the E as a matrix of ABCD.
https://ccrma.stanford.edu/~jos/lattice/Block_matrix_decompositions.html
thank you again for all the help
Best Regards
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
----------------------------------------------------------------------------
Hello
If a function does not exist in Mathematica, you will have write it using
the already defined functions in Mathematica. As I stated in the previous
email, I do not believe the function you want is a predefined function in
Mathematica. However, defining this function in Mathematica should not be
too difficult.
I have attached a notebook showing how to access the different quadrants of
a matrix and recombine them. This should provide a useful starting point
for writing your function.
If you are not familiar with programming in Mathematica, please read the
documentation on lists and functions:
http://reference.wolfram.com/mathematica/howto/WorkWithLists.html
http://reference.wolfram.com/mathematica/tutorial/DefiningFunctions.html
If you have any specific questions about how to use the documented
functions in Mathematica while making this program, please let me know.
Sincerely,
Sean Clarke
Technical Support 能算大型稀疏矩阵的软件就肯定能做分块，你懂的 matlab 本帖最后由 orionsnow 于 2010-10-25 17:27 编辑
matlab
johndoe 发表于 2010-10-25 15:44 http://www.dolc.de/forum/images/common/back.gif
谢了，刚才去网上找了下，看到有个 matrix2cell 命令。
可以把一个矩阵分块。 估计还有其他更多命令。 不知道你有没有其他的例子。
不过现在没有matlab 的lience，回头要去找别的同事试验一下。
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
我的问题和下边这个一模一样，我正在看。
http://www.mapleprimes.com/questions/89722-Calculations-With-Matrix-In-Matrix
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
已解决
matlab矩阵分块与把分块矩阵还原
悬赏分：200 - 解决时间：2008-8-10 17:13
我有一个256*64的矩阵，我要把它分成16*16的矩阵，接着要把这个分好了块的矩阵又可以还原为原256*64的矩阵，请问matlab中该如何编写这个m文件~~谢谢~
提问者： zzyzzm - 三级
最佳答案
A=rand(256,64);
%将A分块
B=mat2cell(A,ones(256/16,1)*16,ones(64/16,1)*16);
%B{i,j}就是所要的分块矩阵
%将分块矩阵合并
C=cell2mat(B)
%C就是合并好的矩阵，即C=A
9 i am not using matlab for about five years. I think your problem is definitely within matlab's shot. Just try /help cell2mat or whatsoever. But I could also remember that cell2mat has to do with the transforming string into numerical values so pls be careful for that.
but honestly, do you really need any matrix with 256*64 dimensions? any matrix operation will be killed with curse of dimensionality.
sorry, i cannot put in any chinese on this computer and has to conclude this way. i am not using matlab for about five years. I think your problem is definitely within matlab's shot. ...
johndoe 发表于 2010-10-25 17:56 http://www.dolc.de/forum/images/common/back.gif
Nein, cell2num sollte den Befehl sein, mit dem eine Umsetzung von String nach numerical erfolgt. 本帖最后由 orionsnow 于 2010-10-25 18:38 编辑
谢谢楼上两位的回贴， matlab 我也是以前写毕业论文的时候用过，主要是做的模拟。
到现在也有5，6年没有碰了，不过原理还是记得。关键知道能解决的包和命令上的名字，回头去google 就是了。
1， 256 *64 是百度上的例子
我的问题是4个或者9 个 最多不超过 1000的满秩非稀松阵，不过在特殊情况下是只有2,3个常数参数的， 比如
diag(1000)+ 0。5 这种的。
然后这4 个阵组成一个2000 的阵。 要求它的逆。目前就在先解决这个简单特例。 因为后边还有其他的运算，觉得不是一下就能行的。 干脆让老板开个项目算了。 要是能保留符号不进行数值结算就最好了，这样可以寻找机会和后边的公式消掉。
2 德语我没有看懂，如果输入不方便的话， 打英语吧。 谢谢楼上两位的回贴， matlab 我也是以前写毕业论文的时候用过，主要是做的模拟。
到现在也有5，6年没有 ...
orionsnow 发表于 2010-10-25 18:37 http://dolc.de/forum/images/common/back.gif
我想你的问题应该可以这样解决
AA = ;
其中AA1是左上角的矩阵，AA2为右上角，BB1为左下角，BB2为右下角，
现在机器上没有matlab了，你可以用小点的矩阵先试试看 本帖最后由 orionsnow 于 2010-10-25 20:44 编辑
我想你的问题应该可以这样解决
AA = ;
其中AA1是左上角的矩阵，AA2为右上 ...
johndoe 发表于 2010-10-25 20:30 http://www.dolc.de/forum/images/common/back.gif
我机器上也没有matlab 呢。。。。。。。。。
然后aa 求逆的符号计算是可行的？
我给同事打电话了，他去做了，先搞个4，4 矩阵，过两天看结果。 我机器上也没有matlab 呢。。。。。。。。。
然后aa 求逆的符号计算是可行的？
我给同事打电话了 ...
orionsnow 发表于 2010-10-25 20:37 http://dolc.de/forum/images/common/back.gif
好啊，事成了刷我们一下 cell是原胞，貌似和分块无关 x2y型命令用于数据类型转换 本帖最后由 aileute 于 2010-10-28 21:36 编辑
分块矩阵不就是种表示方法吗？如果你用到的矩阵维数很高，有很多零，这种矩阵叫稀疏矩阵，稀疏矩阵的存储方式和一般矩阵不一样。因为有很多零，LU，求逆（当然实际应用中没人会去真的先求逆在求解的），有对应的稀疏矩阵算法。因为借此可以通过特殊的算下标的方法来节省很多不必要的运算。
至于lz的分块，把对应的要分的部分提取出来，然后用一般的算法算就可以了，根本不需要什么特殊的函数。你到底要干什么，说具体一点，或许大家会有更多的建议！ 本帖最后由 orionsnow 于 2010-10-28 23:21 编辑
分块矩阵不就是种表示方法吗？如果你用到的矩阵维数很高，有很多零，这种矩阵叫稀疏矩阵，稀疏矩阵的存储方 ...
aileute 发表于 2010-10-28 21:34 http://www.dolc.de/forum/images/common/back.gif
谢谢建议， 不过我想你还是仔细看看我前边的贴子吧，我说的很清楚了。还给了个链接，你可以看那个哈弗大学的例子，他已经给出了分块2阶方矩阵的逆的推导方法.
我的目标矩阵不是稀疏矩阵。 是满秩的相关系数矩阵。
比如 2000 by 2000 的相关矩阵，有100 x 100 个小结构。
我以前课本上看到有很多不错的数值算法。不过我们现在对数值算法不感兴趣，我们在利用矩阵的特殊结构求解析解。
4 x 4 的
R^-1=（1，0.5，0.25，0.25
0.5， 1，0.25，0.25，
0.25，0.25，1，0.25，
0.25，0.25，0.25，1
）^-1
或者
（1，0.5，0.5，0.5
0.5， 1，0.5，0.5，
0.5，0.5，1，0.5，
0.5，0.5，0.5，1
） ^-1
0.5 和0.25 只是举例子。 在实际应用的时候可以取任意定值，或者从实验中直接读取。
样本很大的时候可能会假设一个分布，不过现在还没有到哪一步。 如果数值随机的话，解析解估计很难存在了，最多只能给一个估计。
然后求解方程组
(ZR^-1Z' + lambdc )c + Z W' d =ZR^-1y
WZ' c + (WR^-1W' + lambdd) d=WR^-1 y
这时候要对左边四个块再求一次逆。
here,R is ij by ij , Z isi by ijblock wise diagnal,W is ij by ijidentity matrix, lambdd lambdc is known constant.
y is 1 by ij.c is 1 by i,d is ij by ij.
ifj =1,then all the 4 blocks are squared. else not.
我今天已经用我那可怜的高中手解出来. 然后推了一个通式出来，基本上可以避免超大矩阵求逆了。今天交给老板检查去了，他也是建议我看看线性代数的书。
不过我现在特别不喜欢看书，一看书就打瞌睡。
matlab 那边还没有消息，估计最近开学新生多，我同事泡小姑娘去了。
这个问题基本到这里了，文章已经二审了，就等着个推导，等发出来了我把appendix 传上来。 本帖最后由 aileute 于 2010-10-28 23:47 编辑
回复 14# orionsnow
这样才说清楚了，数值算法要具体问题具体分析的。你是因为2000x2000的矩阵太大了，所以想采用分块的方法来求解。你给出来的网页里的算法是递归的方法，很明显采用的是divide and conquer。不过介绍中没有提及计算复杂度的问题。这个领域不了解，不过我以前只知道对于矩阵乘法O(n^2)，divide and conquer有明显的降低计算复杂度的作用有的文章说可以达到O(n^1.5)，但是对于解方程O(n^3)没有看到过。当然这个领域不熟悉，可能也有改善。
不过大矩阵求逆始终要非常多时间。我的建议就是用matlab直接算，估计3-5个小时内就可以算出来！而且精度病态问题都有保证！
p.s.羡慕一下你老板，你推导个公式都可以让老板帮你检查，挺幸福的！
还有就是一个概念上的问题，如果我理解正确的话，满秩保证矩阵可逆。稀疏矩阵有很多零，但是工程中大部分也是满秩，否则解空间就塌陷了，有多余的线性相关变量在里面了！说明排的方程有冗余。 回复orionsnow
这样才说清楚了，数值算法要具体问题具体分析的。你是因为2000x2000的矩阵太大了，所 ...
aileute 发表于 2010-10-28 23:43 http://www.dolc.de/forum/images/common/back.gif
1
对，满秩是可逆的必要条件。
我说满秩的意思是保证那个不稀松的矩阵一定可逆。 因为一个全满阵看上去张牙舞爪，但是超大的矩阵很可能线性相关行从而导致不可逆。
其实正定随机矩阵已经一定是可逆的了。 博版能人很多，不说的详细了很多人后边会帮你问清楚的。
2
老板一般不亲自上的。 这次要亲自检查，是因为这个学生的文章被reviewer 退回来两次了。就是因为这个公式的证明问题。 而且是小老板先检查，然后大老板再复查。
3
关于时间问题， 数值解是不可行的， 5 个小时，太多5分钟我都觉得多。
因为这个公式是一个打分公式，后边要给simulation 用。 所以必须在几秒种内出结果。否则simulation 估计要化几个世纪。 我更正上面的说法，2000x2000的矩阵求逆用LU（因为先前没有测过，大概估计，发现自己水平还是不够，和真实值差的非常远。因为这些算法在我家电脑上都有，所以我试试还是很容易的。）大概200来秒。测试用的机器是很普通的台机，windows下，AMD 2.6G， 2G内存。LU算法实现来自gsl库。个人猜测matlab还要快！
p.s.还有ls在3中的提法，我的观点是错误的。请问矩阵求逆的解析解是什么？ 本帖最后由 orionsnow 于 2010-10-29 22:44 编辑
我更正上面的说法，2000x2000的矩阵求逆用LU（因为先前没有测过，大概估计，发现自己水平还是不够，和真实值 ...
aileute 发表于 2010-10-29 21:06 http://www.dolc.de/forum/images/common/back.gif
1
c 语言应该是最快的了吧？除了汇编，
除非matlab 还做了其他优化了。
2
解析解就是解析解啊。
你问是要问解析解的定义是什么？
不严格定义的话，就是不算数值用符号表示的那种。
比如说，
我上边那个公式里头,在系数是0.5 和 0的情况下。 对任意维矩阵。
c_i= \bar{y}_i *(nk lambdac) / (nk- (1+lambda_d)(lambda_c))
nk 是子块的维数。 \bar{y}_i 是所有和c_i 有关的观测值求平均。
d_ij 还要长点，我记不住了。
其实为什么要求解析解这个问题很难解释清楚， 有个reviewer 也提到线性加速算法了，也是问我们为什么不用。
这个问题基本上到这里可以告以段落了。 到目前小老板没有查出问题来。
数值算法方面，基本上知道gsl 要用200 秒这个信息就够了。 那估计2000 维矩阵相乘也就是几秒的事情。
那个学生不会c，主函数是用R 写的。 但是这些都不是大问题。
不过这么短的时间真的让我很吃惊。 可能在文章里头加上一句。
“数值算法使用LU 分解，
大概只需要3-5分钟就可以求解出R^-1， 然后保存R^-1
然后再用3-5 分钟求解上边的mme 的逆矩阵 然后保存。
就可以在普通pc 上用数值解模拟这个过程只比使用解析解慢 10分钟+每次模拟多出来5，6 个 矩阵乘法 。内存消耗估计大一些。
原来十五 天的模拟 估计最多就是多用一天。
当问题比较复杂而且分块比较多，解析解推导比较麻烦的时候可以考虑用数值算法。
” 2000x2000的两个矩阵想成大约耗时50秒，这个数据我在验证求逆是是否出现病态问题是一并都算了。我觉得你可能想知道。
LU分解以后，就可以直接求解方程了，其他的我不多说了，再说下去我都不好意思了。如果有空建议复习一下线性代数和数值算法，总会有所收获的！ 本帖最后由 orionsnow 于 2010-10-30 11:47 编辑
2000x2000的两个矩阵想成大约耗时50秒，这个数据我在验证求逆是是否出现病态问题是一并都算了。我觉得你可能 ...
aileute 发表于 2010-10-29 23:08 http://www.dolc.de/forum/images/common/back.gif
不过具体怎么样还是要安排实验下才知道。 书本上来的知识只能做参考。
唉，我昨天还给老板发信说可能有重大突破， 幸亏我信写的保守了一些。 我写的是
400 +matrix multiplication x simulation time seconds =？ days
继续说吧，没啥不好意思的，不要轻视基础知识。 就是细节的地方才会出成果。也许我又把什么基础知识漏掉了呢。
不管怎么说， 数值算法 作为解析算法的验证还是有价值的。 虽然速度上比解析算法要慢和麻烦一点。 但是在推导过程上要快很多。
不过说实话，我们又把书本上的基础知识重复验证了一次。 这个事实貌似某个课本上写过，不记得那本了。 唉今天又出新问题了，我发现我昨天忽略掉一个参数，这个参数每次simulation 都变的，这就是说每次 那个逆矩阵都要重新求。 我印象里头我们的simulation 原来是1-5 秒，这样如果算上求逆的200 秒就是205 秒了。
继续求可以符号计算矩阵运算的软件 怪蜀黍的英语真强悍啊 怪蜀黍，不是鄙视你，你懂的太多了，但一样都不精，standford网页上那个东东是国内大一线数的内容
实在搞不定的话，短我吧，matlab程序，自己写的，200欧，保讲懂，要c++的话，400欧 本帖最后由 庄十三爷 于 2010-11-1 20:32 编辑
怪蜀黍，不是鄙视你，你懂的太多了，但一样都不精，standford网页上那个东东是国内大一线数的内容实在搞不定的话，短我吧，matlab程序，自己写的，200欧，保讲懂，要c++的话，400欧
熊猫羊 发表于 2010-11-1 14:39 http://www.dolc.de/forum/images/common/back.gif
400 欧太多了，叔叔还是短我吧， 8euro 88 cent 就可以了 就可以了 庄13爷，您对一个连 矩阵分解 和 矩阵分块 都没搞清的冤大头也太仁慈了吧？ 不少了，你想阿，如果是在柏林踩三轮车给人家帮忙搬家，一个小时十快钱最多了，写个程序怎么也要一天8个小时吧。所以换算成欧元差不多就是这个价。
运气好还可以卖人家个沙发，赚个50 大欧， 不过这种事情不是每次都能碰上的。 再说还有被人家人肉的风险，我可不认识版主，弹压不下去。 本帖最后由 农大 于 2010-11-11 21:44 编辑
2000*2000 用 “解析" 是不合适的
而且所有所谓的数学软件matlab之类，对待这样的矩阵也不会用所谓的”解析“解
lz 要是会c或者c++ 随便找个solver 解一下就行。。。 本帖最后由 orionsnow 于 2011-1-12 19:59 编辑
上来汇报下。 花了一个多月， 两个超大矩阵在特殊情况下的解析解都被我给推导出来了。
然后我老板在网上搜了几天，居然找到了第一个矩阵解法的推导，那人最后把特例朝他的一般推导里头一代入，就得到了我推出来的解析解。 根据这个事实和数值算法的验证，就是前边无数同学提到的， solver（） inverse（） 这类的函数。 确认第一个矩阵的解法正确。 然后老板自己验证了下第二个公式也是对的。
我现在在模仿那人的一般推导推第二个公式。
论文在这里，感兴趣的可以看了以后讨论下， 公式41-46 是主要证明， 结论在第5章，公式50多。
Searle SR (1996) The matrix handling of BLUE and BLUP in the mixed linear
model. Linear Algebra Appl 264:291-311 232
这个人挺牛的， 普通解推导用的矩阵计算， 特解带入用的矩阵分块运算，主要证明就两页纸。
我为了和老板说明我的思路前后写了40页纸。 现在正在好好学习他的矩阵分块算法的推导。
矩阵分块计算的推导软件应该是不存在的， 我查了这个人paper 的被引用情况，没有人提到更深入的解法了。
估计现在计算机技术还没有那么发达，只能求数值解。 如果计算机能做推导，那估计数学家可以改行了。 再顶起来一下， 下一篇文章还是关于这类公式推导的。现在老板都快把这个当产业了，把以前的老算法都拿出来， 设某个参数为特殊值，然后加上点解析化简，就变成新文章了。
页:
[1]