Sagemath: Solving linear equation symbolically takes a lot of time

Question

I have a set of polynomials $\lbrace p_{ij}|i,j\in\lbrace 1,...,n\rbrace \rbrace$ with $p_{ij}\in \mathbb{Z}[x]$ as well as degree $n$ and a given constant real vector $b\in\mathbb{R}^n$. The polynomials can be seen as a matrix $P$ and I want to solve the equation $Pz=b$ for $n=20$.

I wanted to use sagemath since inverting $20\times 20$ matrices by hand is somewhat tedious. But it turns out that the commands solve and P.solve_right(b) take a lot of time for this kind of problem and did not terminate, finally. Inverting $P$ by P.inverse() did not terminate either.

Why could this happen and what could I do to speed it up? The size is not enormous even against the benchmark for runtime for matrix inversion in $n^3$.

EDIT: Additional information:

$Pz$ is simply the product of $P$ with an element $z$ which gives a constant. $z$ is in my case the restriction of a $n$ dimensional vectorfield $Z$ and you could unterstand this analogously to characteristics for transport equations, $z$ being defined as $Z$ along the characteristic. Each component of $z$ will be a rational function in $x$.

Answer 1

Revised answer

By playing around in SageCell I've discovered that SageMath can calculate the determinant of a $\ 20\times20\ $ matrix over $\ \mathbb{Z}[x]\ $ with moderately sized entries of degree $\ 20\ $ fairly quickly. This means that you should be able to use Cramer's rule to solve your equation in SageMath. Assuming the entries of $\ b\ $ are rational numbers, first multiply each row $\ i\ $ of the equation $$ Pz=b\tag{1}\label{e1} $$ by the denominator of $\ b_i\ $ to obtain an equivalent equation $$ Rz=\beta $$ in which all entries of the right-hand side are now integers. Cramer's rule tells us that the solution of this equation is given by $$ z_j=\frac{\det Q(j)}{\det R} $$ where $\ Q(j)\ $ is the matrix obtained from $\ R\ $ by replacing its $\ j^\text{th}\ $ column with $\ \beta\ .$ Since all the entries of $\ R\ $ and $\ Q(j)\ $ are in $\ \mathbb{Z}[x]\ ,$ SageMath can compute their determinants $\ d=\det R\ $ and $\ \zeta_j=\det Q(j)\ $ reasonably quickly. The solution to your equation is then given by $$ z=\left(\frac{1}{d}\right)\zeta\ . $$

The following SageMath code implements this method to solve equation (\ref{e1}) for randomly generated $\ P\ $ and $\ b\ .$ Hopefully, you will be able to adapt it to solve your problem.

n=20
Id=matrix(ZZ['x'],n,n)
P=matrix(ZZ['x'],n,n)
R=matrix(ZZ['x'],n,n)
ff=FractionField(ZZ['x'])
b=vector(QQ,n)
beta=vector(QQ,n)
z=vector(ff,n)
for i in range(n):
    Id[i,i]=1
    
# Generate random P and b
for i in range(n):
    p=ZZ.random_element(999)
    q=ZZ.random_element(998)+1
    b[i]=p/q
    for j in range(n):
        P[i,j]=(ZZ.random_element(21)-10)
        for k in range(n):
            P[i,j]=P[i,j]*x+(ZZ.random_element(21)-10)
    
# The following code solves the equation P*z=b for z
# First multiply each row i of the equation by the denominator of b[i]
for i in range(n):
    for j in range(n):
        R[i,j]=P[i,j]*b[i].denominator()
    beta[i]=b[i].numerator()
        
d=R.determinant()
for j in range(n):
    Q=R*Id
    for i in range(n):
        Q[i,j]=beta[i]
    z[j]=Q.determinant()/d

# Check that proper solution has been found
print(P*z-b)