Improving the function for calculating eigenvectors and QR factorization

[Sahil Singla]

I am a second year undergraduate ,doing my major in Computer Sciences.
I have done a course in Numerical Linear Algebra and Scientific Computing and it is by far the most interesting topic I have encountered in my life.

I want to improve the code of QR factorization used in sympy by using Householder reflections to find the QR of a matrix. Ideally the naiive algorithm in sympy has complexity of O(m*n*n). Using householder transformation method, this can be made much faster. 

Apart from it, I intend to implement the method for finding the eigenvalues of a real matrix using a different algorithm. The present algorithm which computes the eigenvalues of a matrix by finding the roots of the characteristic polynomial is highly unstable and has a very high condition number. The condition number for calculating an eigenvalue of a matrix(say e) when the i'th coefficient in the matrix (say a(i)) is perturbed is (a(i)*(e^(i-1)))/(p'(e)) where p'(e) is the derivative of the characteristic polynomial at e. Seeing the amount of rounding errors that will be involved in calculation of a particular coefficient, the resulting error in the computation of the eigenvalue will be huge!

The method by which it can be simplified is the power iteration method. The program will have to be a combination of many different methods so that eigenvalues for different matrices . In particular, the method will sort out the issue 3189.

The method I intend to use is  QR Algorithm for real matrices. After applying finite preliminary reductions, the matrix can be converted into its hessenberg form so that the work per iteration (step: A=R*Q) is reduced to O(n^2) from O(n^3). 

Regards,

Sahil Singla

Sophomore,

Computer Sciences,

IIT Delhi

[Sahil Singla]

[Tom Bachmann]

Please note that SymPy is mainly a *symbolic* algebra package, so
iterative and numerical algorithms are not really that useful to us...

> --
> You received this message because you are subscribed to the Google
> Groups "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to
> sympy+un...@googlegroups.com.
> For more options, visit this group at
> http://groups.google.com/group/sympy?hl=en.

[Matthew Rocklin]

Hi Sahil, 

The algorithms you mention are indeed excellent for computing on numerical matrices. They are mindful of numerical rounding errors. 

SymPy focuses on exact symbolic algebra. We don't have issues with rounding errors but instead have issues with building up very large expressions. For this reason algorithms like the power method aren't appropriate for us. The algorithms you mention belong in projects like SciPy. 

I also come from a numerical linear algebra and scientific computing background. It is often surprising to see differences between the symbolic and numeric worlds. 

Of course, there is still a lot of work to do in SymPy and Linear Algebra, and I in particular would be happy to find GSoC students who were interested in this field. 

Best,

-Matt

[Alexey U. Gudchenko]

On 27.03.2012 20:22, Sahil Singla wrote:
> I am a second year undergraduate ,doing my major in Computer Sciences.
> I have done a course in Numerical Linear Algebra and Scientific Computing
> and it is by far the most interesting topic I have encountered in my life.
>
> I want to improve the code of QR factorization used in sympy by using
> Householder reflections to find the QR of a matrix. Ideally the naiive
> algorithm in sympy has complexity of O(m*n*n). Using householder
> transformation method, this can be made much faster.
>
> Apart from it, I intend to implement the method for finding the eigenvalues
> of a real matrix using a different algorithm. The present algorithm which
> computes the eigenvalues of a matrix by finding the roots of the
> characteristic polynomial is highly unstable and has a very high condition
> number. The condition number for calculating an eigenvalue of a matrix(say
> e) when the i'th coefficient in the matrix (say a(i)) is perturbed is
> (a(i)*(e^(i-1)))/(p'(e)) where p'(e) is the derivative of the
> characteristic polynomial at e. Seeing the amount of rounding errors that
> will be involved in calculation of a particular coefficient, the resulting
> error in the computation of the eigenvalue will be huge!
>
> The method by which it can be simplified is the power iteration method. The
> program will have to be a combination of many different methods so that
> eigenvalues for different matrices . In particular, the method will sort
> out the issue 3189.
>

It will be good!
But I have to note, that the main task of the issue 3189 to give a
chance to obtain the eigenvectors (not values) in several cases, when
the nullspace algorithm failed.

-
Alexey Gudchenko

[someone]

Hi,

> Please note that SymPy is mainly a *symbolic* algebra package, so
> iterative and numerical algorithms are not really that useful to us...

Except maybe if we want to provide arbitrary-precision versions.
Such versions are obviously not included in scipy.

But then there is the question if it is possible and fruitful
to just replace floats by mpfloats/bigfloatsor if this gives
bad behaviour.

[Joachim Durchholz]

Am 27.03.2012 21:14, schrieb someone:
> But then there is the question if it is possible and fruitful
> to just replace floats by mpfloats/bigfloatsor if this gives
> bad behaviour.

One bad behaviour would be getting fractions with extremely large
numerators and denominators.
Floats are essentially trading precision for a smaller memory footprint.

[Chris Smith]

The limit_denominator can help with this:

>>> Rational(2**40,3**39)

1099511627776/4052555153018976267

>>> _.limit_denominator(10**15)

244736711/902045048220285

>>> _.limit_denominator(10**12)

77999/287486954567