smallest singular value of large, tall, sparse matrix
in [Matlab]
|
Prev: NI USB-6501
Next: Multiple Tabs with dlmwrite
From: Jakob Heide Jørgensen on 28 Jul 2010 13:31 Thank you all for you suggestions. I am currently running Bruno's #1 s = svds(A,1,0) #2 Afun=@(x) A\(x'/A)'; s = sqrt(eigs(Afun,size(A,2),1,'sm')) and Pat's #3 op = @(x) A'*(A*x); 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); and will report my results here when completed. Steve: I have submitted a bug report on the segmentation fault.
From: Jakob Heide Jørgensen on 1 Aug 2010 04:27 Here are my results for computing the smallest singular value: > Bruno's #1 > s = svds(A,1,0) failed after ~ 8 h 40 mins with the error ??? 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; Error in ==> svds at 140 [W,D,bflag] = eigs(B,bk,bsigma,boptions); > #2 > Afun=@(x) A\(x'/A)'; > s = sqrt(eigs(Afun,size(A,2),1,'sm')) failed after ~ 10 mins with the error: ??? Error using ==> mrdivide Out of memory. Type HELP MEMORY for your options. Error in ==> @(x) A\(x'/A)' Error in ==> eigs>AminusSigmaBsolve at 1192 v = A(u,varargin{afunNargs}); Error in ==> eigs at 257 workd(:,cols(2)) = AminusSigmaBsolve(workd(:,cols(1))); > and Pat's #3 > op = @(x) A'*(A*x); > 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); is still running on the 87th hour, slowly making progress. The last iterations are Iteration 3854: a few Ritz values of the 100-by-100 matrix: 2.3122e-05 Iteration 3855: a few Ritz values of the 100-by-100 matrix: 2.3122e-05 Iteration 3856: a few Ritz values of the 100-by-100 matrix: 2.3121e-05 Iteration 3857: a few Ritz values of the 100-by-100 matrix: 2.3120e-05 I think it is still too soon to conclude that the value is nonzero, but if the iterates level here I am quite happy, since the current value agrees well with smaller scale results.
From: Pat Quillen on 2 Aug 2010 10:06 "Jakob Heide Jørgensen" <jhjspam(a)gmail.com> wrote in message <i33b4o$8ha$1(a)fred.mathworks.com>... <snip for brevity...> > Iteration 3857: a few Ritz values of the 100-by-100 matrix: > 2.3120e-05 > > I think it is still too soon to conclude that the value is nonzero, but if the iterates level here I am quite happy, since the current value agrees well with smaller scale results. Hi Jakob. Bear in mind that the value here (as in Bruno's suggestion) is the square of the smallest singular value, thus your smallest singular value is probably more like 4e-3. Glad to see that this is working (albeit slowly!) for you. Take care. Pat.
From: Jakob Heide Jørgensen on 2 Aug 2010 10:35 > Bear in mind that the value here (as in Bruno's suggestion) is the square of the smallest singular value, thus your smallest singular value is probably more like 4e-3. Glad to see that this is working (albeit slowly!) for you. Hi Pat Thanks for the heads up! I'm aware :) Unfortunately, the computation has now terminated after the maximum 5000 iterations (4 days, 12 hours) without having converged. As you said it would, eigs returned 0, and the last printed output was Iteration 5001: a few Ritz values of the 100-by-100 matrix: 2.2510e-05 I am now trying a run with the default 20-by-20 matrix instead and up to 50000 iterations. As could be expected, each iteration is much faster, but convergence is slower. I am currently at (~ 8 hours) Iteration 3810: a few Ritz values of the 20-by-20 matrix: 1.012732777655090e-04 where for 100-by-100 I had Iteration 3810: a few Ritz values of the 100-by-100 matrix: 2.3156e-05 Perhaps a p somewhere in between 20 and 100 would be a good trade-off. What made you pick p=100?
From: Pat Quillen on 2 Aug 2010 16:09
"Jakob Heide Jørgensen" <jhjspam(a)gmail.com> wrote in message <i36l3f$l4l$1(a)fred.mathworks.com>... > > > Bear in mind that the value here (as in Bruno's suggestion) is the square of the smallest singular value, thus your smallest singular value is probably more like 4e-3. Glad to see that this is working (albeit slowly!) for you. > > Hi Pat > > Thanks for the heads up! I'm aware :) > > Unfortunately, the computation has now terminated after the maximum 5000 iterations (4 days, 12 hours) without having converged. As you said it would, eigs returned 0, and the last printed output was > > Iteration 5001: a few Ritz values of the 100-by-100 matrix: > 2.2510e-05 > > I am now trying a run with the default 20-by-20 matrix instead and up to 50000 iterations. As could be expected, each iteration is much faster, but convergence is slower. I am currently at (~ 8 hours) > > Iteration 3810: a few Ritz values of the 20-by-20 matrix: > 1.012732777655090e-04 > > where for 100-by-100 I had > > Iteration 3810: a few Ritz values of the 100-by-100 matrix: > 2.3156e-05 > > Perhaps a p somewhere in between 20 and 100 would be a good trade-off. What made you pick p=100? Hi Jakob. The larger the p the better. Lanczos tends to resolve eigenvalues on the ends of the spectrum faster than in the middle, so the larger the subspace you choose the better off you are. Of course, bear in mind that too large and you may not be able to stomach the workspace. I have no reason for choosing 100 in particular; I just blindly picked something larger than 20. Basically, Lanczos will pick out p eigenvalues to approximate the 80000 eigenvalues that A'*A has. Perhaps crank it up to 200 and crank maxit up as high as you are willing to wait. Something else you can try, although I'm not sure you'll be able to complete it is: p = colamd(A); % Try to get R as in A(:,p) = QR. p is an ordering which hopefully will encourage % sparser factors and let you finish. If so, then % A(:,p)'*A(:,p) = R'*R, which will have the same eigenvalues as A'*A. R = qr(A(:,p),0); L = R'; op = @(x) R\(L\x); opts.isreal = true; opts.issym = true; opts.disp = 1; opts.p = 100; opts.maxit = 50000; eigs(op, size(A,2), 1, 'sm', opts); Here, you're doing the same thing, but now you are doing shift and invert Lanczos. If you can compute the R from QR then convergence will be much quicker than the 'sa' version (although each step will cost more). I'm afraid, however, that you may not have enough memory to compute R. Good luck and let us know how you fare. Pat. |