Bug with `solve` in `singularities`

[gituliar]

When running the command

fuchsia reduce -e ep -m pap_01_norm.m -t pap_01_norm_t.m -f m pap_01.m

I get the exception

Traceback (most recent call last):
  File "/home/gituliar/opt/sage/sage-7.0/local/lib/python/runpy.py", line 162, in _run_module_as_main
    "__main__", fname, loader, pkg_name)
  File "/home/gituliar/opt/sage/sage-7.0/local/lib/python/runpy.py", line 72, in _run_code
    exec code in run_globals
  File "/home/gituliar/src/fuchsia/fuchsia.py", line 1418, in <module>
    main()
  File "/home/gituliar/src/fuchsia/fuchsia.py", line 1384, in main
    M, T = canonical_form(m, x, epsilon)
  File "/home/gituliar/src/fuchsia/fuchsia.py", line 475, in canonical_form
    m_f, t_f = fuchsify_by_blocks(m_n, b, x, eps)
  File "/home/gituliar/src/fuchsia/fuchsia.py", line 582, in fuchsify_by_blocks
    pts = singularities(bj, x)
  File "/home/gituliar/src/fuchsia/fuchsia.py", line 194, in singularities
    val = sol[x]
KeyError: x

It looks like the reason is that solve is not returning an original equatiion instead of a propper solution. @magv can you take a look at it?

[gituliar]

Here is the initial equation:

-(ep*x^6*(57657600*I*sqrt(1671) - 1852804800) + ep*x^5*(-348957840*I*sqrt(1671) + 20153463120) + ep*x^4*(1008135240*I*sqrt(1671) - 66627047760) + ep*x^3*(-1474351802*I*sqrt(1671) + 77911868874) + ep*x^2*(616853251*I*sqrt(1671) + 22553597843) + ep*x*(970876844*I*sqrt(1671) - 124333776168) + ep*(-1073990918*I*sqrt(1671) + 78010544766))/(902340000*x^7 + x^6*(15039000*I*sqrt(1671) - 17941527000) + x^5*(-277218900*I*sqrt(1671) + 117306706500) + x^4*(1553144370*I*sqrt(1671) - 376231314690) + x^3*(-4018462575*I*sqrt(1671) + 672623577825) + x^2*(5383622230*I*sqrt(1671) - 683725091910) + x*(-3589165965*I*sqrt(1671) + 367738649355) + 924686840*I*sqrt(1671) - 80447755080)

and returned solution:

[{0: 180468000*x^7 + x^6*(3007800*I*sqrt(1671) - 3588305400) + x^5*(-55443780*I*sqrt(1671) + 23461341300) + x^4*(310628874*I*sqrt(1671) - 75246262938) + x^3*(-803692515*I*sqrt(1671) + 134524715565) + x^2*(1076724446*I*sqrt(1671) - 136745018382) + x*(-717833193*I*sqrt(1671) + 73547729871) + 184937368*I*sqrt(1671) - 16089551016}]

You can check explicitely this result in Sage.

[magv]

Looks like a Maxima bug. It can easily factor that expression, but can't solve it for some reason? Ugh.

[magv]

I've committed a workaround for this problem in 847b1177f02764c40e6254713589c94973ad5b49. It's a hack of course, so we'll need to see how it holds up.

This does not however make pap_01.mtx reducible. One problem is that fuchsify_by_blocks enters an infinite loop around block at 48:66.

A separate problem is at cell 68:0 of the matrix: if you'll try to calculate it's residue at x = -1/60*I*sqrt(1671) + 29/20, you'll get an expression with coefficients around 10^400, and any sort of operation on the resulting expression (like, solving a linear equation with it) will hang.

I think Maxima tries to factor numerical coefficients in some of it's operations, and factoring numbers around 10^400 would make any CAS hang.

---

In any case, close this issue if KeyError went away.

[gituliar]

I changed fuchsify_by_blocks function in 511a38283c0216760b5a595e6e3f1b8b77e967ad, so that now it starts from the simplest sectors and proceed to the more complicated ones; it looks to work much faster. I managed to go through pap_01.m example with it (using Maple), however there is a problem with complex points, so the resulting matrix is not Fuchsian in complex points. Look at test_block_fuchsify_6 and test_fuchsify_by_blocks_07, they are related to this problem.

[magv]

fuchsify_by_blocks function [...] looks to work much faster.

The reason it's faster now (on par with 'block_fuchsify') is because previously 'm = transform(mm, x, t).simplify_rational()' was performed for every non-diagonal block, even when no reduction was performed in that block (which is the majority of cases for sparse matrices); now it is performed only after every reduction.

there is a problem with complex points, so the resulting matrix is not Fuchsian in complex points. Look at test_block_fuchsify_6 and test_fuchsify_by_blocks_07, they are related to this problem.

For Maple too?

[magv]

Also, what does Maple say about residue of cell 68:0 at -1/60*I*sqrt(1671) + 29/20?

[magv]

For the record, this problem wasn't really solved. I have examples where solve still fails in the same way. I don't think there's a quick solution here.

One thing we could do is try maintaining matrix cells in a different normal form: currently we use the fully-expanded polynom ratio form (the one produced by simplify_rational), which forces us to constantly expand (to maintain the form) and factor (to find singularities) expressions; it we could use, for example, fully factored x-polynom ratio form, then finding singularities would be trivial, and the problem with solve would mostly go away: if an expression was factored once (to compute the normal form), all the further solve calls would instantly succeed.

PS. I still want to know what Maple says about residue of pap_01[68:0] at x = -1/60*I*sqrt(1671) + 29/20.

[gituliar]

I still want to know what Maple says about residue of

The residue is 3/112*(13151187756684*ep^3 + 184781726331*(18*sqrt(1671)*ep^3 - 27*sqrt(1671)*ep^2 + 13*sqrt(1671)*ep - 2*sqrt(1671))*I - 19726781635026*ep^2 + 9498080046494*ep - 1461243084076)/(1243372988*sqrt(1671)*I*ep^2 + 9747810249*ep^2).

One thing we could do is try maintaining matrix cells in a different normal form: <...>

I am thinking about this idea for some time aleady and it looks promising. What would be extremely nice to support is expressions of the form (x-x1)^-2 (x-x2)^-3, where x1 and x2 are symbols. In fact, what Costas sent us is this case, but with x1 and x2 substituted by numbers. That is a main feature for the second release (if that will happen).

For the moment it looks easier to switch to Fermat or other free CAS.

[magv]

The residue is |3/112_(13151187756684ep^3 + 184781726331(18_sqrt(1671)_ep^3 - 27_sqrt(1671)_ep^2 + 13_sqrt(1671)*ep

• 2_sqrt(1671))_I - 19726781635026_ep^2 + 9498080046494_ep - 1461243084076)/(1243372988_sqrt(1671)_I_ep^2 + 9747810249_ep^2)|.

Yeah, much better than what Maxima calculates. It seems like Maxima's support for complex numbers is rather poor. For example:

sage: ((x^2 + 2*I*x + I*I)/(x + I)).simplify_rational()
(x^2 + 2*I*x - 1)/(x + I)

The correct answer is, of course, 'x + I'.

This sort of simplification failure propagates through expressions, which is the cause of a number of hangs we've seen.

What would be extremely nice to support is expressions of the form |(x-x1)^-2 (x-x2)^-3|, where |x1| and |x2| are symbols.

Strictly speaking, calculating eigenspaces of symbolic matrices is a logical error. Consider this matrix:

[a 0]
[0 b]

Does it have two eigenvalues of multiplicity 1, or a single eigenvalue of multiplicity 2? Depends on the values of 'a' and 'b'. It's impossible to calculate symbolically, in other words.

Our usage of 'epsilon' is sort of acceptable, because it's infinitesimal, so we know that it's never equal to any fixed value. We can't make that assumption with arbitrary named constants, so...

For the moment it looks easier to switch to Fermat or other free CAS.

Fermat, as far as I can see, does not support complex numbers:

$ ./fer64
>Sqrt(-1)
*** Inappropriate symbol: Sqrt
 can't take square root of negative number.

*** Error occurred at this point:
Sqrt(-1);

It's doc does say "[...] items being computed can be rational numbers, real numbers, complex numbers [...]", but I don't see how that is true.

What other free CAS are there?

[gituliar]

On 05/25/2016 08:10 PM, Vitaly Magerya wrote:

It seems like Maxima's support for complex numbers is rather poor.

Right, the reason I think is that Maxima has no support for complex numbers in GCD, which as you understand is crucial for simplify_rational or similar functions, e.g.,

sage: GCD(x**2+1, x-I)
ValueError: gcd: arguments must be polynomials over the rationals

Strictly speaking, calculating eigenspaces of symbolic matrices is a logical error.

[...]

Our usage of 'epsilon' is sort of acceptable, because it's infinitesimal, so we know that it's never equal to any fixed value. We can't make that assumption with arbitrary named constants, so...

Remember, eigenvalues in our case have the form "n + m eps". As you saw, pap_01.m fits into this rule, i.e., eigenvalues have no relation to substituted numbers. Of course, eigenvectors will depend on symbolic parameters, but that is OK.

Fermat, as far as I can see, does not support complex numbers:

I hope it does, otherwise we can workaround this issue by introducing additional parameter I, and then simplify its higher powers.

What can you say about "solve" functionality, can we use Fermat, for example, in factorization step, instead of Maple as we do now for complicated matrices?

What other free CAS are there?

GiNaC, Magma, Singular.

[magv]

sage: GCD(x**2+1, x-I) ValueError: gcd: arguments must be polynomials over the rationals

Yeah, it's bad in all kinds of ways.

Remember, eigenvalues in our case have the form "n + m eps". As you saw, pap_01.m fits into this rule, i.e., eigenvalues have no relation to substituted numbers. Of course, eigenvectors will depend on symbolic parameters, but that is OK.

Can't both 'n' and 'm' depend on the symbols?

Let's say we have this matrix:

[(a+eps)/x 0]
[        0 0]

How many normalization step do you need to perform to normalize it? 'a' steps?

I hope it does, otherwise we can workaround this issue by introducing additional parameter I, and then simplify its higher powers.

Won't we have the same problem as with Maxima's complex number support then?

What can you say about "solve" functionality, can we use Fermat, for example, in factorization step, instead of Maple as we do now for complicated matrices?

Fermat doesn't have explicit support for solving linear equation systems, we'd need to convert our equations into a matrix multiplication problem, and then use Fermat's matrix inversion capabilities to solve the system. I guess it's possible (sans the complex number issue).

[gituliar]

Can't both 'n' and 'm' depend on the symbols?

Though I have no proof, but it can not.

Let's say we have this matrix:

[(a+eps)/x 0] [ 0 0]

We would rather deal with matrices like

[(m+n*eps)/(x+a) 0] [0 0]

I hope it does, otherwise we can workaround this issue by introducing additional parameter I, and then simplify its higher powers.

Won't we have the same problem as with Maxima's complex number support then?

Right, then we need support for multivariate polynomials, which Fermat does to some extent. I still believe that it has decent support for complex numbers and this trick will not be needed, if we decide to use Fermat of course.

[magv]

Can't both 'n' and 'm' depend on the symbols?

Though I have no proof, but it can not.

Take a look at this matrix:

[1/x/(x-a) 0]
[        0 0]

At x=0 it's residue is:

[-1/a 0]
[   0 0]

So, it's eigenvalues are 0 and -1/a.

[...] Fermat [...]

After looking more into it (it it just me, or is Fermat's documentation really hard to parse?), I think Fermat won't help us. Here's the reason:

>Sqrt(2);

        1

Fermat operates in two modes: rational and modular. From the docs:

In rational arithmetic, this function [Sqrt] returns
the largest integer less than or equal to the square
root. In modular mode, this function is disabled.

So, there is no square root function.

>2^(1/2)

*** Inappropriate symbol: ^
 exponent must be integer, not rational or polynomial.

Nope, no square roots. This means we can't pass an irrational expression into Fermat, and we have plenty of those (roots of polynomials), so I think Fermat is disqualified.