matrix diagonalization fails with simple matrix

[Zoltán]

Hi All,

I would like to diagonalize the following matrix

[1.0, 2, 2]
[ 2, 1.0, 0]
[ 2, 0, 1.0]

but sympy fails with it. Here is my traceback:

m = Matrix(3,3,[1., 2, 2, 2, 1., 0, 2, 0, 1.])
(p, d) = m.diagonalize()

NotImplementedError Traceback (most recent call
last)
/home/v923z/<ipython-input-32-d0af2185ab4d> in <module>()
----> 1 (p, d) = m.diagonalize()

/usr/local/lib/python2.7/dist-packages/sympy/matrices/matrices.pyc in
diagonalize(self, reals_only)
2790 if not self.is_square:
2791 raise NonSquareMatrixError()
-> 2792 if not self.is_diagonalizable(reals_only, False):
2793 self._diagonalize_clear_subproducts()
2794 raise MatrixError("Matrix is not diagonalizable")

/usr/local/lib/python2.7/dist-packages/sympy/matrices/matrices.pyc in
is_diagonalizable(self, reals_only, clear_subproducts)
2852 # self._diagonalize_clear_subproducts()

2853 # raise NotImplementedError("Symbolic matrices are
not implemented for diagonalization yet")

-> 2854 self._eigenvects = self.eigenvects(simplify=True)
2855 all_iscorrect = True
2856 for eigenval, multiplicity, vects in self._eigenvects:

/usr/local/lib/python2.7/dist-packages/sympy/matrices/matrices.pyc in
eigenvects(self, **flags)
2522 basis = tmp.nullspace(simplified=True)
2523 if not basis:
-> 2524 raise NotImplementedError("Can't evaluate
eigenvector for eigenvalue %s" % r)
2525 if simplify:
2526 # the relationship A*e = lambda*e will still
hold if we change the


NotImplementedError: Can't evaluate eigenvector for eigenvalue
-1.82842712474619

I have tried the same thing in octave, and the matrix that
diagonalizes the one above is
-7.0711e-01 5.5511e-17 7.0711e-01
5.0000e-01 -7.0711e-01 5.0000e-01
5.0000e-01 7.0711e-01 5.0000e-01

so, I suspect that the problem might come from an underflow in the
(1,2) entry, but I am not sure.

In any case, I wanted to ask if this is a bug, and whether there is a
way to overcome the problem. My sympy version is '0.7.1-git'.
Thanks,
Zoltán

[krastano...@gmail.com]

The error is present in the current master on git.

> --
> 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.
>

[Alexey U. Gudchenko]

It is related with this issue:

http://code.google.com/p/sympy/issues/detail?id=2193

20.01.2012 01:25, krastano...@gmail.com пишет:

--
Alexey U.

[Chris Smith]

>>> I would like to diagonalize the following matrix
>>>
>>> [1.0,   2,   2]
>>> [  2, 1.0,   0]
>>> [  2,   0, 1.0]
>>>
>>> but sympy fails with it. Here is my traceback:

If you change the 1.0 to 1 you get:

>>> d = [nsimplify(w) for w in [1., 2, 2, 2, 1., 0, 2, 0, 1.]]
>>> Matrix(3,3,d).diagonalize()
([ 0, sqrt(2), -sqrt(2)]
[-1, 1, 1]
[ 1, 1, 1], [1, 0, 0]
[0, 1 + 2*sqrt(2), 0]
[0, 0, -2*sqrt(2) + 1])
>>> out = _
>>> out[0]
[ 0, sqrt(2), -sqrt(2)]
[-1, 1, 1]
[ 1, 1, 1]
>>> out[1]
[1, 0, 0]
[0, 1 + 2*sqrt(2), 0]
[0, 0, -2*sqrt(2) + 1]

[Zoltán]

Right, but what if I do this

m = Matrix(3,3,[1., 2, 2, 2, 1., 0, 2, 0, 1.1]) ?

This matrix cannot be simplified in the way you suggested, and it
still fails.
Cheers,
Zoltán

[Zoltán]

Thanks for the reply! As it turns out, since I needed only the
eigenvalues, I could use the eigenvals() method, and that gave the
proper results. But it is still strange that the diagonalisation
routine should fail so miserably...
Cheers,
Zoltán

On Jan 20, 6:22 am, "Alexey U. Gudchenko" <pr...@goodok.ru> wrote:
> It is related with this issue:
>
> http://code.google.com/p/sympy/issues/detail?id=2193
>

> 20.01.2012 01:25, krastanov.ste...@gmail.com пишет:

[Chris Smith]

On Sat, Jan 21, 2012 at 3:31 AM, Zoltán <zvo...@gmail.com> wrote:
> Right, but what if I do this
>
> m = Matrix(3,3,[1., 2, 2, 2, 1., 0, 2, 0, 1.1]) ?
>

I'm not sure exactly what the difficulty is that you are imagining.
There are several routines in SymPy that don't work well when you use
floats, so you have to work around the "features'. For the above,

>>> Matrix(3,3,[nsimplify(w) for w in [1., 2, 2, 2, 1., 0, 2, 0, 1.1]]).diagonalize()

seems to work ok. Perhaps nsimlify could be made to recognize Matrix
so one could do nsimplify(m)...but just simplifying the data before
creating the Matrix isn't too bad.

[Zoltán]

> I'm not sure exactly what the difficulty is that you are imagining.
> There are several routines in SymPy that don't work well when you use
> floats, so you have to work around the "features'.

I don't think that "imagining" is the proper word here: the problem is
very real, as I have demonstrated in two simple cases. This is
certainly not a result of my imagination. I think this whole argument
could have been avoided, if it was stipulated in the documentation
that floats don't work. I have read http://docs.sympy.org/dev/modules/matrices.html
, but nowhere was any reference to problems related to float entries.

> >>> Matrix(3,3,[nsimplify(w) for w in [1., 2, 2, 2, 1., 0, 2, 0, 1.1]]).diagonalize()
>
> seems to work ok. Perhaps nsimlify could be made to recognize Matrix
> so one could do nsimplify(m)...but just simplifying the data before
> creating the Matrix isn't too bad.

Yes, this definitely works, but I could not imagine what could be
simpler than 1.1 as an entry. Besides, a call to nsimplify does not
make execution faster:(
In any case, the bottom line seems to be that sympy should not be used
for diagonalising float matrices.
Cheers,
Zoltán

[Alexey U. Gudchenko]

21.01.2012 16:06, Zoltán пишет:

I agreed that this problem exists: the problem with the `diagonalize`
method, which is based on the problem of the calculation of the
eigenvectors for some eigenvalues when the cases are not trivial. (And
more deeply, in common cases, in particular, it is (or was) related with
complex numbers, and with the so called zero-test problem: the complex
theoretical problem how to determine whether the symbolic expression is
zero or not). Now you added an example for the floats.

It is somehow described only in that issue [1]. Although only an
approximate location of the problem is described there.
I have accepted, and I am marked for this issue as owner, and fill
myself to continue with it.

And I agreed that the `nsimplify` can be considered only as a temporary
roundabout way.

As this issue is almost year old we must decide indeed, either to
describe it in the documentation or to increase the priority of the issue.

For my part, I am going to deal with it in a few weeks only (when I have
free time).

Now I attached your examples to that issue.

Thank you and others for raising the interest in this problem.

I hope it will be resolved soon.

[1] http://code.google.com/p/sympy/issues/detail?id=2193
--
Alexey U.

[Chris Smith]

On Sat, Jan 21, 2012 at 5:51 PM, Zoltán <zvo...@gmail.com> wrote:
>
>> I'm not sure exactly what the difficulty is that you are imagining.
>> There are several routines in SymPy that don't work well when you use
>> floats, so you have to work around the "features'.
>
> I don't think that "imagining" is the proper word here:

In our part of the country "imagining" means more along the lines of
"thinking". To you there is a problem. You posted it, I sent a
workaround, and then you sent another problem which was essentially
the same, so I didn't know what is what that you were thinking the
problem was.

To me it's just the common problem of having to use Rational instead
of floats. But this is not a special problem in the sense that python
addresses this in a variety of ways: int, long, float, decimal and I
showed a way to work around it in SymPy. What threw me off is the
"another problem" that you gave which is (to me) just the same problem
-- you have to be careful with how and where you use floats. There's
no easy, unambiguous way around that problem. Of course SymPy could
totally disallow the use of floats, but that's not pythonic. You can
use them at your own risk.

Anyway, we're always watching for ways to make SymPy more usable, and
do appreciate the feedback, so please don't imagine -- and now I do
mean "imagine" -- that what you have posted is unwelcome in anyway or
that I was trying to downplay it.

Best regards,
Chris

[Zoltán]

What threw me off is the
> "another problem" that you gave which is (to me) just the same problem
> -- you have to be careful with how and where you use floats.

I think the problem was that my original example used integers that
happened to be represented as floats, and to me it seemed that you
were trying to address that particular issue, so I attempted to
diagonalise another matrix, which had a real float, and that failed
again. What you pointed out about the representation of numbers in
python is reasonable.

> no easy, unambiguous way around that problem. Of course SymPy could
> totally disallow the use of floats, but that's not pythonic. You can
> use them at your own risk.

Hm. I think it would be a huge limitation on the use of sympy. I
really like the idea that one can change between symbolic and
numerical calculations.

> Anyway, we're always watching for ways to make SymPy more usable, and
> do appreciate the feedback,

I guess my only grievance was that this particular issue was not
mentioned in the documentation. As Alexey pointed out, it is OK, if
this issue is not fixed but then there should be a warning, so that
people are not left wondering about why something trivial does not
work.

On a more constructive note, would it be a viable option to add a non-
mandatory argument to diagonalize, and whenever that argument is set,
the diagonalisation routine would simplify all entries before
attempting to calculate anything. Something like this:

matrix.diagonalize(True)

If I am not mistaken, this would more or less be identical to the
solution you proposed above, with the exception that the
simplification would be hidden from the user.
Cheers,
Zoltán

[Joachim Durchholz]

Am 21.01.2012 15:18, schrieb Zolt�n:
> I really like the idea that one can change between symbolic and
> numerical calculations.

I don't think this can even be done. Not in the general case anyway;
simple cases like 1+1.0 will probably work as expected :-)

It's numeric instabilities that force the numeric solution algorithms
into fundamentally different approaches than their symbolic counterparts.
Matrix calculations are a case in point; computing determinants tends to
give rise to instabilities, and there's an entire lore around avoiding
having to compute determinants, and if it's unavoidable, how to
condition your equations so that the determinants are stable (i.e. the
denominators don't get near to zero).

I'm not sure that SymPy can do anything better than convert Python
floats to Rationals as soon as it tries to do arithmetic with them.

Regards,
Jo