做矩阵分块解析运算（证明也解决了，软件是不存在的，更新在3页）

熊猫羊

cell是原胞，貌似和分块无关

熊猫羊

x2y型命令用于数据类型转换

aileute

分块矩阵不就是种表示方法吗？如果你用到的矩阵维数很高，有很多零，这种矩阵叫稀疏矩阵，稀疏矩阵的存储方式和一般矩阵不一样。因为有很多零，LU，求逆（当然实际应用中没人会去真的先求逆在求解的），有对应的稀疏矩阵算法。因为借此可以通过特殊的算下标的方法来节省很多不必要的运算。

至于lz的分块，把对应的要分的部分提取出来，然后用一般的算法算就可以了，根本不需要什么特殊的函数。你到底要干什么，说具体一点，或许大家会有更多的建议！

orionsnow

分块矩阵不就是种表示方法吗？如果你用到的矩阵维数很高，有很多零，这种矩阵叫稀疏矩阵，稀疏矩阵的存储方 ...
aileute 发表于 2010-10-28 21:34

谢谢建议， 不过我想你还是仔细看看我前边的贴子吧，我说的很清楚了。还给了个链接，你可以看那个哈弗大学的例子，他已经给出了分块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^-1  y
WZ'                           c + (WR^-1W' + lambdd) d  =WR^-1 y

这时候要对左边四个块再求一次逆。

here,  R is ij by ij ,   Z is  i by ij  block wise diagnal,  W is   ij by ij  identity matrix, lambdd lambdc is known constant.
y is 1 by ij.  c is 1 by i,  d is ij by ij.

if  j =1,  then all the 4 blocks are squared.   else not.


我今天已经用我那可怜的高中手解出来. 然后推了一个通式出来，基本上可以避免超大矩阵求逆了。今天交给老板检查去了，他也是建议我看看线性代数的书。

不过我现在特别不喜欢看书，一看书就打瞌睡。

matlab 那边还没有消息，估计最近开学新生多，我同事泡小姑娘去了。

这个问题基本到这里了，文章已经二审了，就等着个推导，等发出来了我把appendix 传上来。

aileute

回复 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

回复  orionsnow


这样才说清楚了，数值算法要具体问题具体分析的。你是因为2000x2000的矩阵太大了，所 ...
aileute 发表于 2010-10-28 23:43

1
对，满秩是可逆的必要条件。

我说满秩的意思是保证那个不稀松的矩阵一定可逆。 因为一个全满阵看上去张牙舞爪，但是超大的矩阵很可能线性相关行从而导致不可逆。

其实正定随机矩阵已经一定是可逆的了。 博版能人很多，不说的详细了很多人后边会帮你问清楚的。

2
老板一般不亲自上的。 这次要亲自检查，是因为这个学生的文章被reviewer 退回来两次了。就是因为这个公式的证明问题。 而且是小老板先检查，然后大老板再复查。

3
关于时间问题， 数值解是不可行的， 5 个小时，太多5分钟我都觉得多。

因为这个公式是一个打分公式，后边要给simulation 用。 所以必须在几秒种内出结果。否则simulation 估计要化几个世纪。

aileute

我更正上面的说法，2000x2000的矩阵求逆用LU（因为先前没有测过，大概估计，发现自己水平还是不够，和真实值差的非常远。因为这些算法在我家电脑上都有，所以我试试还是很容易的。）大概200来秒。测试用的机器是很普通的台机，windows下，AMD 2.6G， 2G内存。LU算法实现来自gsl库。个人猜测matlab还要快！
p.s.还有ls在3中的提法，我的观点是错误的。请问矩阵求逆的解析解是什么？

orionsnow

我更正上面的说法，2000x2000的矩阵求逆用LU（因为先前没有测过，大概估计，发现自己水平还是不够，和真实值 ...
aileute 发表于 2010-10-29 21:06

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 个 矩阵乘法 。  内存消耗估计大一些。

原来十五 天的模拟 估计最多就是多用一天。

当问题比较复杂而且分块比较多，解析解推导比较麻烦的时候可以考虑用数值算法。
”

aileute

2000x2000的两个矩阵想成大约耗时50秒，这个数据我在验证求逆是是否出现病态问题是一并都算了。我觉得你可能想知道。
LU分解以后，就可以直接求解方程了，其他的我不多说了，再说下去我都不好意思了。如果有空建议复习一下线性代数和数值算法，总会有所收获的！

orionsnow

2000x2000的两个矩阵想成大约耗时50秒，这个数据我在验证求逆是是否出现病态问题是一并都算了。我觉得你可能 ...
aileute 发表于 2010-10-29 23:08

不过具体怎么样还是要安排实验下才知道。 书本上来的知识只能做参考。

唉，我昨天还给老板发信说可能有重大突破， 幸亏我信写的保守了一些。 我写的是

400 +  matrix multiplication x simulation time   seconds =？ days

继续说吧，没啥不好意思的，不要轻视基础知识。 就是细节的地方才会出成果。也许我又把什么基础知识漏掉了呢。


不管怎么说， 数值算法 作为解析算法的验证还是有价值的。 虽然速度上比解析算法要慢和麻烦一点。 但是在推导过程上要快很多。

不过说实话，我们又把书本上的基础知识重复验证了一次。 这个事实貌似某个课本上写过，不记得那本了。