Use SymPy.solve instead of generic lu for A\b

[olof3]

@vars a1 a2

n = 7

# Just a made-up example
A = diagm(0 => ones(Int, n), 1 => fill(a1, n-1), 2 => fill(1, n-2), -1 => fill(a2, n-1))
b = Vector{Sym}(1:n)

x = A \ b

This doesn't return on my machine, at least not before I got tired of waiting. Seems to work for smaller n though, but the solution is very long and not simplified. The problems get worse if more symbolic variables are included. The generic lu that is called by \ doesn't simplify intermediate expressions, which could be one reason for what seems like exponential complexity blow up.

SymPy.solve seems to work better. I guess one should really use linsolve, but I couldn't get this to work?! I guess one has to think about what to return, or if to throw an error, if there is no / multiple solutions. However, my feeling is that anything would be better than what is currently done by the generic lu.

function Base.:\(A::AbstractMatrix{Sym}, b::AbstractVecOrMat)
    n = size(A, 2)
    X = [symbols(prod(("x$k," for k=1:n)))...] # FIXME: This could almost surely be done in a better way
    sol = solve(A*X - b, X) # FIXME: should use linsolve

    if length(sol) == 0 # FIXME: Add better error handling
        error("No solution found")
    elseif length(sol) < n
        error("Multiple solutions found")
    end
    [sol[x] for x in X]
end

[jverzani]

This seems like a job for linsolve. I thought it great how generic julia algorithms can be utilized, but having it stall out so long is definitely something to work around. I think the proper fix is something like

function  \(A,b) 
    M = Sym.([A -b])
    m,n = size(M)
    linsolve(M,  NTuple{m,Sym}("x_i"  for   i in  1:m))
end

I should have time this week. Thanks for all the feedback.

[jverzani]

Thanks again. I ended using your approach in #356.

[olof3]

Thanks a lot for sorting things out!