Here is a small example. It surprised me at how slow it is.
import numpy as np
import sympy as sp
X = np.empty(shape=(10,3), dtype=object)
for i in range(10):
for j in range(3):
X[i,j] = sp.Symbol(f'X{i}{j}')
X = sp.Matrix(X.T @ X)
print(X.inv()) # This is really slow...
I would have thought that taking an inverse of a relatively small matrix would be fast. Am I missing something about the intended use case?
Here is your matrix:
In [43]: X
Out[43]:
⎡ 2 2 2 2 2 2 2 2 2 2 ⎤
⎢ X₀₀ + X₁₀ + X₂₀ + X₃₀ + X₄₀ + X₅₀ + X₆₀ + X₇₀ + X₈₀ + X₉₀ X₀₀⋅X₀₁ + X₁₀⋅X₁₁ + X₂₀⋅X₂₁ + X₃₀⋅X₃₁ + X₄₀⋅X₄₁ + X₅₀⋅X₅₁ + X₆₀⋅X₆₁ + X₇₀⋅X₇₁ + X₈₀⋅X₈₁ + X₉₀⋅X₉₁ X₀₀⋅X₀₂ + X₁₀⋅X₁₂ + X₂₀⋅X₂₂ + X₃₀⋅X₃₂ + X₄₀⋅X₄₂ + X₅₀⋅X₅₂ + X₆₀⋅X₆₂ + X₇₀⋅X₇₂ + X₈₀⋅X₈₂ + X₉₀⋅X₉₂⎥
⎢ ⎥
⎢ 2 2 2 2 2 2 2 2 2 2 ⎥
⎢X₀₀⋅X₀₁ + X₁₀⋅X₁₁ + X₂₀⋅X₂₁ + X₃₀⋅X₃₁ + X₄₀⋅X₄₁ + X₅₀⋅X₅₁ + X₆₀⋅X₆₁ + X₇₀⋅X₇₁ + X₈₀⋅X₈₁ + X₉₀⋅X₉₁ X₀₁ + X₁₁ + X₂₁ + X₃₁ + X₄₁ + X₅₁ + X₆₁ + X₇₁ + X₈₁ + X₉₁ X₀₁⋅X₀₂ + X₁₁⋅X₁₂ + X₂₁⋅X₂₂ + X₃₁⋅X₃₂ + X₄₁⋅X₄₂ + X₅₁⋅X₅₂ + X₆₁⋅X₆₂ + X₇₁⋅X₇₂ + X₈₁⋅X₈₂ + X₉₁⋅X₉₂⎥
⎢ ⎥
⎢ 2 2 2 2 2 2 2 2 2 2 ⎥
⎣X₀₀⋅X₀₂ + X₁₀⋅X₁₂ + X₂₀⋅X₂₂ + X₃₀⋅X₃₂ + X₄₀⋅X₄₂ + X₅₀⋅X₅₂ + X₆₀⋅X₆₂ + X₇₀⋅X₇₂ + X₈₀⋅X₈₂ + X₉₀⋅X₉₂ X₀₁⋅X₀₂ + X₁₁⋅X₁₂ + X₂₁⋅X₂₂ + X₃₁⋅X₃₂ + X₄₁⋅X₄₂ + X₅₁⋅X₅₂ + X₆₁⋅X₆₂ + X₇₁⋅X₇₂ + X₈₁⋅X₈₂ + X₉₁⋅X₉₂ X₀₂ + X₁₂ + X₂₂ + X₃₂ + X₄₂ + X₅₂ + X₆₂ + X₇₂ + X₈₂ + X₉₂ ⎦
Computing the symbolic inverse of a matrix is usually not a good idea unless the matrix is very small, sparse or has some special symmetry or something. If the matrix is fully symbolic and there is no reason for the expressions in the inverse to simplify then you should expect that symbolic expressions for the inverse will be very complicated. SymPy's inv method will try to simplify these for you because the expression is unlikely to be much use if it cannot be simplified. That simplification is also needed in order to determine whether the matrix actually is invertible. It is the simplification part that is slow here.
If you don't want the simplification and you are happy to assume that the matrix is invertible then you can calculate an expression for the inverse quickly using the standard inverse formula:
In [34]: %time Xinv = X.adjugate() / X.det()
CPU times: user 543 ms, sys: 957 µs, total: 544 ms
Wall time: 543 ms
Even here it is only taking 0.5 seconds because det attempts some simplification.
I won't show the result that is returned for Xinv because it is enormous. Instead I'll just show how long its string representation is:
In [38]: len(str(Xinv))
Out[38]: 618867
In [39]: len(str(Xinv[0,0]))
Out[39]: 68219
Even just the first element of this matrix has a string representation that is too long to be pasted into an answer on SO.