Fix calculation of rank of matrix.

[rpisarev]

Fix calculation of rank of a matrix with complicated elements.

This PR able to close #9480 and #11238.

I did some performance test for 3 implementations of _iszero()
So, 1st function (from current master):

def _iszero(x):
    """Returns True if x is zero."""
    return x.is_zero

We can see next:

In [23]: def check_rank():
    time_all = 0
    for k in xrange(1000):
        A = Matrix([[7, 6, 5, 4, 3, 2],
                        [9, 7, 8, 9, 4, 3],
                        [7, 4, 9, 7, 0, 0],
                        [5, 3, 6, 1, 0, 0],
                        [0, 0, 5, 6, 0, 0],
                        [0, 0, 6, 8, 0, sqrt(k)]])
        begin_t = time.time()
        b = A.rank()
        time_all += time.time() - begin_t
    return time_all
   ....: 

In [24]: check_rank()
Out[24]: 5.189665794372559

In [25]: check_rank()
Out[25]: 5.310587406158447

In [26]: check_rank()
Out[26]: 5.2525904178619385

In [27]: check_rank()
Out[27]: 5.224452018737793

In [28]: check_rank()
Out[28]: 5.226367950439453

In [29]: %paste
from sympy import Point, tan, sqrt, pi, Matrix, simplify
x = 8*tan(13*pi/45)/(tan(13*pi/45) + sqrt(3))
y = (-8*sqrt(3)*tan(13*pi/45)**2 + 24*tan(13*pi/45))/(-3 + tan(13*pi/45)**2)
p1 = Point(0, 0)
p2 = Point(1, -sqrt(3))
p0 = Point(x,y)
m1 = Matrix([p1 - simplify(p0), p2 - simplify(p0)])
m2 = Matrix([p1 - p0, p2 - p0])
m3 = Matrix([simplify(p1 - p0), simplify(p2 - p0)])

print(m1.rank(simplify=True))
print(m2.rank(simplify=True))
print(m3.rank(simplify=True))

## -- End pasted text --
2
2
2

If we can use function from #10650,

def _iszero(x):
      """Returns True if x is zero."""
    if hasattr(x, "free_symbols") and not x.free_symbols:
        # if we do not have free symbols, we will call _simplify()
        # without risks to lose some confines for domain of an expression,
        # because it will be some number
        x = _simplify(x)
    return x.is_zero

we will get

In [1]: import time

In [2]: from sympy import *

In [3]: def check_rank():
    time_all = 0
    for k in xrange(1000):
        A = Matrix([[7, 6, 5, 4, 3, 2],
                        [9, 7, 8, 9, 4, 3],
                        [7, 4, 9, 7, 0, 0],
                        [5, 3, 6, 1, 0, 0],
                        [0, 0, 5, 6, 0, 0],
                        [0, 0, 6, 8, 0, sqrt(k)]])
        begin_t = time.time()
        b = A.rank()
        time_all += time.time() - begin_t
    return time_all
   ...: 

In [4]: check_rank()
Out[4]: 23.70626139640808

In [5]: check_rank()
Out[5]: 23.882262468338013

In [6]: check_rank()
Out[6]: 23.626940965652466

In [7]: %paste
from sympy import Point, tan, sqrt, pi, Matrix, simplify
x = 8*tan(13*pi/45)/(tan(13*pi/45) + sqrt(3))
y = (-8*sqrt(3)*tan(13*pi/45)**2 + 24*tan(13*pi/45))/(-3 + tan(13*pi/45)**2)
p1 = Point(0, 0)
p2 = Point(1, -sqrt(3))
p0 = Point(x,y)
m1 = Matrix([p1 - simplify(p0), p2 - simplify(p0)])
m2 = Matrix([p1 - p0, p2 - p0])
m3 = Matrix([simplify(p1 - p0), simplify(p2 - p0)])

print(m1.rank(simplify=True))
print(m2.rank(simplify=True))
print(m3.rank(simplify=True))

## -- End pasted text --
1
2
2

So, the second implementation is slower as first and as wrong, as first.

Third function (this PR):

def _iszero(x):
    """Returns True if x is zero."""
    if hasattr(x, "free_symbols") and not x.free_symbols:
        return x.equals(0)
    else:
        return x.is_zero

And we have output:

In [1]: from sympy import *

In [2]: import time

In [3]: def check_rank():
    time_all = 0
    for k in xrange(1000):
        A = Matrix([[7, 6, 5, 4, 3, 2],
                        [9, 7, 8, 9, 4, 3],
                        [7, 4, 9, 7, 0, 0],
                        [5, 3, 6, 1, 0, 0],
                        [0, 0, 5, 6, 0, 0],
                        [0, 0, 6, 8, 0, sqrt(k)]])
        begin_t = time.time()
        b = A.rank()
        time_all += time.time() - begin_t
    return time_all
   ...: 

In [4]: check_rank()
Out[4]: 25.26923680305481

In [5]: check_rank()
Out[5]: 25.200393676757812

In [6]: check_rank()
Out[6]: 25.13577675819397

In [7]: %paste
from sympy import Point, tan, sqrt, pi, Matrix, simplify
x = 8*tan(13*pi/45)/(tan(13*pi/45) + sqrt(3))
y = (-8*sqrt(3)*tan(13*pi/45)**2 + 24*tan(13*pi/45))/(-3 + tan(13*pi/45)**2)
p1 = Point(0, 0)
p2 = Point(1, -sqrt(3))
p0 = Point(x,y)
m1 = Matrix([p1 - simplify(p0), p2 - simplify(p0)])
m2 = Matrix([p1 - p0, p2 - p0])
m3 = Matrix([simplify(p1 - p0), simplify(p2 - p0)])

print(m1.rank(simplify=True))
print(m2.rank(simplify=True))
print(m3.rank(simplify=True))

## -- End pasted text --
1
1
1

Third implementation is slower as first too, and even as second. But, it is absolutly right.

[asmeurer]

When does x not have free_symbols?

[asmeurer]

At any rate, you can use getattr(x, 'free_symbols', set()).

[rpisarev]

Actually, I think, if x does not have free_symbols, value of hasattr(x, "free_symbols") and not x.free_symbols will False and we use else-branch. But, yes, I can use getattr

[asmeurer]

This should also fix some issues with rref.

[asmeurer]

I think using equals in _iszero is a good fix for this issue (at least until we write matrices with real domains). Raw simplify is too slow, as you saw (it's actually worse, because unlike your benchmark, it's really easy for simplify to just hang forever). is_zero has simple errors. equals should be a good middle ground, that misses most errors, but is performant as well.

[siefkenj]

PR #11554 resolves this issue, but doesn't include any benchmarks. Could this PR be reworked to include benchmarks?

[moorepants]

@ rpisarev Do you want to rework this as benchmarks? If not, we will close the PR.

[rpisarev]

Hi. I don't know how I should add benchmark

[siefkenj]

@rpisarev I think benchmarks are in the separate repository https://github.com/sympy/sympy_benchmarks Can you send a PR there?

[gxyd]

@siefkenj did PR sympy/sympy_benchmarks#34 made all the necessary benchmarks? Should this PR be closed or kept open?

[siefkenj]

The benchmarks added in sympy/sympy_benchmarks#34 are slightly different, but they should be sufficient.

[gxyd]

So I think this PR can be closed then. The necessary code got merged in PR #11554 and the necessary benchmarks were added in sympy/sympy_benchmarks#34 . Thanks @siefkenj for the explanation.