Python Sympy : Check equality or inequality between 2 projected matrices

I reformulate more precisely this post. I try to check the equality or the inequality between 2 matrices. Each of these 2 matrices is computed slightly differently.

Actually, this is the computation of changing parameters between initial parameters for each row/colum and final parameters for the final matrix. That's why in both computations, I am using the Jacobian J formulating the derivatives between initial and final parameters :

The formula is : F_final = J^T F_initial J

The first matrix has size 31x31 and the second one has 32x32 size.

The Jacobian applied on the first matrix is 31x7 and on the second one is 32x8.

My issue is that I want to check if the 2 matrix projected 8x8 and computed by the formula :

F_final = J^T F_initial J

is equal to the matrix projected 7x7 by removing the 3th row/column of the final matrix 8x8.

I have coded the following Python Sympy script to check that :

import os, sys
import numpy as np
from sympy import *
from sympy import symbols, Matrix, Symbol, Transpose, eye, solvers


# Big_31 Fisher : symmetry by inversinf (i,j) in min(i,j) max(i,j)
FISH_Big_1_SYM = Matrix([[Symbol(f'sp_{min(i,j)}_{max(j,i)}') for i in range(1, 31+1)] for j in range(1, 31+1)])

# Big_32 Fisher : symmetry by inversinf (i,j) in min(i,j) max(i,j)
FISH_Big_2_SYM = Matrix([[Symbol(f'sp_{min(i,j)}_{max(j,i)}') for i in range(1, 32+1)] for j in range(1, 32+1)])

# Introduce neutrino
NEU_row_SYM = Matrix([[Symbol(f'neu_{min(i,j)}_{max(i,j)}') for i in range(1,2)] for j in range(1, 32+1)])
NEU_col_SYM = Matrix([[Symbol(f'neu_{min(i,j)}_{max(i,j)}') for i in range(1, 32+1)] for j in range(1, 2)])

# Introduce neutrino
J_NEU_row_SYM = Matrix([[Symbol(f'j_neu_{min(i,j)}_{max(i,j)}') for i in range(1, 8+1)] for j in range(1, 32+1)])
J_NEU_col_SYM = Matrix([[Symbol(f'j_neu_{min(i,j)}_{max(i,j)}') for i in range(1, 32+1)] for j in range(1, 8+1)])

# Temporary matrix to build FISH_Big_2_SYM 
mat_temp = Matrix([[Symbol(f'sp_{min(i,j)}_{max(j,i)}') for i in range(1,32+1)] for j in range(1, 32+1)])

mat_temp[0:2,:] = FISH_Big_2_SYM[0:2,:]
mat_temp[:,0:2] = FISH_Big_2_SYM[:,0:2]
mat_temp[2,:] = NEU_col_SYM
mat_temp[:,2] = NEU_row_SYM
mat_temp[3:32,3:32] = FISH_Big_2_SYM[2:31,2:31]

# Copy built FISH_Big_2_SYM
FISH_Big_2_SYM = np.copy(mat_temp)

# Jacobian 1 : not symmetric
J_1_SYM = Matrix([[Symbol(f'j_{min(i,j)}_{max(i,j)}') for i in range(1, 7+1)] for j in range(1, 31+1)])
print('shape_1', np.shape(J_1_SYM))

# Jacobian 2 : not symmetric
J_2_SYM = Matrix([[Symbol(f'j_{min(i,j)}_{max(i,j)}') for i in range(1, 8+1)] for j in range(1, 32+1)])
print('shape_2', np.shape(J_2_SYM))

# Jacobian 2 : row and column supplementary
J_2_row_SYM = Matrix([[Symbol(f'j_{min(i,j)}_{max(i,j)}') for i in range(1,8+1)] for j in range(1, 32+1)])
J_2_col_SYM = Matrix([[Symbol(f'j_{min(i,j)}_{max(i,j)}') for i in range(1, 32+1)] for j in range(1, 8+1)])

# Temporary matrix to build J_2_SYM 
j_temp = np.copy(J_2_SYM)

# Add row/col into J_2_SYM
j_temp = np.insert(j_temp, 2, J_NEU_row_SYM[0,:], axis=0)
j_temp = np.insert(j_temp, 2, J_NEU_col_SYM[:,0], axis=1)

# Copy built J_2_SYM
J_2_SYM = np.copy(j_temp)

print('Big_shape_1', np.shape(FISH_Big_1_SYM))
print('Big_shape_2', np.shape(FISH_Big_2_SYM))

# Projection 31x31
FISH_proj_1 = np.dot(np.dot(J_1_SYM.T,FISH_Big_1_SYM),J_1_SYM)

# Projection 32x32
FISH_proj_2 = np.dot(np.dot(J_2_SYM.T,FISH_Big_2_SYM),J_2_SYM)

# Test equality between 2 matrices
print('compare = ', FISH_proj_1.compare(FISH_proj_2))

But whereas the first matricial product seems to be valid ( FISH_proj_1 = np.dot(np.dot(J_1_SYM.T,FISH_Big_1_SYM),J_1_SYM)), the second one FISH_proj_2 = np.dot(np.dot(J_2_SYM.T,FISH_Big_2_SYM),J_2_SYM) generates an error :

shape_1 (31, 7)
shape_2 (32, 8)
Big_shape_1 (31, 31)
Big_shape_2 (32, 32)
Traceback (most recent call last):
  File "demo_projection_sympy.py", line 62, in <module>
    FISH_proj_2 = np.dot(np.dot(J_2_SYM.T,FISH_Big_2_SYM),J_2_SYM)
  File "<__array_function__ internals>", line 6, in dot
ValueError: shapes (9,33) and (32,32) not aligned: 33 (dim 1) != 32 (dim 0)

I tried to do many tests but impossible to find the error.

If anyone coud see at first sight what's wrong. The dimensions seems to be correct.

from Python Sympy : Check equality or inequality between 2 projected matrices