Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

matrices: Allow inverse call on symbolic matrix #19217

Closed

Conversation

vanillafrosty
Copy link

Fixes #19162
A symbolic matrix A whose elements can be retrieved by indexing into A (A[i], where i is an integer),
is now invertible after A is made explicit with A.as_explicit(). This will return a new subclass of
ImmutableDenseMatrix, where only the inverse method is overwritten.

For example:
A = MatrixSymbol('A', 2, 2)
exA = A.as_explicit()
exA.inv()

Now works.

Release Notes

  • matrices
    • fixed a bug where explicit, symbolic matrices could not be inverted.

@sympy-bot
Copy link

sympy-bot commented Apr 28, 2020

Hi, I am the SymPy bot (v158). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.7.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it.

Click here to see the pull request description that was parsed.

Fixes #19162 
A symbolic matrix A whose elements can be retrieved by indexing into A (A[i], where i is an integer),
is now invertible after A is made explicit with A.as_explicit(). This will return a new subclass of
ImmutableDenseMatrix, where only the inverse method is overwritten.

For example:
A = MatrixSymbol('A', 2, 2)
exA = A.as_explicit()
exA.inv()

Now works. 

#### Release Notes

<!-- BEGIN RELEASE NOTES -->
* matrices
  * fixed a bug where explicit, symbolic matrices could not be inverted. 

<!-- END RELEASE NOTES -->

@oscarbenjamin
Copy link
Contributor

Thanks but I think this misses the actual problem.

I don't think that we want a new class to solve this. We just need to fix the bug that stops inverse working for the previous class.

@sylee957
Copy link
Member

I don't think that it should return unevaluated inverse for explicit matrices, computing the actual inverse should be its own design goal.
If there are some other intermediate computations that are just expensive but can be displayed in some lazy form, it can be a valid option, but I don't think that Inverse can serve for that purpose because it is one of the end goal.

@vanillafrosty
Copy link
Author

@sylee957 I agree, and have made updates to what I believe was the root of the issue. See my latest comment in issue #19162. The actual inverse is now computed.

@oscarbenjamin
Copy link
Contributor

The root of the issue is b.equals(a) raising. It would be better to fix that than catch the exception.

@oscarbenjamin
Copy link
Contributor

As noted in #19162 the MatrixElement constructor needs to check its arguments. I think if the exception was raised from subs then .equals would catch it already.

@vanillafrosty
Copy link
Author

vanillafrosty commented Apr 30, 2020

@oscarbenjamin What's the reasoning behind adding checks in the constructor? I'm not seeing that .subs creates new instances of the instance it gets called on so I don't think we can raise an exception from the constructor? My Python experience (or lack thereof) may be holding me back here.

I can see adding a check for name not being a Number in __new__ in MatrixElement, is that what you mean? For this issue in particular shouldn't we be doing something like overriding _eval_subs in MatrixElement to fail if we pass a Number in for substitution?

@oscarbenjamin
Copy link
Contributor

When you call subs a new object is created with the new substituted arguments so it goes back through __new__:

In [7]: class Thing(Basic): 
   ...:     def __new__(self, arg): 
   ...:         if isinstance(arg, Symbol): 
   ...:             return super().__new__(Thing, arg) 
   ...:         else: 
   ...:             raise ValueError("Bad arguments") 
   ...:                                                                                                                           

In [8]: a = Thing(x)                                                                                                              

In [9]: a                                                                                                                         
Out[9]: Thing(x)

In [10]: a.subs(x, y)                                                                                                             
Out[10]: Thing(y)

In [11]: a.subs(x, 3)                                                                                                             
---------------------------------------------------------------------------
ValueError

I think that for MatrixElement the first arg should be a MatrixSymbol.

@sylee957
Copy link
Member

sylee957 commented May 1, 2020

MatrixElement has a wider scope than MatrixSymbol. For example, it can be used with MatAdd or indexing MatrixBase

A = MatrixSymbol('A', 2, 2)
B = MatrixSymbol('B', 2, 2)
MatrixElement(A+B, 0, 1)

If you need an atomic class only for the MatrixSymbol, I don't think that MatrixElement should be modified,

@vanillafrosty
Copy link
Author

vanillafrosty commented May 2, 2020

@sylee957 are you suggesting we create a class similar to MatrixSymbol with checks in the constructor for the name arg? The issue with modifying the constructor in MatrixElement or MatrixSymbol is that the .inv() will error out with the exception you raise. For simplicity I added an else condition in the MatrixElement constructor:
if isinstance(name, str): name = Symbol(name)
else: raise ValueError('bad element in MatrixElement constructor')

This will raise an exception if I do this:
A = MatrixSymbol("A",2,2)
A.as_explicit().inv()

The behavior is the same if I modify MatrixSymbol's constructor. Of course this means we should be catching the exception somewhere. But is that "somewhere" better than the place I've caught the TypeError in .is_constant?

Catching the TypeError in expr.py in my latest commit does not change the behavior of is_constant or .equals in expr.py. These methods are meant to try a bunch of different expressions and catch exceptions, trying to determine whether the argument passed to it is a constant or whether two things equal each other. The .inv() properly computes the matrix inverse now, with row operations computed.

Let me know what you guys think

@oscarbenjamin
Copy link
Contributor

The problem is that the exception is being raised in the wrong place and the diff here tries to compensate for that by catching it in the wrong place. The right place to raise this is in the call to subs (MatrixElement.__new__). The right place to catch it is around the call to subs. The original issue is fixed with this diff:

diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index 388595f428..9c7e0159d8 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -675,7 +675,7 @@ def check_denominator_zeros(expression):
                 if a is S.NaN:
                     # evaluation may succeed when substitution fails
                     a = expr._random(None, 0, 0, 0, 0)
-            except ZeroDivisionError:
+            except (ZeroDivisionError, TypeError):
                 a = None
             if a is not None and a is not S.NaN:
                 try:
@@ -684,7 +684,7 @@ def check_denominator_zeros(expression):
                     if b is S.NaN:
                         # evaluation may succeed when substitution fails
                         b = expr._random(None, 1, 0, 1, 0)
-                except ZeroDivisionError:
+                except (ZeroDivisionError, TypeError):
                     b = None
                 if b is not None and b is not S.NaN and b.equals(a) is False:
                     return False
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index 56f080da5c..229da0e82b 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -712,6 +712,9 @@ def __new__(cls, name, n, m):
         if isinstance(name, str):
             name = Symbol(name)
         name = _sympify(name)
+        if not isinstance(name, (Symbol, MatrixExpr)):
+            msg = "name should be Symbol/MatrixExpr not %s"
+            raise TypeError(msg % (type(name),))
         obj = Expr.__new__(cls, name, n, m)
         return obj

That gives:

In [1]: import sympy as sp                                                                                                                     

In [2]: X = sp.Matrix(sp.MatrixSymbol('X', 2, 2))                                                                                              

In [3]: X.inv()                                                                                                                                
Out[3]: 
⎡       X₁₁               -X₀₁       ⎤
⎢─────────────────  ─────────────────⎥
⎢X₀₀X₁₁ - X₀₁X₁₀  X₀₀X₁₁ - X₀₁X₁₀⎥
⎢                                    ⎥
⎢      -X₁₀                X₀₀       ⎥
⎢─────────────────  ─────────────────⎥
⎣X₀₀X₁₁ - X₀₁X₁₀  X₀₀X₁₁ - X₀₁X₁₀⎦

However I dislike catching generic TypeError there. We should make a new BadArguments error or something like that in basic.py and catch that instead. It could be a subclass of TypeError. Then other Basic subclasses could raise that from __new__.

@vanillafrosty
Copy link
Author

I may have pushed the latest commit too soon, I see I'm failing some of the tests. Will go check.

@sylee957 please feel free to weigh in if you think MatrixElement shouldn't be modified. I'm not sure what you mean by creating an atomic class for MatrixSymbol.

@vanillafrosty
Copy link
Author

There is still an issue with my latest commit - KroneckerDelta computation seems to be broken. In line 553 in matexpr.py, we're setting identity to S.One. This will not work with the new changes when we construct a MatrixElement on the next line.

I don't know enough about KroneckerDelta to set identity to something that will make the test in line 259 of sympy/matrices/expressions/tests/test_indexing.py work. And I am short on free time. If someone could look at this that would be great, we are so close to fixing the original issue!

@oscarbenjamin
Copy link
Contributor

Does this work?

diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index 56f080da5c..959474cb77 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -548,9 +548,11 @@ def recurse_expr(expr, index_ranges={}):
             elif isinstance(expr, KroneckerDelta):
                 i1, i2 = expr.args
                 if dimensions is not None:
-                    identity = Identity(dimensions[0])
+                    size  = dimensions[0]
                 else:
-                    identity = S.One
+                    range1 = index_ranges[i1]
+                    size = range1[1] - range1[0] + 1
+                identity = Identity(size)
                 return [(MatrixElement(identity, i1, i2), (i1, i2))]
             elif isinstance(expr, MatrixElement):
                 matrix_symbol, i1, i2 = expr.args

I'm not sure I understand that function fully though and it does seem like there could be some other problems there...

@vanillafrosty
Copy link
Author

Yes that works! It looks correct to me, but then again I'm not sure how to write comprehensive tests for that KroneckerDelta computation.

@sylee957
Copy link
Member

sylee957 commented May 4, 2020

Although substituting MatrixSymbol with scalar constants should be disallowed in any other patches, I don't think that it should rely on error catching.

The actual problem is is_constant is forcing zero or one regardless of the types of the original argument.
So I think that the better patch should be connecting them with respective zero or one matrix.

diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index 388595f428..694d2e7372 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -665,13 +665,24 @@ def check_denominator_zeros(expression):
         if expr.is_zero:
             return True

-        # try numerical evaluation to see if we get two different values
+        def zero(wrt):
+            from sympy.matrices.immutable import ImmutableDenseMatrix as Matrix
+            if wrt.is_Matrix:
+                return Matrix.zeros(wrt.rows, wrt.cols)
+            return S.Zero
+
+        def one(wrt):
+            from sympy.matrices.immutable import ImmutableDenseMatrix as Matrix
+            if wrt.is_Matrix:
+                return Matrix.ones(wrt.rows, wrt.cols)
+            return S.One
+
         failing_number = None
         if wrt == free:
             # try 0 (for a) and 1 (for b)
             try:
-                a = expr.subs(list(zip(free, [0]*len(free))),
-                    simultaneous=True)
+                subs_dict = {x: zero(x) for x in free}
+                a = expr.subs(subs_dict, simultaneous=True)
                 if a is S.NaN:
                     # evaluation may succeed when substitution fails
                     a = expr._random(None, 0, 0, 0, 0)
@@ -679,8 +690,8 @@ def check_denominator_zeros(expression):
                 a = None
             if a is not None and a is not S.NaN:
                 try:
-                    b = expr.subs(list(zip(free, [1]*len(free))),
-                        simultaneous=True)
+                    subs_dict = {x: one(x) for x in free}
+                    b = expr.subs(subs_dict, simultaneous=True)
                     if b is S.NaN:
                         # evaluation may succeed when substitution fails
                         b = expr._random(None, 1, 0, 1, 0)

@vanillafrosty
Copy link
Author

@sylee957 that approach works for me. I've kept the new BadArgumentsError for usage going forward, and took out the place we were catching it in is_constant. However, even with your diff above, the new BadArgumentsError is raised in the MatrixElement constructor when .subs tries to replace a symbol with the multiplication of two Dummy instances, in line 963 of basic.py.

I've added a conditional check for MatrixSymbol, and .inv works now. The tests pass as well. However my update to the Dummy multiplication may not be optimal, feel free to take a look.

                if isinstance(old, MatrixSymbol):
                    rv = rv._subs(old, Symbol(str(d*m)))
                else:
                    rv = rv._subs(old, d*m, **kwargs)

@oscarbenjamin
Copy link
Contributor

line 963 of basic.py

I don't know why that code is there in basic.py but I think it's been problematic in other situations so maybe there's a way to get rid of it.

@oscarbenjamin
Copy link
Contributor

The tests have failed because the newline was removed at the end of basic.py

sympy/core/expr.py Outdated Show resolved Hide resolved
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

MatrixSymbol inversion fails.
6 participants