headers
Nathan Bell | 1 May 05:47
Picon
Gravatar

Re: Sparse factorisation

On Thu, Apr 30, 2009 at 11:54 AM, Gael Varoquaux
<gael.varoquaux <at> normalesup.org> wrote:
>
> I have dug a bit further. It seems that all that needs to be done is
> expose the L, U, perm_c and perm_r attributes of the SciPyLUObject (AKA
> factored_lu object) returned by splu. Of course this easier said than
> done, as exposing these objects requires creating scipy sparse matrices
> from the inner SuperLU representation of these objects.
>
> I'd really appreciate if someone (most probably Nathan) could give me a
> hand with this. I realise that I am asking people to do my work, and I
> know exactly what my reaction is when someone comes around to me with
> this request (ie, not happy), but I am not sure I have to time required
> to learn the library and the tools to do this myself, and would hate to
> have to find an ugly workaround (this does sound like a bad excuse).
>

Hi Gael,

I hate to disappoint, but I don't have enough free time right now.  If
you or someone else wants to give it a try I can probably provide some
support.

I don't know what's "broken" in the umfpack wrappers, but I suspect
they may be the path of least resistance.

--

-- 
Nathan Bell wnbell <at> gmail.com
http://graphics.cs.uiuc.edu/~wnbell/
(Continue reading)

Permalink | Reply | 3 comments |
headers
Ondrej Certik | 1 May 05:52
Picon
Gravatar

Re: Sparse factorisation

On Thu, Apr 30, 2009 at 6:56 AM, Gael Varoquaux
<gael.varoquaux <at> normalesup.org> wrote:
> Hi there,
>
> This is a rather general question about sparse factorisation and Scipy. I
> no nothing about sparse linear algebra, but I happen to need to do a
> sparse Cholesky.
>
> By sparse Cholesky, I don't mean a sparse implementation of Cholesky, my
> matrices are small, I don't actually care about sparse linear algebra,
> but I mean finding P, permutation matrix, and L, lower-only matrix, such
> that the sparse positive definite matrix A can be written:
>        A = P L L' P',
> where P is a reordering of the rows and column of A used to maximise the
> sparsity of L.
>
> What is really of interest to me is P. P can be computed using the colamd
> routine, present in UMFPACK. Scipy sparse does make use of the colamd
> routine (by setting the permc_spec argument of scipy.sparse.splu to 3).
> However, I cannot figure out a way of extracting P.
>
> I have been looking around in the sparse code source, as well as the
> scikits.umfpack code source, and I must admit I am a bit at loss to what
> is the best way to achieve my goals. The umfpack.py module, both in the
> scipy code base and in the scikit, is broken. Scipy seems to hint that it

Is umfpack in scipy really broken? That is very disappointing. I was
planning to use umfpack through scipy. I hope it will not be difficult
for me to fix it.
(Continue reading)

Permalink | Reply | 22 comments |
headers
David Warde-Farley | 1 May 06:10
Picon
Favicon
Gravatar

Re: Sparse factorisation


On 30-Apr-09, at 11:52 PM, Ondrej Certik wrote:

> Is umfpack in scipy really broken? That is very disappointing. I was
> planning to use umfpack through scipy. I hope it will not be difficult
> for me to fix it.

I don't use it, but a friend of mine was having trouble building SciPy  
from svn recently and UMFPACK turned out to be the culprit. I have  
been planning to open a ticket, but I haven't had a chance to try it  
out myself and confirm the behaviour.

David
Permalink | Reply | 0 comments |
headers
Gael Varoquaux | 1 May 10:17
Favicon
Gravatar

Re: Sparse factorisation

On Thu, Apr 30, 2009 at 11:47:38PM -0400, Nathan Bell wrote:
> On Thu, Apr 30, 2009 at 11:54 AM, Gael Varoquaux
> <gael.varoquaux <at> normalesup.org> wrote:

> > I have dug a bit further. It seems that all that needs to be done is
> > expose the L, U, perm_c and perm_r attributes of the SciPyLUObject (AKA
> > factored_lu object) returned by splu. Of course this easier said than
> > done, as exposing these objects requires creating scipy sparse matrices
> > from the inner SuperLU representation of these objects.

> > I'd really appreciate if someone (most probably Nathan) could give me a
> > hand with this. I realise that I am asking people to do my work, and I
> > know exactly what my reaction is when someone comes around to me with
> > this request (ie, not happy), but I am not sure I have to time required
> > to learn the library and the tools to do this myself, and would hate to
> > have to find an ugly workaround (this does sound like a bad excuse).

> I hate to disappoint, but I don't have enough free time right now.  If
> you or someone else wants to give it a try I can probably provide some
> support.

Fair enough. I can understand that.

In the mean time, I figured out that I only needed to expose thee
permutation vectors from the factored_lu object. These permutation vector
are not sparse objects, so it is very easy to expose them to numpy (just
a question of PyArray_SimpleNewFromData and using the 'base' trick (
http://blog.enthought.com/?p=62 ) to do the reference counting. 

This is a bit ugly, because it does not expose all the information.
(Continue reading)

Permalink | Reply | 2 comments |
headers
Gael Varoquaux | 1 May 10:22
Favicon
Gravatar

Re: Sparse factorisation

On Thu, Apr 30, 2009 at 08:52:28PM -0700, Ondrej Certik wrote:
> Is umfpack in scipy really broken? That is very disappointing. 

There are some signs of it being depreciated in the source code.

> I was planning to use umfpack through scipy. I hope it will not be
> difficult for me to fix it.

I had an email exchange with one of the pysparse authors, and it does
seem that the pysparse is currently well-maintained, and the latest
release may be exposing the best wrappers to umfpack. In the interest of
reducing duplication and exposing a coherent view to our users, it might
be good to focus on pysparse for sparse matrix functionnality not exposed
in scipy.

This brings to an important point: how do we list semi-official extension
to scipy that we know are well QAed and complete well the functionnality
of scipy?

Gaël
Permalink | Reply | 20 comments |
headers
David Cournapeau | 1 May 10:08
Picon
Picon

Re: Sparse factorisation

Gael Varoquaux wrote:
> On Thu, Apr 30, 2009 at 08:52:28PM -0700, Ondrej Certik wrote:
>   
>> Is umfpack in scipy really broken? That is very disappointing. 
>>     
>
> There are some signs of it being depreciated in the source code.
>   

umfpack is under the GPL. My understanding is that new contributions to
umfpack should go to the corresponding scikit.

cheers,

David
Permalink | Reply | 19 comments |
headers
Gael Varoquaux | 1 May 10:33
Favicon
Gravatar

Re: Sparse factorisation

On Fri, May 01, 2009 at 05:08:37PM +0900, David Cournapeau wrote:
> Gael Varoquaux wrote:
> > On Thu, Apr 30, 2009 at 08:52:28PM -0700, Ondrej Certik wrote:

> >> Is umfpack in scipy really broken? That is very disappointing. 

> > There are some signs of it being depreciated in the source code.

> umfpack is under the GPL. My understanding is that new contributions to
> umfpack should go to the corresponding scikit.

Right, new versions of UMFPACK are indeed under GPL. I just noticed that.
Have you tried contacting Tim Davis about that. The web page does suggest
to contacting him for distributing UMFPACK under a different license. My
guess is that distirbuting UMFPACK under BSD with scipy would be contrary
to distributing it under GPL in the main download page, so our chances
are small.

Now I understand better the switch to SuperLU.

Gaël
Permalink | Reply | 18 comments |
headers
Robert Cimrman | 1 May 11:07
Picon

Re: Sparse factorisation

Gael Varoquaux wrote:
> On Fri, May 01, 2009 at 05:08:37PM +0900, David Cournapeau wrote:
>> Gael Varoquaux wrote:
>>> On Thu, Apr 30, 2009 at 08:52:28PM -0700, Ondrej Certik wrote:
> 
>>>> Is umfpack in scipy really broken? That is very disappointing. 
> 
>>> There are some signs of it being depreciated in the source code.
> 
>> umfpack is under the GPL. My understanding is that new contributions to
>> umfpack should go to the corresponding scikit.
> 
> Right, new versions of UMFPACK are indeed under GPL. I just noticed that.
> Have you tried contacting Tim Davis about that. The web page does suggest
> to contacting him for distributing UMFPACK under a different license. My
> guess is that distirbuting UMFPACK under BSD with scipy would be contrary
> to distributing it under GPL in the main download page, so our chances
> are small.
> 
> Now I understand better the switch to SuperLU.

Yes, we have tried to contact Tim Davis (me and Nathan Bell, if I 
remeber correctly?), but he refused to provide an exception for scipy. 
That is why the umfpack wrappers (written quite a time ago by me) were 
changed into the umfpack scikit. Now there is a situation when some 
version of the wrappers still survives in scipy proper, while another 
version is in the scikit. This is not good, I know, but I have been very 
busy recently. But I use the wrappers on a daily basis without problem 
myself - if you have problems, open a ticket, please, and notify me.
(Continue reading)

Permalink | Reply | 17 comments |
headers
Lubos Vrbka | 1 May 10:56

fourier transform in numpy/scipy

hi guys,

are the discrete fourier transform pairs in numpy/scipy self adjoint
(np.fft.fft - np.fft.ifft)? i read somewhere, that this property doesn't
necessarily have to be fulfilled. i performed some tests and it seems
that it indeed is so, but i wanted to be sure...

best,

--

-- 
Lubos _ <at> _"
http://www.lubos.vrbka.net
Permalink | Reply | 0 comments |
headers
Gael Varoquaux | 1 May 11:43
Favicon
Gravatar

Re: Sparse factorisation

On Fri, May 01, 2009 at 11:07:26AM +0200, Robert Cimrman wrote:
> if you have problems, open a ticket, please, and notify me.

http://projects.scipy.org/scipy/ticket/935

Stefan is looking at fixing a few things with the scikit. IMHO the
UMFPACK code in scipy should either be killed or fixed.

By the way, you need to ad an 'import nose' at the end of
'test_umfpack.py' in the scikit, because nose is no longer imported in
numpy.testing.

Thanks for all your work,

Gaël
Permalink | Reply | 16 comments |

Gmane