matlab可计算多大阶大型线性稀疏矩阵Ax=b

newdiver

算过十万阶的 Ax=b求解器是  x=A\b;  速度特别慢，一看内存是满的，A 和b本身没占多少内存啊，但是8G的电脑配置内存一下子满了，想知道matlab可计算多大阶大型线性稀疏矩阵，如果是非常大的稀疏矩阵，matlab是否不适用了 ，是不是要换其他的软件来计算呢
求高人指点

messenger

计算所用的内存，与矩阵所占的内存是两码事。稍微有点数值计算概念的人，也不会问这个问题。

newdiver

不好意思让斑竹见笑了，不过是问题总归要问的吧不论大小，那有没有解决内存不足的方法呢，比如想算20万的，窗口显示Out of memory 了，除了增加内存外还有更好的办法吗

messenger

最好的办法当然是增加内存了。
此外，你还可以试试矩阵分块，http://forum.simwe.com/thread-962367-1-1.html。

newdiver

矩阵分块是不是只针对图像的呢，我只是想求出 Ax=b 中的x，A是稀疏对称n*n矩阵，如果把A分块了，怎么去计算呢，谢谢指教！

pasuka

不好意思让斑竹见笑了，不过是问题总归要问的吧不论大小，那有没有解决内存不足的方法呢，比如想算20万的，窗口显示Out of memory 了，除了增加内存外还有更好的办法吗
newdiver 发表于 2011-4-14 19:08

用symamd或symrcm把矩阵重排一下，也许能够节约一点点内存
加内存肯定是最方便快捷的办法，实在不行的话，只有换C/Fortran，再不行动话，要么修改算法，把一部分计算数据转移到硬盘上来节约内存

newdiver

6# pasuka
谢谢指教，考虑C语言混编，不知道会不会好一点

pasuka

6# pasuka
谢谢指教，考虑C语言混编，不知道会不会好一点
newdiver 发表于 2011-4-15 19:30

安装ivf专业版，里面自带pardiso，主流稀疏矩阵求解器，可以用C/Fortran混编
如果还是内存不足的话，还是加内存吧

rocwoods

20万阶的稀疏矩阵，看稀疏程度如何，你可以生成一下看看占多少内存，如果还有一半以上物理内存空余的话，MATLAB搞定一般没问题。这个不是哪个语言的问题，而是你的算法和如何使用MATLAB的问题。另外虽然blockproc是用来处理大图像的，但大图像也是矩阵，blockproc当然可以处理矩阵。

newdiver

9# rocwoods 谢谢版主指教
blockproc可以处理矩阵，只是不知道blockproc如何应用到 Ax=b 的计算中，版主能不能解释具体一点呢？
另外，我看了一下20万阶的稀疏矩阵占内存不超过总物理内存的百分之十，只是一运行 x=A\b 的时候，内存消耗突然增加到满格，然后出现内存不足的信号，我怀疑是不是x=A\b这个命令太耗内存，应该不可能是中病毒吧，因为其他时候都是正常的

penultimate

9# rocwoods 谢谢版主指教
blockproc可以处理矩阵，只是不知道blockproc如何应用到 Ax=b 的计算中，版主能不能解释具体一点呢？
另外，我看了一下20万阶的稀疏矩阵占内存不超过总物理内存的百分之十，只是一运行  ...
newdiver 发表于 2011-4-18 20:00

A 确实太大了，如果不考虑内/外存交换的 LDLT 类求解器，或 sparse matrix，基本上无法求解

penultimate

9# rocwoods 谢谢版主指教
blockproc可以处理矩阵，只是不知道blockproc如何应用到 Ax=b 的计算中，版主能不能解释具体一点呢？
另外，我看了一下20万阶的稀疏矩阵占内存不超过总物理内存的百分之十，只是一运行  ...
newdiver 发表于 2011-4-18 20:00

A 确实太大了，如果不考虑内/外存交换的 LDLT 类求解器，或 sparse matrix，基本上无法求解

rocwoods

10# newdiver
从最一般情况考虑，假设A非常大，存放在硬盘上。
先利用blockproc分块A，假设A = [B11 B12;B21 B22];（其中B11,B22为方子阵）分块后的矩阵也放在硬盘上。根据下面公式求出A的逆矩阵A^(-1) = [C11 C12;C21 C22];
C22 = (B22-B21*B11^(-1)*B12)^(-1);
C12 = -B11^(-1)*B12*C22;
C21 = -C22*B21*B11^(-1);
C22 = B11^(-1) - C12*B21*B11^(-1);
如果某个A的对角子阵Bii还大的超过内存话，在对其分块求出逆矩阵，以此类推。这样可以把A的逆矩阵求出来。求出来后再和b相乘就得到解了。

假设求出来的C11 C12 C21 C22都很大，放在硬盘上，只要其文件名命名规范，就可以利用blockproc，来求其以b的乘积，得到解。
特别值得说明的是，这种方法可以处理任意大的矩阵，不受制于内存。只要你有足够的时间和硬盘存储分块的矩阵。

pasuka

10# newdiver
从最一般情况考虑，假设A非常大，存放在硬盘上。
先利用blockproc分块A，假设A = ;（其中B11,B22为方子阵）分块后的矩阵也放在硬盘上。根据下面公式求出A的逆矩阵A^(-1) = [C11 C12;C21 C22];
C22  ...
rocwoods 发表于 2011-4-19 17:17

一般大型矩阵是不直接求逆的，此外，20w的实对称稀疏矩阵不能算很大，自己编程的话Fortran/C/C++都能够应付

rocwoods

15# pasuka
上面只是给了楼主一个方法，毕竟从楼主发帖初衷推断，还是想在MATLAB下搞定的。大矩阵的话，也有近似求逆的方法，就看楼主对问题的精度需要了。还是那个观点，这个不是用哪个语言的问题，而是能用MATLAB处理的话干嘛非得用C/C++/Fortran，而且对大多数人来讲用MATLAB处理矩阵问题，效率未必比自己用C/C++吭哧半天写出来的要低。
以前项目上处理遥感影像数据，动不动就是几G几十G的原始影像，重采样一幅大的影像，调用blockproc封装成的dll处理比自己写C++处理，要快上好多倍。当然这和开发人员的C++水平有关系，但不可能每个人都是C/C++专家。

iomega

最近我也做了些求解大型线性稀疏矩阵AX=b的工作，我的体会是解AX=b无非是直接法和迭代法两种。
对于直接法，如果是对称稀疏矩阵，可以先用LDL^T,Cholesky或QR进行matrix factorization, 比如： A = QR, 然后用x = R\Q\b来求解。如果是一般nxn稀疏矩阵，可用LU factorization, 建议用 [ L, U, P, Q, R] = lu(A), 然后 x = Q * ( U \ (L \ (P * (R \b))));  计算速度取决于稀疏矩阵的number of non-zeros terms。如果用32位系统的话， 阶数大的话，会出现out of memory, 那只有上64位系统了。

给个例子，我的64位Win 7, 36G DDR3 memory,  12 core/24threads, 10万阶稀疏矩阵Ax=b,  Matlab用LU factorization求解，大约花30秒左右, 用的是自带的UMFPACK.

如果内存不够， 可以用out-of-core的求解器。 我现在用的是GNU的MUMPS+GotoBlas+Metis ordering。 在Windows 下，用microsoft C 和 g95 fortran 编译生成lib文件用C code调用, 对于上面同样的问题，大约12秒左右.

我现在算的问题，阶数大概40万， matlab用500秒， MUMPS C code 用120秒左右。



 

newdiver

16# iomega 谢谢指点，非常受用，揣测了许久，求解大型稀疏矩阵两种方法 1.用matlab 自带求解器求
  2. 内存不够的话，根据GNU的MUMPS+GotoBlas+Metis ordering，用c语言调用fortran和c编译的lib文件计算，（这个没有在matlab里调用吧）

本人对计算机是个外行，不知是否正确，还望版主纠正

另外 out of core 求解器  还有 GNU的MUMPS+GotoBlas+Metis ordering，看的一头雾水，能不能具体解释一下呢
非常感谢

iomega

你的理解正确。 这里是我在32 位winodws下用GNU 的gcc 和g95编译的MUMPS lib 文件， 并且有个例子c_example.c.     comp.bat 是用来在DOS环境下用命令行调用cl.exe （visual C++)编译和link c_example.c。所需要的lib文件在lib目录下面。

如果你安装了visual C, 只需解压rar到给定目录下，然后运行comp, 生成c_example.exe.
你看懂的话，自己改写， 解50万阶的AX=b应该没有问题的。


17# newdiver

chinasnowlee

1# newdiver
请问A是不是m*n的矩阵  并且m>n  
也就是所求的是一个矛盾方程组!
有三种方法吧  
左除
lsqnonneg
lsqlin

我也有计算!