RFC on solve being used to solve for undetermined coefficients

[Chris Smith]

Nearly from the beginning, `solve` has recognized the `undetermined coefficients` case wherein a single equation is solved for a set of symbols which will simultaneously set it equal to zero. This is different than the normal use of `solve`, however, and might be considered a bug or feature. As an example, consider:

[in1] solve(a*x + b - (2*x + 4), (a, b))  # solve for a and b
{a: 2, b: 4}
[in2] solve(a*x + b - (2*x + 4), exclude=[x])  # solve for as many as possible, but not x
{a: 2, b: 4}
[in3] solve(a*x + b - (2*x + 4), a)  # solve for a
[(-b + 2*x + 4)/x]

[in4] solve(a*x + b - (2*x + 4))  # solve for anything and tell me what you did

[{a: (-b + 2*x + 4)/x}]

The "undetermined coefficients" case is returned from inputs 1 and 2 while the more usual "solution" is returned from 3 and 4.

The request for comment is whether this feature should be removed from `solve` so a better defined equation set is provided rather than building that equation set automatically when this case is detected as one-equation-many-symbols. 

To me it feels natural to get the undetermined coefficients solution when I ask for many symbols from one equation. If I want sympy to pick any variable other than some I would use `exclude=some` (where `some` is a list or set of things to ignore). But not everyone feels this way and refusing the temptation to guess is important, too.

Since this behavior is described in the docstring and has been present for many years it feels like a regression of this feature to me (though to others it feels like a bug fix). If you use this feature and/or have comments, please let us know your thoughts on this: 

1) should it be kept?

2) if removed, should there be deprecation?

Thanks,

 Chris

[Aaron Meurer]

How does it decide when to take coefficients? I tried this

>>> solve(a*cos(x) + b*sin(x), [a, b])
[(-b*tan(x), b)]

I would expect an undetermined coefficients solution to be [(0, 0)].
OTOH, it would need to be careful about this. This is only valid if
the terms in question are linearly independent (this can be checked by
seeing if the Wronskian is nonzero).

Personally, I would prefer if solve() were more explicit and did less
guessing about what the user wants. The problem with guessing is that
if you do know what you want, it becomes very difficult to tell
solve() to do that specific thing. It's especially problematic that
the actual behavior of solve depends on seemingly insignificant things
like how the symbols are specified. Just as a simple example
(unrelated to undetermined coefficients), say I have an
underdetermined system like solve([x + y + z, x - y - z], [x, y, z]).
I would expect specifying [x, y, z] to mean that I don't want the
solutions to be in terms of any of those variables (even preferring an
error to that), but I can't figure out how to make it do that.

There is a function solve_undetermined_coefficients, although it
currently only supports polynomials. Its behavior seems to be more
well-specified for problems like this. One can also extract the system
manually with collect(), although that's a bit advanced for the
average user.

It would be better for the solvers to follow the Unix philosophy of
"each function does one thing and does it well". So there should be
one function just for solving a function algebraically, one function
for solving a system, one for solving undetermined coefficients, and
so on. For a function like simplify(), there is a benefit to the "just
do something, even if it might not be exactly what the user wants"
philosophy. But I don't see the benefit of that for solve(). If you
have an equation or system of equations, you want to solve it in a
very specific way, and solving it in some other way is almost always
unhelpful (or at least that's been my experience).

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/b3553af5-0df2-49cb-bfaa-b3c354429864n%40googlegroups.com.

[Chris Smith]

Do you have any feelings about making the change without deprecation?

[Aaron Meurer]

How long has the feature been there?

Aaron Meurer

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/beaefb88-28cb-4661-ae09-a17778bf0bf0n%40googlegroups.com.

[Chris Smith]

since 2011 (2ff42e6).

[Aaron Meurer]

Ah I missed your "nearly from the beginning" bit. I would say it
should be deprecated, although if there are other kinds of cleanups
that can be done with solve, it might make sense to do them all
together (but this could be a big project).

Aaron Meurer

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/826aaa9c-19c3-49e9-95ac-70f71cc4a76fn%40googlegroups.com.

[Chris Smith]

Perhaps we could use an interim keyword `coefficients` which defaults to True. If an undetermined coefficients case is detected, a warning can be raised and a message given before returning that solution which says that if a non-coefficients solution is desired (during deprecation period) the keyword `coefficients=False` can be used. Once the deprecation period is over, the `coefficients` keyword and the recognition of undetermined coefficients can be removed.

/c

[Oscar Benjamin]

On Thu, 30 Jun 2022 at 23:27, Chris Smith <smi...@gmail.com> wrote:
>
> Nearly from the beginning, `solve` has recognized the `undetermined coefficients` case wherein a single equation is solved for a set of symbols which will simultaneously set it equal to zero. This is different than the normal use of `solve`, however, and might be considered a bug or feature. As an example, consider:
>
> [in1] solve(a*x + b - (2*x + 4), (a, b)) # solve for a and b
> {a: 2, b: 4}
> [in2] solve(a*x + b - (2*x + 4), exclude=[x]) # solve for as many as possible, but not x
> {a: 2, b: 4}
> [in3] solve(a*x + b - (2*x + 4), a) # solve for a
> [(-b + 2*x + 4)/x]
> [in4] solve(a*x + b - (2*x + 4)) # solve for anything and tell me what you did
> [{a: (-b + 2*x + 4)/x}]
>
> The "undetermined coefficients" case is returned from inputs 1 and 2 while the more usual "solution" is returned from 3 and 4.
>
> The request for comment is whether this feature should be removed from `solve` so a better defined equation set is provided rather than building that equation set automatically when this case is detected as one-equation-many-symbols.
>
> To me it feels natural to get the undetermined coefficients solution when I ask for many symbols from one equation. If I want sympy to pick any variable other than some I would use `exclude=some` (where `some` is a list or set of things to ignore). But not everyone feels this way and refusing the temptation to guess is important, too.

It does not feel natural to me at all that solve should arbitrarily
decide to get an undetermined coefficients solution. Also solve is
just flat out inconsistent:

In [7]: a, b, x = symbols('a, b, x')

In [8]: eq = a + x*b

In [9]: solve([eq], [a, b])
Out[9]: {a: -b⋅x}

In [10]: nonlinsolve([eq], [a, b])
Out[10]: {(-b⋅x, b)}

In [11]: linsolve([eq], [a, b])
Out[11]: {(-b⋅x, b)}

Okay so far so good but what happens to solve if I change [eq] to eq:

In [12]: solve([eq], [a, b])
Out[12]: {a: -b⋅x}

In [13]: solve(eq, [a, b])
Out[13]: {a: 0, b: 0}

Firstly there should be no distinction between eq and [eq] here.
Secondly in Out[13] solve has arbitrarily decided to completely
reinterpret the problem that was set. The problem is: compute the set
of solutions for a and b to the equation system involving a and b that
is parametrised by a symbol x and give expressions for the solutions
that are valid for almost all possible values of the parameter x.

That is what {a: -b*x} is. The return value {a:0, b:0} implies that
there is a unique solution for a and b when it is clear that there is
not a unique solution for this underdetermined system. If you
substitute values for x then you get different results:

In [17]: solve(eq.subs(x, 0), [a, b])
Out[17]: [(0, b)]

In [18]: solve(eq, [a, b])
Out[18]: {a: 0, b: 0}

Why are these even different types? One is a list of tuples and the
other is a dict...

The return from solve(eq, [a, b]) should be valid for different values
of the parameter x i.e. ideally these should be equivalent:

solve(eq.subs(x, value), [a, b])
solve(eq, [a, b]).subs(x, xvalue)

(Although it isn't actually possible to use subs directly on the
returned output like that.)

> Since this behavior is described in the docstring and has been present for many years it feels like a regression of this feature to me (though to others it feels like a bug fix). If you use this feature and/or have comments, please let us know your thoughts on this:

I don't think this is a feature:

1. I doubt that anyone even realises that this is intentional behaviour.
2. The docstring of solve is impossible to understand.
3. This feature hardly works because there isn't any way to make solve
consistently behave like this.

> 1) should it be kept?

The output should just be changed to be consistent with the
mathematical definition of computing a parametrised solution set.

> 2) if removed, should there be deprecation?

It isn't possible to deprecate this. The output just needs to be made
correct. What is possible is adding a separate function that does
undetermined coefficients solving more directly.

--
Oscar

[gu...@uwosh.edu]

+1 for separating out undetermined coefficients.

I am sure this will break somebody's code. However, I agree that the behavior of solve is so inconsistent that I do not think there is a reasonable way of deprecating this. I generally only use solve in an interactive environment and as little as possible because of how mysterious it is about what choices it is making. I advise my students to avoid using it.

Jonathan

[Oscar Benjamin]

On Fri, 1 Jul 2022 at 13:01, gu...@uwosh.edu <gu...@uwosh.edu> wrote:
>
> +1 for separating out undetermined coefficients.
> I am sure this will break somebody's code. However, I agree that the behavior of solve is so inconsistent that I do not think there is a reasonable way of deprecating this. I generally only use solve in an interactive environment and as little as possible because of how mysterious it is about what choices it is making. I advise my students to avoid using it.

To be clear it is possible to make solve much more consistent both in
terms of type of return values and the solutions it computes if you
always pass a list of solutions (even for a single equation) and also
pass the dict=True flag. Then the output is guaranteed to a be a list
of dicts which is the most useful form because you can use it for
substitution:

In [21]: eq = x**2 - 1

In [22]: solve([eq], [x], dict=True)
Out[22]: [{x: -1}, {x: 1}]

In [23]: [s1, s2] = solve([eq], [x], dict=True)

In [24]: eq.subs(s1)
Out[24]: 0

Internally in solve there are two different codepaths:

https://github.com/sympy/sympy/blob/8a5296712e30cc5fc29b6a519797be737da51b38/sympy/solvers/solvers.py#L1118-L1121

The weirdness as referred to in this thread generally lies down the
bare_f path. Passing a list of equations i.e. [eq] takes the other
path.

For downstream library usage I would recommend to always give a list
of equations and always use dict=True.

--
Oscar

[Chris Smith]

A deprecation is possible. Here is a potential docstring addition:

** Automatic detection of undetermined coefficients**


When solving a single equation that is not passed in a list and
is passed with more than 1 variable for which to solve a check is made
to see whether the equation and variables represent a *linear*
system of equations that can determine coefficients on independent
symbols not passed, This behavior will likely be disabled in the
future; to disable it now, use the keyword `coeffs=False`.

    >>> solve(Eq(a*x - b, 2*x + 3), (a, b))
    {a: 2, b: -3}
    >>> solve(Eq(a*x - b, 2*x + 3), (a, b), coeffs=False)
    [{b: a*x - 2*x - 3}]

At the end of the deprecation period, the first solution would no longer be returned and the `coeffs` keyword would enter a deprecation cycle so that those using it would have time to remove it. Finally, the keyword would be removed.

The user will be directed to use `solve_undetermined_coeffs` and that function will be upgraded to handle linear coefficients on potentially non-polynomial generators on more than one symbol, e.g.

    solve_undetermined_coeffs(Eq(a*sin(x) + b*x*y, 3*sin(x) + 4*x*y), (a,b))

    solve_undetermined_coeffs(Eq(a*x**2 + b*x*y + c*y**2 + d, 3*x**2+4), (a,b,c,d))

/c