Handling ssqrt and scbrt returned from symbolic_solve

[akirakyle]

I was playing around with the exciting developments landed in #1192 and noticed some peculiarities.

First cosider this example:

using Symbolics, Groebner, Nemo, Latexify, LinearAlgebra

@variables A[1:2, 1:2], λ
char_poly = det(Num.(collect(A-λ*one(A))))
sol = symbolic_solve(char_poly, λ);
simplify(*((λ .- sol)...) - char_poly; expand=true)

which gives me

$$\frac{1}{4} \left( A_{1,1}^{2} + A_{2,2}^{2} \right) - \frac{1}{2} A_{1,1} A_{2,2} + A_{1,2} A_{2,1} - \frac{1}{4} \left( ssqrt\left( A_{1,1}^{2} - 2 A_{1,1} A_{2,2} + 4 A_{1,2} A_{2,1} + A_{2,2}^{2} \right) \right)^{2}$$

As I understand (although would be good to put somewhere in the docs), the ssqrt function was introduced to handle negative arguments without raising errors. Thus I would expect sqrt(x)^2 = x would be a valid rule to have applied by simplify? If I add a rule to rewrite ssqrt to be a power instead I get the expected simplification so the following gives zero as expected.

@variables A[1:2, 1:2], λ
char_poly = det(Num.(collect(A-λ*one(A))))
sol = symbolic_solve(char_poly, λ);
r = @rule Symbolics.ssqrt(~x) => (~x)^(1//2)
sol = simplify.(sol; rewriter=RuleSet([r]))
simplify(*((λ .- sol)...) - char_poly; expand=true)

I then went to try this strategy of checking symbolic_solve with a generic third order characteristic polynomial and it seemed to fail to simplify

@variables A[1:3, 1:3], λ
char_poly = det(Num.(collect(A-λ*one(A))))
sol = symbolic_solve(char_poly, λ);
r1 = @rule Symbolics.ssqrt(~x) => (~x)^(1//2)
r2 = @rule Symbolics.scbrt(~x) => (~x)^(1//3)
sol = simplify.(sol; rewriter=RuleSet([r1, r2]))
simplify(*((λ .- sol)...) - char_poly; expand=true)

The output is far to long to paste here, but it looks like there are floating point coefficients introduced at some point which might be why it fails to simplify to zero?

I was hoping to use the new symbolic_solve to find the eigenvalues of small symbolic matrices so this was my first test to make sure I was understanding how to work with the output of symbolic_solve.

[sumiya11]

Thus I would expect sqrt(x)^2 = x would be a valid rule to have applied by simplify?

I think you are correct. symbolic_solve assumes that ssqrt is the single principal root. We could document it better.

Note that this simplification already happens as postprocessing in symbolic_solve (but it is not a part of simplify)

Symbolics.jl/src/solver/postprocess.jl

Line 91 in 195c95d

# (sqrt(N))^M => N^div(M, 2)*sqrt(N)^(mod(M, 2))

As I understand (although would be good to put somewhere in the docs), the ssqrt function was introduced to handle negative arguments without raising errors

I agree; yes.

The output is far to long to paste here, but it looks like there are floating point coefficients introduced at some point which might be why it fails to simplify to zero?

Hm. If I am not mistaken the option expand=true in simplify may convert numbers to floats.

[akirakyle]

Thanks for the reply! I dug into what is going on a bit more and it does not appear to be expand which causes the floats to appear since inspecting (the really long) output of

@variables A[1:3, 1:3], λ
expr = det(Num.(collect(A-λ*one(A))))
sol = symbolic_solve(expr, λ);
simp = expand(*((λ .- sol)...) - expr)

shows no floating points.
Investigating a bit more it appears they are introduced due to unstable_pow in SymbolicUtils.jl which is illustrated with the example

(2*λ^2)^(1//4)

which outputs $1.1892 \lambda^{\frac{1}{2}}$.

Is there some way to symplify powers of the output of symbolic_solve without them being converted to floats? This seems to be what prevents the above check on a generic characteristic polynomial of order 3 from succeeding and would limit the usefulness for further symbolic computations involving expressions with powers of the corresponding matrix's eigenvalues.

[akirakyle]

If I'm understanding the situation with respect to symbolic representations of purely numeric coefficients which may be irrational, it seems like perhaps landing the PR in JuliaSymbolics/SymbolicUtils.jl#617 might help the situation for symbolic_solve as the numeric coefficents could be wrapped in a Const then Pow wouldn't apply unstable_pow?