Revision history [back]

The solution for me was to remove the free variables from solve:

def sagemath_solve(augmented_matrix):

    A = augmented_matrix[:, :-1]
    Y = augmented_matrix[:, -1]

    m, n = A.dimensions()
    p, q = Y.dimensions()

    if m!=p:
        raise RuntimeError("The matrices have different numbers of rows")
    X = vector([var("x_{}".format(i)) for i in [1..n]])

    # don't include the free variables in solve
    X_pivots = vector(X[:max(A.pivots())+1])

    sols = []
    for j in range(q):
        system = [A[i]*X==Y[i,j] for i in range(m)]
        sols += solve(system, *X_pivots)
    return sols


A = matrix([
    [1,1,1,-1,1],
    [0,1,-1,1,-1],
    [3,0,6,-6,6],
    [0,-1,1,-1,1]
])
print(sagemath_solve(A))

# [[
#   x_1 == -2*x_3 + 2*x_4 + 2, 
#   x_2 == x_3 - x_4 - 1
# ]]

The solution for me was to remove the free variables from solve:solve.

(I also fixed the vars to start from x_0).

def sagemath_solve(augmented_matrix):

    A = augmented_matrix[:, :-1]
    Y = augmented_matrix[:, -1]

    m, n = A.dimensions()
    p, q = Y.dimensions()

    if m!=p:
        raise RuntimeError("The matrices have different numbers of rows")
    X = vector([var("x_{}".format(i)) for i in [1..n]])
[0..n-1]])

    # don't include the free variables in solve
    X_pivots = vector(X[:max(A.pivots())+1])

    sols = []
    for j in range(q):
        system = [A[i]*X==Y[i,j] for i in range(m)]
        sols += solve(system, *X_pivots)
    return sols


A = matrix([
    [1,1,1,-1,1],
    [0,1,-1,1,-1],
    [3,0,6,-6,6],
    [0,-1,1,-1,1]
])
print(sagemath_solve(A))

# [[
#   x_0 == -2*x_2 + 2*x_3 + 2, 
#   x_1 == -2*x_3 + 2*x_4 + 2, 
#   x_2 == - x_3 - x_4 - 1
# ]]

The solution for me was to remove the free variables from solve.

(I also fixed the vars to start from x_0).

def sagemath_solve(augmented_matrix):

    A = augmented_matrix[:, :-1]
    Y = augmented_matrix[:, -1]

    m, n = A.dimensions()
    p, q = Y.dimensions()

    if m!=p:
        raise RuntimeError("The matrices have different numbers of rows")
    X = vector([var("x_{}".format(i)) for i in [0..n-1]])

    # don't include the free variables in solve
    X_pivots = vector(X[:max(A.pivots())+1])
vector([var("x_{}".format(i)) for i in [0..n-1] if i in A.pivots()])

    sols = []
    for j in range(q):
        system = [A[i]*X==Y[i,j] for i in range(m)]
        sols += solve(system, *X_pivots)
    return sols


A = matrix([
    [1,1,1,-1,1],
    [0,1,-1,1,-1],
    [3,0,6,-6,6],
    [0,-1,1,-1,1]
])
print(sagemath_solve(A))

# [[
#   x_0 == -2*x_2 + 2*x_3 + 2, 
#   x_1 == x_2 - x_3 - 1
# ]]