transformation matrix for variable matrix of given jordan type

savecancel

We can set up the problem as follows:

n = 3
jordan_type = (2,1)

assert sum(jordan_type) == n
import itertools
A_coeff_names = {(i,j) : 'a_{}_{}'.format(i,j) for (i,j) in itertools.combinations(range(n),2)}
R = PolynomialRing(QQ, names=A_coeff_names.values())
A_coeff = {idx : R(A_coeff_names[idx]) for idx in A_coeff_names}
A = Matrix(R, n, A_coeff)
J = block_diagonal_matrix([jordan_block(0, s) for s in jordan_type], subdivide=False)

A necessary condition for $gAg^{-1} = J$ is that the rank of $A^k$ equals the rank of $J^k$ for $1 \leq k < n$.

(And since $A$ has all zeros as eigenvalues, this will give the right Jordan type.)

In particular it is necessary that the $(\operatorname{rank}(J^k)+1)$-minors of $A^k$ vanish; these are polynomial equations upon the coefficients of $A$:

A_coeff_eqns = []
J_power = identity_matrix(QQ, n)
A_power = identity_matrix(R, n)
for k in range(n):
    r = J_power.rank()
    # need (r+1) x (r+1) minors to vanish
    A_coeff_eqns.extend([eqn for eqn in A_power.minors(r+1) if eqn != 0])
    J_power *= J
    A_power *= A

We can solve these equations symbolically, substitute a solution into $A$, change to an exact ring (the fraction field of the polynomial ring in the new variables), calculate the Jordan form (to be sure), and change the names of the variables back to the originals where possible:

A_coeff_symb = map(SR, A_coeff_names.values())
A_coeff_eqns_symb = map(SR, A_coeff_eqns)
for sol in solve(A_coeff_eqns_symb, A_coeff_symb):
    # substitute symbolic solution
    A_new = A.change_ring(SR).apply_map(lambda x: x.subs(sol))
    A_new_vars = list(set(sum([list(x.variables()) for x in A_new.list()], [])))
    # change ring to fraction field of polynomial ring in the new variables
    A_new_ring = PolynomialRing(QQ, names=A_new_vars).fraction_field()
    A_new_rat = A_new.change_ring(A_new_ring)
    if A_new_rat.jordan_form(subdivide=False) != J:
        continue
    # change names of variables back to the originals where possible
    A_final_subs = {A_new_rat[i,j] : A_coeff[(i,j)] for (i,j) in A_coeff if A_new_rat[i,j] in A_new_vars}
    A_final = A_new_rat.apply_map(lambda x: x.subs(A_final_subs))
    print A_final
    print 
    print A_final.jordan_form(subdivide=False, transformation=True)[1]
    print
    print '---'
    print

Output:

[    0 a_0_1 a_0_2]
[    0     0     0]
[    0     0     0]

[         a_0_1              0              0]
[             0              1              1]
[             0              0 a_0_1/(-a_0_2)]

---

[    0     0 a_0_2]
[    0     0 a_1_2]
[    0     0     0]

[a_0_2     0     1]
[a_1_2     0     0]
[    0     1     0]

---

Beware that this assumes Sage can solve the polynomial equations symbolically, which may not be true for large n. To give Sage an easier time (though the solutions may become less pretty), you can replace the list of equations by a Groebner basis for the same ideal:

A_coeff_eqns = R.ideal(A_coeff_eqns).groebner_basis()