From: <do...@vi...> - 2002-03-26 14:33:59
|
Hi, I have some more questions to vnl_sparse_matrix_linear_system. I have an impression that this class is only an adaptor to vnl_linear_system. So the question is: - does vnl_sparse_matrix_linear_system adress the sparseness? I need to solve a sparse matrix like 50,000 x 50,000. Will it handle it? - What is the safe way to get x from Ax=b and/or how do I get inverse of sparse A in order to use mult from vnl_sparse_matrix to get x=A^{-1}*b. thank you Dominique |
From: Ian S. <ian...@st...> - 2002-03-26 18:26:28
|
> I have some more questions to vnl_sparse_matrix_linear_system. > I have an impression that this class is only an adaptor to > vnl_linear_system. So the question is: > > - does vnl_sparse_matrix_linear_system adress the sparseness? > I need to > solve a sparse matrix like 50,000 x 50,000. Will it handle it? You'll need to try - it depends on the amount of memory, your matrix density, and the particular implementation details. > > - What is the safe way to get x from Ax=b and/or how do I get > inverse of > sparse A in order to use mult from vnl_sparse_matrix to get > x=A^{-1}*b. I'd suggest using the VXL documentation search engine to answer questions like this. http://www.isbe.man.ac.uk/search_vxl.html http://www.isbe.man.ac.uk/cgi-bin/htsearch?config=htdig-vxl;words=sparse%20v nl; seems to suggest a few possibilities, e.g. vnl_lsqr Ian. |
From: <do...@vi...> - 2002-03-26 18:57:33
|
> I'd suggest using the VXL documentation search engine to answer questions > like this. I'm sorry, but there is no clear answer there to "questions like this". I see a lot of classes, but none is giving me either an inverse of a big *sparse* matrix or a corresponding solution - not to mention other expectations/doubts I wrote about :( Dominique |
From: Ian S. <ian...@st...> - 2002-03-27 09:35:38
|
> I see a lot of classes, but none is giving me either an > inverse of a big > *sparse* matrix or a corresponding solution - not to mention other > expectations/doubts I wrote about :( I think vnl_lsqr is the class you are looking for. From the documentation http://www.isbe.man.ac.uk/public_vxl_doc/vxl/vnl/html/class_vnl_lsqr.html "vnl_lsqr implements an algorithm for large, sparse linear systems and sparse, linear least squares" As I undestand it, linear least squares applied to the equation Ax=b will give you a least squares solution for x, which is as good as you are likely to get for a massive sparse matrix A. > I'm sorry, but there is no clear answer there to "questions > like this". The place for an overview of numerical algorithms is a book like "Numerical Recipes." There is no point in us replicating this sort of information in VXL's documentation. However, the search engine and documentation do provide a means of finding algorithms. vnl_lsqr comes 4th on a search for "vnl sparse". Ok, the vnl_lsqr documentation doesn't mention the word "inverse" (because it doesn't actually do an inverse,) but this is exactly the process I used to answer your question. I readily acknowledge that VXL's documentation is far from perfect. If you can think of some better documentation (either for the vnl chapter in the book, or the doxygen output,) we would be grateful. Improvements to the VXL code or documentation are always welcome. Ian. |
From: <do...@vi...> - 2002-03-27 12:50:43
|
> I think vnl_lsqr is the class you are looking for. From the documentation Ian, yes, you are right. I was there before, but I didnt know I can use it (well, matix inverse and least squares are kind of two different jobs, arent they?). I checked - it works for small systems, but if it works good enough for huge sparse systems I will know in a day or two. > I readily acknowledge that VXL's documentation is far from perfect. If you > can think of some better documentation (either for the vnl chapter in the > book, or the doxygen output,) we would be grateful. Improvements to the VXL > code or documentation are always welcome. I dont need a perfect documentation (although a little bit less messy would be great). Instead I prefer the idea of getting in touch with someone who had similar problem (and asking him for solution). I can offer the same. thank you Dominique |
From: <do...@vi...> - 2002-03-27 13:26:31
|
vnl_sqrt doesnt solve Ax=b for systems aleardy around 30. For 20 it gives good results (for a naive linear problem) and: vnl_lsqr.cxx : x = 0 is the exact solution. No iterations were performed. vnl_lsqr.cxx : The iteration limit ITNLIM was reached. vnl_lsqr.cxx : residual norm estimate = 8.40484e-06 vnl_lsqr.cxx : result norm estimate = 14.6304 vnl_lsqr.cxx : condition no. estimate = 67.683 vnl_lsqr.cxx : iterations = 20 whereas for 30 equations it goes banana (wrong solution): vnl_lsqr.cxx : x = 0 is the exact solution. No iterations were performed. vnl_lsqr.cxx : The iteration limit ITNLIM was reached. vnl_lsqr.cxx : residual norm estimate = 3.24731 vnl_lsqr.cxx : result norm estimate = 14.2702 vnl_lsqr.cxx : condition no. estimate = 106.047 vnl_lsqr.cxx : iterations = 30 I guess residual norm estimate holds the key. Any ideas? thanks Dominique |