mldivide for sparse matrix runs out of memory
in [Matlab]
|
From: Pascal Schulthess on 25 Feb 2010 10:57 Hi everybody, I've read almost all threads here on sparse matrices, lu, qr and so on and it really helped me move forward with my problem. Everyone says, that if you have to solve something like A*x=b, use the backslash operator. Unfortunately, this results in an out of memory error for my matrices. A has full rank and is a (12729 x 9303) sparse matrix with a density of 3.6847e-04. Furthermore, it might be important to mention that A = [ I ; A0 ] where A0 is a (3426 x 9303) matrix with a density of 0.0011. b is rank deficient and has the dimensions (12729 x 11395) sparse matrix with a density of 4.4472e-04. Now I know that x has to be a (9303 x 11395) matrix. But if I perform: x = A\b my 64bit 2x2.8 Ghz Quad-Core Intel Xeon (16 GB RAM) runs out of memory. Any suggestions or ideas how to tackle this? Thanks in advance, Pascal
From: Matt J on 25 Feb 2010 11:49 "Pascal Schulthess" <news(a)pascalschulthess.de> wrote in message <hm66ke$hr7$1(a)fred.mathworks.com>... > Hi everybody, > > I've read almost all threads here on sparse matrices, lu, qr and so on and it really helped me move forward with my problem. Everyone says, that if you have to solve something like A*x=b, use the backslash operator. > Unfortunately, this results in an out of memory error for my matrices. > > A has full rank and is a (12729 x 9303) sparse matrix with a density of 3.6847e-04. Furthermore, it might be important to mention that A = [ I ; A0 ] where A0 is a (3426 x 9303) matrix with a density of 0.0011. > > b is rank deficient and has the dimensions (12729 x 11395) sparse matrix with a density of 4.4472e-04. > > Now I know that x has to be a (9303 x 11395) matrix. > > But if I perform: x = A\b my 64bit 2x2.8 Ghz Quad-Core Intel Xeon (16 GB RAM) runs out of memory. ====== The problem you have to get around is that x will be stored as type sparse, but inherently will not be sparse. So, it will consume 3 times as much memory as if you stored it full. You would really like for A\b to return a non-sparse result, however, there is no way to do that without converting A or b to full, which you don't want to do. The only solution I can think of is to partition the columns of b into ten or so smaller matrices and solve each as a sub-problem. Each time you solve a sub-problem, convert the resulting x to full before proceeding to the next one, e.g. bchunks=mat2cell(b,12729,(1:11)*1085); Xcell=cellfun(@(b), full(A\b), bchunks,'uni',0); clear A b x=[Xcell{:}];
From: Matt J on 25 Feb 2010 11:59 "Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message <hm69m0$7tt$1(a)fred.mathworks.com>... > Xcell=cellfun(@(b), full(A\b), bchunks,'uni',0); ===== In fact, you could save even more memory if you also converted the results to singles Xcell=cellfun(@(b), single(full(A\b)), bchunks,'uni',0);
From: Pascal Schulthess on 26 Feb 2010 03:41 "Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message <hm6a8p$h4q$1(a)fred.mathworks.com>... > "Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message <hm69m0$7tt$1(a)fred.mathworks.com>... > > > Xcell=cellfun(@(b), full(A\b), bchunks,'uni',0); > ===== > > In fact, you could save even more memory if you also converted the results to singles > > Xcell=cellfun(@(b), single(full(A\b)), bchunks,'uni',0); Hi Matt, this works great, although I hoped to complete this task a little faster ;) (it takes about 50 minutes) Another idea I had was to benefit from the form of of A since there's identity matrix in there. Perhaps something like inv(A) = inv([ I ; A0 ]) = [ I , inv(A0) ]
From: Pascal Schulthess on 26 Feb 2010 05:25
"Pascal Schulthess" <news(a)pascalschulthess.de> wrote in message <hm81f1$sp5$1(a)fred.mathworks.com>... > "Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message <hm6a8p$h4q$1(a)fred.mathworks.com>... > > "Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message <hm69m0$7tt$1(a)fred.mathworks.com>... > > > > > Xcell=cellfun(@(b), full(A\b), bchunks,'uni',0); > > ===== > > > > In fact, you could save even more memory if you also converted the results to singles > > > > Xcell=cellfun(@(b), single(full(A\b)), bchunks,'uni',0); > > Hi Matt, > > this works great, although I hoped to complete this task a little faster ;) (it takes about 50 minutes) > > Another idea I had was to benefit from the form of of A since there's identity matrix in there. Perhaps something like > > inv(A) = inv([ I ; A0 ]) = [ I , inv(A0) ] Ok, I found something called block matrix pseudoinverse at the wikipedia (http://en.wikipedia.org/wiki/Block_matrix_pseudoinverse). Unfortunately, after simplifications due to identity matrix in A, I still have to inverse A0 a couple of times, and since A0 is still (3426 x 9303), this is hard to do. Thanks again for you help Matt. |