Jordan form hangs for a 3x3 matrix

[ferenci-tamas]

Minimal reproducible example:

A = Matrix([[-3, 1, 2], [1, -1, 0], [1, 0, -2]])
A.jordan_form()

This unfortunately hangs from some reason... (Despite the fact that the matrix is only 3x3, so this could be even solved analytically.) Wolfram Alpha gives a result.

[oscarbenjamin]

This can be improved in SymPy especially since the matrix is diagonalisable and SymPy can already compute that:

A.diagonalize()

The Jordan form calculation doesn't use the same techniques as the eigenvalue, eigenvector calculation.

The Wolfram Alpha result is only approximate and can also be given by SymPy:

In [25]: P, D = A.evalf().diagonalize()

In [26]: P.evalf(5)
Out[26]: 
⎡-0.87681  0.53495  0.35049⎤
⎢                          ⎥
⎢0.27278   0.79236  -0.7606⎥
⎢                          ⎥
⎣0.39597   0.31935  0.65004⎦

In [27]: D.evalf(5)
Out[27]: 
⎡-4.2143     0         0   ⎤
⎢                          ⎥
⎢   0     -0.32487     0   ⎥
⎢                          ⎥
⎣   0        0      -1.4608⎦

Unfortunately that doesn't work with the Jordan form method because internally it tries to remove the floats and replace them with exact rational numbers.

[ferenci-tamas]

Thank you very much @oscarbenjamin! This is perhaps a room for improvement itself (perhaps jordan_form could first check whether the matrix is diagonalizable? - it is somewhat strange that it can diagonalize the matrix, but it can't produce a Jordan normal form, while they're exactly identical in this case), but much more importantly: what to do if the matrix is not diagonalizable? Are there any general solutions, even if approximate...?

[oscarbenjamin]

The eigenvector code was significantly improved in gh-20614.

The Jordan normal form routine should be rewritten along the same lines as the eigenvects code.

[ferenci-tamas]

Looking forward to this improvement, thank you very much @oscarbenjamin!

[vprayag2005]

it can find Jordan block matrix

 >>> A = Matrix([[-3, 1, 2], [1, -1, 0], [1, 0, -2]])
>>> A.jordan_form(calc_transform=False)
>>>Matrix([[-2 - 4/((-1/2 - sqrt(3)*I/2)*(27 + 3*sqrt(111)*I)**(1/3)) - (-1/2 - sqrt(3)*I/2)*(27 + 3*sqrt(111)*I)**(1/3)/3, 0, 0], [0, -2 - (-1/2 + sqrt(3)*I/2)*(27 + 3*sqrt(111)*I)**(1/3)/3 - 4/((-1/2 + sqrt(3)*I/2)*(27 + 3*sqrt(111)*I)**(1/3)), 0], [0, 0, -2 - (27 + 3*sqrt(111)*I)**(1/3)/3 - 4/(27 + 3*sqrt(111)*I)**(1/3)]])

but it takes time to find M ie, P-1* J * P
i think to find p it is better to use eigenvects method since it diagonalizable matrix

[vprayag2005]

This line make it slower. ie nullspace() https://github.com/sympy/sympy/blob/9b7d671ea15131aab0fe80c8cd73ed8d0cca7be5/sympy/matrices/eigen.py#L1194-L1195

I used domainmatrix then compute nullpace().so it fixes. previously it hangs and noteven geting result.

--- a/sympy/matrices/eigen.py
+++ b/sympy/matrices/eigen.py
@@ -1187,6 +1187,7 @@ def pick_vec(small_basis, big_basis):
     # most matrices have distinct eigenvalues
     # and so are diagonalizable.  In this case, don't
     # do extra work!
+    jordan_basis=[]
     if len(eigs.keys()) == mat.cols:
         blocks     = sorted(eigs.keys(), key=default_sort_key)
         jordan_mat = mat.diag(*blocks)
@@ -1194,8 +1195,11 @@ def pick_vec(small_basis, big_basis):
         if not calc_transform:
             return restore_floats(jordan_mat)

-        jordan_basis = [eig_mat(eig, 1).nullspace()[0]
-                for eig in blocks]
+        for eig in blocks:
+            DOM = DomainMatrix.from_Matrix(eig_mat(eig, 1))
+            null_dom = DOM.nullspace()
+            null_mat = null_dom.to_Matrix().T
+            jordan_basis.append(null_mat)
         basis_mat    = mat.hstack(*jordan_basis)

         return restore_floats(basis_mat, jordan_mat)

Output

>>> P, J = A.jordan_form()
>>> P.evalf()
Matrix([
[-0.921622254378222 + 0.e-22*I, 1.35026174113329 - 0.e-23*I,  6.42863948675507 + 0.e-22*I],
[                          2.0,                         2.0,                         -2.0],
[ -1.70927535943692 - 0.e-25*I, 0.80606343352537 + 0.e-23*I, -2.90321192591155 + 0.e-23*I]])

>>> J.evalf()
Matrix([
[-1.46081112718911 + 0.e-22*I,                             0,                            0],
[                           0, -0.324869129433354 - 0.e-23*I,                            0],
[                           0,                             0, -4.21431974337754 + 0.e-21*I]])

[oscarbenjamin]

The whole calculation should be rewritten to use only DomainMatrix.

[vprayag2005]

As per https://github.com/sympy/sympy/issues/28097#issuecomment-2916610664 We can use for rational coefficients. If it is not rational coefficients we can use exiting techniques with domain matrix right?

[oscarbenjamin]

If the coefficients are not rational then it is not necessarily possible to express the roots using RootOf. For algebraic coefficients we can sometimes although it can be slow:

In [5]: p = sqrt(2)*x**2 + x + 1

In [6]: all_roots(p, extension=True)
Out[6]: 
⎡       ⎛   4    2             ⎞         ⎛   4    2             ⎞⎤
⎣CRootOf⎝2⋅x  - x  - 2⋅x - 1, 2⎠, CRootOf⎝2⋅x  - x  - 2⋅x - 1, 3⎠⎦

The eigenvects code uses FiniteMonogenicExtension for domains that are not rational coefficients and that can also be used here.

[oscarbenjamin]

I think it is better to focus on providing an implementation for rationals first though.

[vprayag2005]

This is a rough implementaion for rational.

@@ -1060,10 +1060,13 @@ def _jordan_form(M, calc_transform=True, *, chop=False):

     jordan_block
     """
-
+    from sympy import all_roots, eye, AlgebraicNumber
+    from collections import Counter
     if not M.is_square:
         raise NonSquareMatrixError("Only square matrices have Jordan forms")

+    a = None
+    eigenvalues={}
     mat        = M
     has_floats = M.has(Float)

@@ -1099,12 +1102,9 @@ def eig_mat(val, pow):

         if (val, pow) in mat_cache:
             return mat_cache[(val, pow)]
-
-        if (val, pow - 1) in mat_cache:
-            mat_cache[(val, pow)] = mat_cache[(val, pow - 1)].multiply(
-                    mat_cache[(val, 1)], dotprodsimp=None)
         else:
-            mat_cache[(val, pow)] = (mat - val*M.eye(M.rows)).pow(pow)
+            eigen_scalar = a if a is not None else val
+            mat_cache[(val, pow)] = (mat - eigen_scalar*eye(M.rows)).to_DM(extension=True).pow(pow)

         return mat_cache[(val, pow)]

@@ -1117,7 +1117,10 @@ def nullity_chain(val, algebraic_multiplicity):
         # so use the rank-nullity theorem
         cols    = M.cols
         ret     = [0]
-        nullity = cols - eig_mat(val, 1).rank()
+        mat_expr = eig_mat(val, 1).to_Matrix()
+        if a:
+            mat_expr = mat_expr.subs(a, val)
+        nullity = cols - mat_expr.rank()
         i       = 2

         while nullity != ret[-1]:
@@ -1126,7 +1129,10 @@ def nullity_chain(val, algebraic_multiplicity):
             if nullity == algebraic_multiplicity:
                 break

-            nullity  = cols - eig_mat(val, i).rank()
+            mat_expr = eig_mat(val, i).to_Matrix()
+            if a:
+                mat_expr = mat_expr.subs(a, val)
+            nullity = cols - mat_expr.rank()
             i       += 1

             # Due to issues like #7146 and #15872, SymPy sometimes
@@ -1173,37 +1179,44 @@ def pick_vec(small_basis, big_basis):
     if has_floats:
         from sympy.simplify import nsimplify
         mat = mat.applyfunc(lambda x: nsimplify(x, rational=True))
-
-    # first calculate the jordan block structure
-    eigs = mat.eigenvals()
-
-    # Make sure that we have all roots in radical form
-    for x in eigs:
-        if x.has(CRootOf):
-            raise MatrixError(
-                "Jordan normal form is not implemented if the matrix have "
-                "eigenvalues in CRootOf form")
+
+    # check if matrix has rational coefficients
+    has_rationals = all(x.is_rational for x in mat)

     # most matrices have distinct eigenvalues
     # and so are diagonalizable.  In this case, don't
     # do extra work!
-    if len(eigs.keys()) == mat.cols:
-        blocks     = sorted(eigs.keys(), key=default_sort_key)
+
+    if has_rationals:
+        char_eq = mat.charpoly().as_expr()
+        roots = all_roots(char_eq)
+        eigenvalues = dict(Counter(roots))
+        for i in eigenvalues.keys():
+            if not isinstance(i,int) and a is not None:
+                a = AlgebraicNumber(i, alias='a')
+                break
+
+
+    if len(eigenvalues.keys()) == mat.cols:
+        jordan_basis = []
+        blocks     = sorted(eigenvalues.keys(),key=default_sort_key)
         jordan_mat = mat.diag(*blocks)

         if not calc_transform:
             return restore_floats(jordan_mat)

-        jordan_basis = [eig_mat(eig, 1).nullspace()[0]
-                for eig in blocks]
+        for eig in blocks:
+            if a is not None:
+                jordan_basis.append(eig_mat(eig, 1).nullspace().to_Matrix().subs(a, eig).T)
+            else:
+                jordan_basis.append(eig_mat(eig, 1).nullspace().to_Matrix().T)
         basis_mat    = mat.hstack(*jordan_basis)

         return restore_floats(basis_mat, jordan_mat)

     block_structure = []
-
-    for eig in sorted(eigs.keys(), key=default_sort_key):
-        algebraic_multiplicity = eigs[eig]
+    for eig in sorted(eigenvalues.keys(), key=default_sort_key):
+        algebraic_multiplicity = eigenvalues[eig]
         chain = nullity_chain(eig, algebraic_multiplicity)
         block_sizes = blocks_from_nullity_chain(chain)

@@ -1245,23 +1258,33 @@ def pick_vec(small_basis, big_basis):
     # are linearly independent.
     jordan_basis = []

-    for eig in sorted(eigs.keys(), key=default_sort_key):
+    for eig in sorted(eigenvalues.keys(), key=default_sort_key):
         eig_basis = []

         for block_eig, size in block_structure:
             if block_eig != eig:
                 continue
+            big_vec = eig_mat(eig, size).nullspace().to_Matrix().T
+            small_vec = eig_mat(eig, size - 1).nullspace().to_Matrix().T
+
+            if a:
+                big_vec = big_vec.subs(a, eig)
+                small_vec = small_vec.subs(a, eig)
+
+            null_big = [big_vec.col(i) for i in range(big_vec.cols)]
+            null_small = [small_vec.col(i) for i in range(small_vec.cols)]

-            null_big   = (eig_mat(eig, size)).nullspace()
-            null_small = (eig_mat(eig, size - 1)).nullspace()

             # we want to pick something that is in the big basis
             # and not the small, but also something that is independent
             # of any other generalized eigenvectors from a different
             # generalized eigenspace sharing the same eigenvalue.
             vec      = pick_vec(null_small + eig_basis, null_big)
-            new_vecs = [eig_mat(eig, i).multiply(vec, dotprodsimp=None)
-                    for i in range(size)]
+            new_vecs = [
+                eig_mat(eig, i).to_Matrix().subs(a, eig) if a else eig_mat(eig, i).to_Matrix()
+                .multiply(vec, dotprodsimp=None)
+                for i in range(size)
+            ]

             eig_basis.extend(new_vecs)
             jordan_basis.extend(reversed(new_vecs))

[oscarbenjamin]

That looks good, although I would write a separate function rather than modifying the existing one.

[vprayag2005]

This is the seperate function

    def eig_nullspace_matrix(eig, pow, nullspace = True):
        char_mat = eig_mat(eig, pow)
        char_mat = char_mat.nullspace().to_Matrix() if nullspace else char_mat.to_Matrix()
        if a is not None:
            char_mat = char_mat.subs(a, eig)
        return char_mat.T if nullspace else char_mat

and modified eig_mat function

@@ -1099,12 +1102,11 @@ def eig_mat(val, pow):

         if (val, pow) in mat_cache:
             return mat_cache[(val, pow)]
-
         if (val, pow - 1) in mat_cache:
-            mat_cache[(val, pow)] = mat_cache[(val, pow - 1)].multiply(
-                    mat_cache[(val, 1)], dotprodsimp=None)
+            mat_cache[(val, pow)] = mat_cache[(val, pow - 1)] * mat_cache[(val, 1)]
         else:
-            mat_cache[(val, pow)] = (mat - val*M.eye(M.rows)).pow(pow)
+            eigen_scalar = a if a is not None else val
+            mat_cache[(val, pow)] = (mat - eigen_scalar*eye(M.rows)).to_DM(extension=True).pow(pow)

         return mat_cache[(val, pow)]

We can replace eig_mat where it calls with this new function

[vprayag2005]

The above code works fine with rational coefficients

The eigenvects code uses FiniteMonogenicExtension for domains that are not rational coefficients and that can also be used here.

I think eigenvects code uses FiniteMonogenicExtension only when domain are not EX.

[vprayag2005]

If the coefficients are not rational, should we follow the dom_eigenvects code. But it is applicable only if the domain is not EX. If is EX domain we have to use existing techniques. We can’t use DomainMatrix because the current design can’t find the nullspace with DomainMatrix, so we have to rely on existing methods.

[oscarbenjamin]

We can’t use DomainMatrix because the current design can’t find the nullspace with DomainMatrix

I don't understand what you mean:

In [1]: M = randMatrix(3)

In [2]: r1, r2, r3 = M.charpoly().all_roots()

In [3]: K = QQ.algebraic_field(r1)

In [5]: ((r1*eye(3)).to_DM(domain=K) - M.to_DM()).nullspace()
Out[5]: DomainMatrix({0: {2: ANP([1, -104, 333], [1, -119, -9387, -145125], QQ), 0: ANP([87, -5736], [1, -119, -9387, -145125], QQ), 1: ANP([99, 6099], [1, -119, -9387, -145125], QQ)}}, (1, 3), QQ<CRootOf(x**3 - 119*x**2 - 9387*x - 145125, 0)>)

You have to use a domain that understands the algebraic constraint on the eigenvalue r1 i.e. K here.

[vprayag2005]

I don't understand what you mean:

What I mean is that the existing method using DomainMatrix can’t be used to compute the nullspace if domain is EX.

>>> M
Matrix([
[sqrt(2), pi],
[sqrt(3),  E]])
>>> val = sorted(M.eigenvals().keys(),key=default_sort_key)[0]
>>> val
sqrt(2)/2 + E/2 + sqrt(-2*sqrt(2)*E + 2 + exp(2) + 4*sqrt(3)*pi)/2
>>> (M - val*eye(M.rows)).to_DM(extension=True)
DomainMatrix({0: {0: -sqrt(-2*sqrt(2)*E + 2 + exp(2) + 4*sqrt(3)*pi)/2 - E/2 + sqrt(2)/2, 1: pi}, 1: {1: -sqrt(-2*sqrt(2)*E + 2 + exp(2) + 4*sqrt(3)*pi)/2 - sqrt(2)/2 + E/2, 0: sqrt(3)}}, (2, 2), EXRAW)
>>> (M - val*eye(M.rows)).to_DM(extension=True).nullspace()
TypeError: cannot determine truth value of Relational: -sqrt(3)*pi + (-sqrt(-2*sqrt(2)*E + 2 + exp(2) + 4*sqrt(3)*pi)/2 - sqrt(2)/2 + E/2)*(-sqrt(-2*sqrt(2)*E + 2 + exp(2) + 4*sqrt(3)*pi)/2 - E/2 + sqrt(2)/2) < 0

In [2]: r1, r2, r3 = M.charpoly().all_roots()

If the domain is EX, we can’t use all_roots(), right?

NotImplementedError: sorted roots not supported over EX

[oscarbenjamin]

TypeError: cannot determine truth value of Relational: -sqrt(3)pi + (-sqrt(-2sqrt(2)E + 2 + exp(2) + 4sqrt(3)pi)/2 - sqrt(2)/2 + E/2)(-sqrt(-2*sqrt(2)E + 2 + exp(2) + 4sqrt(3)*pi)/2 - E/2 + sqrt(2)/2) < 0

That is a bug that should be fixed.

I think that you are overcomplicating things by trying to consider all domains at once. It is best just to start by implementing something that is only for matrices of rational numbers.

[vprayag2005]

Actually this is what i implemented for rational matrix https://github.com/sympy/sympy/issues/24942#issuecomment-3007548748

[oscarbenjamin]

I suggest opening a pull request that can use DomainMatrix for the case of rational coefficients only. I think that it would be better to write new separate functions rather than trying to adapt the existing functions. Then there can be a check to see if the matrix is all rational numbers and if so use the DomainMatrix calculation.