New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Wrong Jordan form: complex eigenvalues w/ geo. mult. > alg. mult. #9274

Closed
cbm755 opened this Issue Apr 9, 2015 · 5 comments

Comments

Projects
None yet
2 participants
@cbm755
Copy link
Contributor

cbm755 commented Apr 9, 2015

Consider the following matrix:

In [172]: L = Matrix([[2, 4], [-4, 2]])

In [173]: Z = zeros(2,2);  I = Identity(2)

In [174]: A = Matrix(BlockMatrix([[L, I], [Z, L]]))

In [175]: pprint(A)
⎡2   4  1   0⎤
⎢            ⎥
⎢-4  2  0   1⎥
⎢            ⎥
⎢0   0  2   4⎥
⎢            ⎥
⎣0   0  -4  2⎦

Compute its Jordan form:

In [177]: pprint(A.jordan_form())
⎛⎡0  ⅈ  ⅈ  0⎤, ⎡2 - 4⋅ⅈ     1        0        0        0        0        0        0        0        0   ⎤⎞
⎜⎢          ⎥  ⎢                                                                                        ⎥⎟
⎜⎢0  1  1  0⎥  ⎢   0     2 - 4⋅ⅈ     0        0        0        0        0        0        0        0   ⎥⎟
⎜⎢          ⎥  ⎢                                                                                        ⎥⎟
⎜⎢0  0  0  ⅈ⎥  ⎢   0        0     2 - 4⋅ⅈ     1        0        0        0        0        0        0   ⎥⎟
⎜⎢          ⎥  ⎢                                                                                        ⎥⎟
⎜⎣0  0  0  1⎦  ⎢   0        0        0     2 - 4⋅ⅈ     0        0        0        0        0        0   ⎥⎟
⎜              ⎢                                                                                        ⎥⎟
⎜              ⎢   0        0        0        0     2 - 4⋅ⅈ     0        0        0        0        0   ⎥⎟
⎜              ⎢                                                                                        ⎥⎟
⎜              ⎢   0        0        0        0        0     2 + 4⋅ⅈ     1        0        0        0   ⎥⎟
⎜              ⎢                                                                                        ⎥⎟
⎜              ⎢   0        0        0        0        0        0     2 + 4⋅ⅈ     0        0        0   ⎥⎟
⎜              ⎢                                                                                        ⎥⎟
⎜              ⎢   0        0        0        0        0        0        0     2 + 4⋅ⅈ     1        0   ⎥⎟
⎜              ⎢                                                                                        ⎥⎟
⎜              ⎢   0        0        0        0        0        0        0        0     2 + 4⋅ⅈ     0   ⎥⎟
⎜              ⎢                                                                                        ⎥⎟
⎝              ⎣   0        0        0        0        0        0        0        0        0     2 + 4⋅ⅈ⎦⎠

This is wrong. The J matrix is too big. And there is a zero (generalized) eigenvector.

So far, I've traced it back to here:

In [178]: pprint(A._jordan_block_structure())
⎧2 - 4⋅ⅈ: ⎧1: ⎡⎡⎡ⅈ⎤⎤⎤, 2: ⎡⎡⎡0⎤, ⎡ⅈ⎤⎤, ⎡⎡ⅈ⎤, ⎡0⎤⎤⎤⎫, 2 + 4⋅ⅈ: ⎧1: ⎡⎡⎡-ⅈ⎤⎤⎤, 2: ⎡⎡⎡0⎤, ⎡-ⅈ⎤⎤, ⎡⎡-ⅈ⎤, ⎡0 ⎤⎤⎤⎫⎫
⎪         ⎪   ⎢⎢⎢ ⎥⎥⎥     ⎢⎢⎢ ⎥  ⎢ ⎥⎥  ⎢⎢ ⎥  ⎢ ⎥⎥⎥⎪           ⎪   ⎢⎢⎢  ⎥⎥⎥     ⎢⎢⎢ ⎥  ⎢  ⎥⎥  ⎢⎢  ⎥  ⎢  ⎥⎥⎥⎪⎪
⎪         ⎪   ⎢⎢⎢1⎥⎥⎥     ⎢⎢⎢0⎥  ⎢1⎥⎥  ⎢⎢1⎥  ⎢0⎥⎥⎥⎪           ⎪   ⎢⎢⎢1 ⎥⎥⎥     ⎢⎢⎢0⎥  ⎢1 ⎥⎥  ⎢⎢1 ⎥  ⎢0 ⎥⎥⎥⎪⎪
⎨         ⎨   ⎢⎢⎢ ⎥⎥⎥     ⎢⎢⎢ ⎥  ⎢ ⎥⎥  ⎢⎢ ⎥  ⎢ ⎥⎥⎥⎬           ⎨   ⎢⎢⎢  ⎥⎥⎥     ⎢⎢⎢ ⎥  ⎢  ⎥⎥  ⎢⎢  ⎥  ⎢  ⎥⎥⎥⎬⎬
⎪         ⎪   ⎢⎢⎢0⎥⎥⎥     ⎢⎢⎢0⎥  ⎢0⎥⎥  ⎢⎢0⎥  ⎢ⅈ⎥⎥⎥⎪           ⎪   ⎢⎢⎢0 ⎥⎥⎥     ⎢⎢⎢0⎥  ⎢0 ⎥⎥  ⎢⎢0 ⎥  ⎢-ⅈ⎥⎥⎥⎪⎪
⎪         ⎪   ⎢⎢⎢ ⎥⎥⎥     ⎢⎢⎢ ⎥  ⎢ ⎥⎥  ⎢⎢ ⎥  ⎢ ⎥⎥⎥⎪           ⎪   ⎢⎢⎢  ⎥⎥⎥     ⎢⎢⎢ ⎥  ⎢  ⎥⎥  ⎢⎢  ⎥  ⎢  ⎥⎥⎥⎪⎪
⎩         ⎩   ⎣⎣⎣0⎦⎦⎦     ⎣⎣⎣0⎦  ⎣0⎦⎦  ⎣⎣0⎦  ⎣1⎦⎦⎦⎭           ⎩   ⎣⎣⎣0 ⎦⎦⎦     ⎣⎣⎣0⎦  ⎣0 ⎦⎦  ⎣⎣0 ⎦  ⎣1 ⎦⎦⎦⎭⎭
@jksuom

This comment has been minimized.

Copy link
Member

jksuom commented Apr 9, 2015

It seems that the bug is in the implementation of what these lines explain:

# We want the vectors in `Kernel((self-lI)^s)` (**),
# but without those in `Kernel(self-lI)^s-1` so we will add  these as additional equations
# to the sytem formed by `S`

The idea of adding new equations is to restrict the null space of (self-lI)^s to those vectors that are orthogonal to Kernel(self-lI)^s-1. This works well in the real case, but vectors in the complex domain can be 'orthogonal' to themselves unless hermitian metric is applied. So, instead of
adding the vectors of Kernel(self-lI)^s-1 as such, their complex conjugates should be used here:

for k in range(0, a[s-1]):
    S = S.col_join((exclude_vectors[k]).transpose())
@cbm755

This comment has been minimized.

Copy link
Contributor

cbm755 commented Apr 18, 2015

I changed to S = S.col_join((exclude_vectors[k]).H) and that does look better, although not quite right yet:

In [13]: pprint(A.jordan_cells())
⎛⎡ⅈ  0  ⅈ  -ⅈ⎤, ⎡⎡2 - 4⋅ⅈ     1   ⎤, [2 - 4⋅ⅈ], ⎡2 + 4⋅ⅈ     1   ⎤, [2 + 4⋅ⅈ]⎤⎞
⎜⎢           ⎥  ⎢⎢                ⎥             ⎢                ⎥           ⎥⎟
⎜⎢1  0  1  1 ⎥  ⎣⎣   0     2 - 4⋅ⅈ⎦             ⎣   0     2 + 4⋅ⅈ⎦           ⎦⎟
⎜⎢           ⎥                                                                ⎟
⎜⎢0  ⅈ  0  0 ⎥                                                                ⎟
⎜⎢           ⎥                                                                ⎟
⎝⎣0  1  0  0 ⎦                                                                ⎠

In [14]: pprint(A.jordan_form())
⎛⎡ⅈ  0  ⅈ  -ⅈ⎤, ⎡2 - 4⋅ⅈ     1        0        0        0        0   ⎤⎞
⎜⎢           ⎥  ⎢                                                    ⎥⎟
⎜⎢1  0  1  1 ⎥  ⎢   0     2 - 4⋅ⅈ     0        0        0        0   ⎥⎟
⎜⎢           ⎥  ⎢                                                    ⎥⎟
⎜⎢0  ⅈ  0  0 ⎥  ⎢   0        0     2 - 4⋅ⅈ     0        0        0   ⎥⎟
⎜⎢           ⎥  ⎢                                                    ⎥⎟
⎜⎣0  1  0  0 ⎦  ⎢   0        0        0     2 + 4⋅ⅈ     1        0   ⎥⎟
⎜               ⎢                                                    ⎥⎟
⎜               ⎢   0        0        0        0     2 + 4⋅ⅈ     0   ⎥⎟
⎜               ⎢                                                    ⎥⎟
⎝               ⎣   0        0        0        0        0     2 + 4⋅ⅈ⎦⎠

It is still giving erroneous multiplicity-1 blocks for each eigenvalue.

@jksuom

This comment has been minimized.

Copy link
Member

jksuom commented Apr 19, 2015

Those vectors are added in two places, here and here. After changing transpose() to adjoint() in both lines I seem to get quite reasonable looking results.

In [5]: pprint(A.jordan_cells())
⎛⎡ⅈ  0  -ⅈ  0 ⎤, ⎡⎡2 - 4⋅ⅈ     1   ⎤, ⎡2 + 4⋅ⅈ     1   ⎤⎤⎞
⎜⎢            ⎥  ⎢⎢                ⎥  ⎢                ⎥⎥⎟
⎜⎢1  0  1   0 ⎥  ⎣⎣   0     2 - 4⋅ⅈ⎦  ⎣   0     2 + 4⋅ⅈ⎦⎦⎟
⎜⎢            ⎥                                          ⎟
⎜⎢0  ⅈ  0   -ⅈ⎥                                          ⎟
⎜⎢            ⎥                                          ⎟
⎝⎣0  1  0   1 ⎦                                          ⎠

In [6]: pprint(A.jordan_form())
⎛⎡ⅈ  0  -ⅈ  0 ⎤, ⎡2 - 4⋅ⅈ     1        0        0   ⎤⎞
⎜⎢            ⎥  ⎢                                  ⎥⎟
⎜⎢1  0  1   0 ⎥  ⎢   0     2 - 4⋅ⅈ     0        0   ⎥⎟
⎜⎢            ⎥  ⎢                                  ⎥⎟
⎜⎢0  ⅈ  0   -ⅈ⎥  ⎢   0        0     2 + 4⋅ⅈ     1   ⎥⎟
⎜⎢            ⎥  ⎢                                  ⎥⎟
⎝⎣0  1  0   1 ⎦  ⎣   0        0        0     2 + 4⋅ⅈ⎦⎠
@cbm755

This comment has been minimized.

Copy link
Contributor

cbm755 commented Apr 22, 2015

Thanks! are you going to do a PR or shall I?

@jksuom

This comment has been minimized.

Copy link
Member

jksuom commented Apr 22, 2015

It would suit me fine if you could do that. The comments would require some editing, at least to fix the typos. It seems that this paragraph could be removed altogether as the vectors added in the complex case are conjugates.

@cbm755 cbm755 referenced this issue Apr 24, 2015

Open

[WIP] dsolve fixes for systems #9235

5 of 13 tasks complete

@cbm755 cbm755 closed this in de08ee5 May 2, 2015

skirpichev added a commit to diofant/diofant that referenced this issue Jul 1, 2015

jordan_structure: use conjugate transpose, fixes complex eigenvalues
Add a test of a real matrix which needs this fix.

Fixes sympy/sympy#9274.

skirpichev added a commit to diofant/diofant that referenced this issue Jul 3, 2015

jordan_structure: use conjugate transpose, fixes complex eigenvalues
Add a test of a real matrix which needs this fix.

Fixes sympy/sympy#9274.

skirpichev added a commit to diofant/diofant that referenced this issue Jul 3, 2015

jordan_structure: use conjugate transpose, fixes complex eigenvalues
Add a test of a real matrix which needs this fix.

Fixes sympy/sympy#9274.

skirpichev added a commit to skirpichev/diofant that referenced this issue Jun 10, 2016

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment