cholesky decomposition,

[Soroosh]

Hi all,

I was trying to implement Fincke-Pohst in SAGE recently. In there, I need to calculate the cholesky decomposition of a real symmetric matrix of order n. In my case, this matrix is a real matrix. When I run

sage: A.cholesky()

it tells me that this is inexact field, and hence it fails. When I try

sage: A.cholesky_decomposition()

it warns that this is deprecated, although it does return an answer. My question is that, is there a non-deprecated function that can do this? (I'm using sage 5.9. I'm not sure what the status of the function is in 5.11)

Cheers,

Soroosh

PS. Here is my example.

<code>

M=diagonal_matrix([log(RR(2)),log(RR(3)),log(RR(5)),log(RR(7))])

M[-1]=[log(RR(2)),log(RR(3)),log(RR(5)),log(RR(7))]

A=M.transpose()*M

A.cholesky() <---- fails

A.cholesky_decomposition()   <---- gives warnings

</code>

[William Stein]

Without any thought, I recommend "use numpy", e.g.,

http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.cholesky.html

import numpy as np

A = np.array([[1,-2],[2,5]])

L = np.linalg.cholesky(A)

 

[Soroosh]

Cool. That will probably work.

Out of curiosity, can numpy handle multiprecision real numbers?

Soroosh

[William Stein]

On Thu, Aug 15, 2013 at 5:46 PM, Soroosh <syaz...@gmail.com> wrote:

Cool. That will probably work.

Out of curiosity, can numpy handle multiprecision real numbers?

No, unfortunately not.

 

--

You received this message because you are subscribed to the Google Groups "sage-nt" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-nt+u...@googlegroups.com.
To post to this group, send an email to sag...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-nt.
For more options, visit https://groups.google.com/groups/opt_out.

--
William Stein
Professor of Mathematics
University of Washington
http://wstein.org

[Martin Raum]

Hi Soroosh,

Early this year I implemented Fincke-Pohst in C++ and wrapped it in Cython. I have just pushed the result to my psage repository on github. https://github.com/martinra/psage/tree/shortvectorenum. You will then find the implementation in psage/lattice. The whole thing is intended for an ongoing project, so I have built in caching facilities, which you find in the folder precomputed_short_vectors. I attached the "specification" of the file format, which, I guess, is currently not included in the file.

The whole thing is quite extensively tested (using nosetest). When reading the code, be aware that the there is a C++ implementation (which I also use for a standalone program running on clusters). On top of this, there is a python wrapper, written in C++. On top of this, there is a Cython wrapper. I guess it would be best to move the Cython wrapper into the C++ file (templates and qualifiers in Cython are painful), but I haven't yet taken my time to do this.

The implementation is still slower than Magma's, but actually not too much. Since it used interval arithmetics, it will give provably correct results (up to potential bugs) - contrasting PARI, who have decided to compute at least all short vectors for given length, but possibly more.

Btw: As soon as we have a git development workflow, I will consider integrating this into Sage.

Martin

short vectors file:

File header:
  (unsighed int64) rank of lattice, which we denote by N
  (signed int64[N][N]) the entries of the lattice, where the necessary
                       symmetries must hold
  (unsigned int64) number of entries in the header
  (unsigned int64) NULL; For future use as pointer to continued header

  header entry:
    (unsigned int64) lenght of vectors
    (unsigned int64) number of vectors
    (unsigned int64) position of corresponding data in file

File body consists of vector data.
Vector data is a list of vectors, the number of which
  is determined by the header.
Vector:
  (signed int64[N]) Entries of vectors.
 
Vectors are stored up to multiplication by -1.

[Soroosh Yazdani]

Thanks Martin! That's awesome. Your code is exactly what I was looking for.

I did make a naive implementation of Fincke-Pohst in cython, but my implementation is running much, much slower than pari's implementation. It makes me feel I missed a step or two.

Cheers,

Soroosh

--