smallest singular value of large, tall, sparse matrix
in [Matlab]
|
Prev: NI USB-6501
Next: Multiple Tabs with dlmwrite
From: Bruno Luong on 28 Jul 2010 02:47 "Jakob Heide Jørgensen" <jhjspam(a)gmail.com> wrote in message <i2nqdt$r7l$1(a)fred.mathworks.com>... > > The documentation of mldivide says that the rank is determined from a QR factorization with pivoting. Using R = qr(A,0) I get the same segmentation fault. To have pivoting, I believe you must call with three outputs (column permutation, available for sparse matrix since v 2008, not sure 'a' or 'b'). I believe the Q matrix is still huge regardless and might result to crash of the memory is not large enough. Another way is doing incomplete Cholesky factorization on A'*A (two-output call), but that is also memory eater. The SVDs and EIGs (when calling correctly) use deflation and invert iteration methods, and are surely less demanding. Bruno
From: Jakob Heide Jørgensen on 28 Jul 2010 03:22 "Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <i2ojpb$iur$1(a)fred.mathworks.com>... > "Jakob Heide Jørgensen" <jhjspam(a)gmail.com> wrote in message <i2nqdt$r7l$1(a)fred.mathworks.com>... > > > > > The documentation of mldivide says that the rank is determined from a QR factorization with pivoting. Using R = qr(A,0) I get the same segmentation fault. > > To have pivoting, I believe you must call with three outputs (column permutation, available for sparse matrix since v 2008, not sure 'a' or 'b'). I believe the Q matrix is still huge regardless and might result to crash of the memory is not large enough. > > Another way is doing incomplete Cholesky factorization on A'*A (two-output call), but that is also memory eater. > > The SVDs and EIGs (when calling correctly) use deflation and invert iteration methods, and are surely less demanding. > > Bruno s = eigs(A'*A,1,0) crashes after about half an hour with the error message ??? Error using ==> lu Sparse lu with 4 outputs (UMFPACK) failed Error in ==> eigs>LUfactorAminusSigmaB at 1079 [L,U,P,Q] = lu(AsB); Error in ==> eigs at 132 [L,U,P,permAsB] = LUfactorAminusSigmaB; I have s = svds(A,1,0) running on the 7'th hour now. When it finishes, I will report the result here.
From: Bruno Luong on 28 Jul 2010 08:41 "Jakob Heide Jørgensen" <jhjspam(a)gmail.com> wrote in message <i2olqu$hj$1(a)fred.mathworks.com>... > > s = eigs(A'*A,1,0) crashes after about half an hour with the error message Not sure if it helps, but it's probably more clever to code A'*A as function form so as to preserve sparsity rather than fill large matrix A'*A , example: % Small example data A=sparse(rand(10,8)) Afun=@(x) A\(x'/A)'; s = sqrt(eigs(Afun,size(A,2),1,'sm')) % Check min(svd(full(A))) Bruno
From: Pat Quillen on 28 Jul 2010 10:04 "Jakob Heide Jørgensen" <jhjspam(a)gmail.com> wrote in message <i2olqu$hj$1(a)fred.mathworks.com>... > "Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <i2ojpb$iur$1(a)fred.mathworks.com>... > > "Jakob Heide Jørgensen" <jhjspam(a)gmail.com> wrote in message <i2nqdt$r7l$1(a)fred.mathworks.com>... > > > > > > > > The documentation of mldivide says that the rank is determined from a QR factorization with pivoting. Using R = qr(A,0) I get the same segmentation fault. > > > > To have pivoting, I believe you must call with three outputs (column permutation, available for sparse matrix since v 2008, not sure 'a' or 'b'). I believe the Q matrix is still huge regardless and might result to crash of the memory is not large enough. > > > > Another way is doing incomplete Cholesky factorization on A'*A (two-output call), but that is also memory eater. > > > > The SVDs and EIGs (when calling correctly) use deflation and invert iteration methods, and are surely less demanding. > > > > Bruno > > > s = eigs(A'*A,1,0) crashes after about half an hour with the error message > > ??? Error using ==> lu > Sparse lu with 4 outputs (UMFPACK) failed > > Error in ==> eigs>LUfactorAminusSigmaB at 1079 > [L,U,P,Q] = lu(AsB); > > Error in ==> eigs at 132 > [L,U,P,permAsB] = LUfactorAminusSigmaB; > > I have s = svds(A,1,0) running on the 7'th hour now. When it finishes, I will report the result here. Hello Jakob. You can also try op = @(x) A'*(A*x); lam = eigs(op, size(A,2), 1, 'sa', struct('issym',true,'isreal',true)) which will use much less memory than requesting 'sm' but will also most likely converge much more slowly than using 'sm' or 0 as the target. To improve convergence you could also try increasing maxit and p. See help eigs for more details, but for example, try opts.issym = true; opts.isreal = true; opts.p = 100; opts.maxit = 5000; opts.disp = 1; % to turn on display messages and watch the eigenvalue converge. lam = eigs(op, size(A,2), 1, 'sa', opts); Be aware that eigs will return 0 if convergence did not occur and I don't think it will warn (something to fix!), which is why I recommend turning on the display (which was turned off by default in R2009a, I think). Pat.
From: Steven_Lord on 28 Jul 2010 11:29
"Jakob Heide Jørgensen" <jhjspam(a)gmail.com> wrote in message news:i2nqdt$r7l$1(a)fred.mathworks.com... > "Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message > <i2ngh4$ki2$1(a)fred.mathworks.com>... >> "Jakob Jørgensen" <jhjspam(a)gmail.com> wrote in message >> <i2nes3$18v$1(a)fred.mathworks.com>... >> > I have a large, tall, sparse M-by-N matrix A, where for instance >> > M=100000 and N=80000. I want to determine whether or not rank(A) = N. >> > As far as I can see, the problem is equivalent to determining whether >> > or not the N'th singular value of A is nonzero. Commands rank and svd >> > fail since A is sparse, and svds tries to do a factorization, and runs >> > out of memory or takes forever. sprank gives an upper bound of the >> > rank, but returns N for my problem, which is not helpful. Any other >> > suggestions on how to determine whether or not rank(A) = N? Thank you. >> ================== >> >> I don't know how reliable this is for what you're trying to do, but if >> you do A\zeros(N,1) you should get a warning if the solver perceives A to >> be rank deficient. > > Thank you for the quick reply. > > For small A, your suggestion seems to work. Unfortunately Matlab returns a > segmentation fault, when I try to run it for my large A on our high memory > server. I will try to determine why. If possible, please send a MAT-file containing your A matrix, the segmentation fault log, the code you ran to produce that segmentation fault, and the output of the VER function to Technical Support. If the MAT-file is too large to send, you could instead describe how you construct A (and any structure it has that you know about.) > The documentation of mldivide says that the rank is determined from a QR > factorization with pivoting. Using R = qr(A,0) I get the same segmentation > fault. -- Steve Lord slord(a)mathworks.com comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ To contact Technical Support use the Contact Us link on http://www.mathworks.com |