Make all matrix properties accessible as @property

[sylee957]

The properties in MatrixProperties class in common.py are not accessed in consistent manner right now.
Though it may have been a very old issue.

Properties	Functions
is_hermitian	is_symmetric
is_Identity	is_anti_symmetric
is_square	is_diagonal
is_lower_hessenberg
is_upper_hessenberg
is_lower
is_upper

And also,

Properties	Functions
free_symbols	atoms
	values
	is_symbolic

These may not have been consistent,
though But has These may be leaved as-is like before.

One of the problem it had introduced, it that some optional arguments like simplifywould not work once it is set to property,

from sympy import *
a = Matrix.eye(10)
a.is_hermitian(simplify=False)

unlike the code advertises like it should be

sympy/sympy/matrices/common.py

Lines 1155 to 1156 in cb8a524

	@property
	def is_hermitian(self, simplify=True):

because a.is_hermitian as a whole becomes a boolean object after called if set as property, and it would not be callable.
I think there would not be a possible way to overload such thing in python nature, as something like a.is_hermitian to be used as property while also as a.is_hermitian() or a.is_hermitian(simplify=False),
so literally it would be no different than making the function hardcoded with simplify.

It would also mean that it will be tricky to make compatible with legacies, if we make any changes.
Because it will simply not be usable as functions once we add @property decorator with it.

If it should be made changes with, I think there should be some compatibility layer added with some global variables, to make an option to choose whether to use Old_MatrixProperty class or New_MatrixProperty class,
so that it would be possible to throw some deprecation warning that you should change some global variable, and then it may be possible to use a.is_symmetric instead of a.is_symmetric() or such.

[asmeurer]

Generally properties should be cheap to compute, and if it isn't cheap it should be a function. That doesn't seem to be the case here though. However, the fact that there are arguments to the function is another reason it has to be a function.

I agree there's no way to do a real deprecation here if we wanted to change this.

Also, atoms comes from Basic, where it is a function (it takes an optional parameter).

[Teut2711]

I would like to contribute.

[asmeurer]

This isn't a good issue to contribute on. We first need to decide if it is worth making this change. I am thinking it might not, as it would be a backwards incompatible change and we cannot do a deprecation.

[sylee957]

Generally properties should be cheap to compute,

I would personally have some question about, for which criteria should a computation to be considered cheap or expensive?

For example, almost every matrix properties would have at least O(n^2) time complexity, for obvious reason that there should be elementwise investigation to search for a certain property.
Though some properties like is_square can obviously be considered cheap, because it can take advantage of comparing some 'metadata' instead of counting the whole 'data'.

However, to generalize that anything that needs elementwise investigation on 'data' should be considered expensive in SymPy, I think a lot of SymPy's assumption systems on Basic objects are already using recursive sequence to test out assumptions, because they may have to investigate properties of every tree leafs before making a conclusion, and such computation may also have reached the time complexity comparable to matrix computations, though many of they may already have been defined as @property.

[sylee957]

An idea I have is that is_symmetric, is_hermitian should use Eq instead of ==, even if it can slow down in some aspect. And in this way, it can use 3-valued logic to return uncertainty value.

I would make everything wrapped under @property and throw out simplify keyword

I would also split is_diagonalizable into real and complex parts, since there should be no keyword for @property.

I'd try to finish this discussion before next release, since it introduces a lot of mess for working with matrix methods.

Some users may face some abrupt changes, but I'd expect that abruptly changing to using @property can be better since it can the users who were using it as a callable will notice the exceptions like bool is not a callable, but detaching @property will silently introduce some wrong results.

I would expect some other questions like how should things be simplified that I have not answered, but I will tell you the reasons. Any questions?

[tom-pytel]

Instead of splitting is_diagonalizable into two why not add the property is_real which would basically be all(e.is_real for e in self) so that is_diagonalizable would always function as is_diagonalizable(reals_only=False) and if you wanted to restrict to reals then you would check for M.is_real and M.is_diagonalizable?

[sylee957]

Instead of splitting is_diagonalizable into two why not add the property is_real which would basically be all(e.is_real for e in self) so that is_diagonalizable would always function as is_diagonalizable(reals_only=False) and if you wanted to restrict to reals then you would check for M.is_real and M.is_diagonalizable?

Unfortunately, it's a different matter since is_diagonalizable have to check the realness of the eigenvalues, which can not be judged by looking at the matrix elements if the matrix is not hermitian.

[tom-pytel]

Good point, mixed it up with the quick out on symmetric real matrices. How about M.has_real_eigenvals and m.is_diagonalizable, that property could be useful in other contexts?

[oscarbenjamin]

I would prefer something like is_real to be is_real_matrix. At the moment M.is_real means literally that M is a real number which can never be true of a matrix.

We already have bugs from M.is_zero being abused in this way:

In [19]: Matrix([[0, 0], [0, 0]]).is_zero                                                                                                      
Out[19]: True

In [20]: Integral(Matrix([[0, 0], [0, 0]]), x)                                                                                                 
Out[20]: 
⌠          
⎮ ⎡0  0⎤   
⎮ ⎢    ⎥ dx
⎮ ⎣0  0⎦   
⌡          

In [21]: Integral(Matrix([[0, 0], [0, 0]]), x).doit()                                                                                          
Out[21]: 0

[sylee957]

Regarding the diagonalizability, I'd sum up to following 4 cases

1. Orthogonally-diagonalizable with real D matrix
   There is also a proof that every orthogonally-diagonalizable matrix is symmetric.
2. Unitarily-diagonalizable with real D matrix
   There is also a proof that every unitarily-diagonalizable matrix is hermitian.

And some other cases are:

3. Diagonalizable with real D matrix
   
4. Diagonalizable with complex D matrix
   
5. Diagonalizable with real U matrix
   
6. Diagonalizable with complex U matrix
   

which doesn't necessarily need U to be orthogonal or unitary matrix.
I'm not sure that these things can be easily determined before computing eigenvalues though.

So whatever, I had confused about reals_only and orthogonally-diagonalizable before, and reals_only should also be True for hermitian matrices for some sense,
but as @Pristine-Cat 's idea, I'd consider generalizing it as a new property to test if all eigenvalues are real.

[tom-pytel]

On another note, and please disregard if you are sick of hearing of caching, but if you guys don't like the idea of custom matrix caching - @sylee957 because of memory usage for memoizing and @oscarbenjamin for not wanting to add more crud to an already messy matrix interface - then consider adding one or both of the following to the immutable matrix classes:

eigenvals = cacheit(DenseMatrix.eigenvals)
eigenvects = cacheit(DenseMatrix.eigenvects)

This will enable caching of these but only for the immutable classes (which I understand is the direction the matrix module is heading in anyway). For @oscarbenjamin it uses the current SymPy caching mechanism and for @sylee957 it is a fixed size cache so memory usage will not get out of control. And it has a decent impact on a single call to diagonalize due to its double call to eigenvects via is_diagonalizable, regardless of the fact that is_diagonalizable is already cached in this way since the first call to is_diagonalizable has to be evaluated anyway:

M = 
⎡              45   37⋅ⅈ   ⎤
⎢   -3/4       ── - ────   ⎥
⎢              32    16    ⎥
⎢                          ⎥
⎢49⋅z   149    177   1369⋅ⅈ⎥
⎢──── - ───  - ─── - ──────⎥
⎣ 32     64    128    128  ⎦

M.diagonalize()

Without caching this call takes 1.34s, with caching it is 0.81s. Of the two the eigenvects caching is the important one but I included eigenvals since it is so commonly used and relatively small to cache and would certainly help with a M.has_real_eigenvals and m.is_diagonalizable, but is not so critical so if you only want to cache one it would be the eigenvects.

EDIT: I am referring to the behavior of is_diagonalizable before the addition of _is_diagonalizable_with_eigen but I think none of us likes that function.

EDIT2: I mention caching in a theoretical discussion because if all these things become properties and the old functions are deprecated somehow then it makes sense to cache the expensive ones.

[oscarbenjamin]

We should separate out functions like eigenvals. They should be separate functions that sympify their args and use cacheit.

[tom-pytel]

You mean standalone function as opposed to matrix method. Would cacheit work if calling eigenvals(MutableMatrix, ...? And what other functions would you want to separate out, or maybe most/all of them and make the Matrix classes simple interfaces? (By all I mean the substantive ones like in MatrixReductions, MatrixEigen, etc... and not shaping or other matrix utilities.)

[oscarbenjamin]

I (personally) think that all substantial code should be split out as separate functions. There is no need to deprecate the methods but they should be simple wrappers and the docs should suggest the functions:

class Matrix:
    def eigenvals(self, *args, **kwargs):
        '''Compute the eigenvalues of a matrix. See eigenvals function for details'''
        from sympy.matrices.eigen import eigenvals
        return eigenvals(self, *args, **kwargs)

@cacheit
def eigenvals(M):
    M = _sympify(M)
    # Actual code here

Then internally all calculations are done with immutable matrices so anything useful can be cached.

Then we can have a whole module full of different functions for computing eigenvalues/vectors which provides a natural place to document the different methods and when they might be useful.

[tom-pytel]

I understand what you mean, I was just curious if cacheit would work with mutable matrix arguments, I understood we would be returning immutable, but had a look through the lru_cache code and it would handle it ... obviously ... as it does mutable lists or dicts.

On the other hand a design question, is it wise to always force operations to return immutable matrices? Would there be cases where someone would be working with large matrices and would want to change some elements of some of these returned matrices in place - if so they would have to create a new mutable matrix from the returned immutable even if the original argument was mutable.

[oscarbenjamin]

We should soft-deprecate mutable matrices and encourage users to use immutable matrices. Having Matrix as an alias for MutableDenseMatrix was a poor decision that is hard to recover from: most of our users are using mutable matrices without even wanting to.

By "soft-deprecate" I mean that we should make it straight-forward to do everything with immutable matrices and emphasise ways of doing that in the docs. Mutable matrices can be mentioned as a side-note at the end with text that clearly discourages using them.

In terms of the implementation the mutable matrices code should be completely separate from the everything else. There should be no need to worry about the possibility of having mutable matrices in any of the core algorithms etc.

[oscarbenjamin]

if so they would have to create a new mutable matrix from the returned immutable even if the original argument was mutable.

The alternative would be creating a new mutable matrix at the point of return every time. I think that's fine though as long as it's clear how to avoid using mutable matrices in the first place.

[tom-pytel]

So have the new mutable matrix implement wrappers for these new standalone functions which would cast the returned immutable to a mutable?

[tom-pytel]

Well it might not be so bad since many of these functions don't do so well on very large matrices anyway...

[tom-pytel]

The good thing about what you are suggesting is that it can be done incrementally, and then when the matrix classes reach a point where they are essentially wrappers they can be reorganized.

[oscarbenjamin]

The first thing to fix is we need a recommended way to create an immutable matrix and then we should recommend that in the docs:
https://docs.sympy.org/1.5.1/tutorial/matrices.html
Right now every new user is being encouraged to use mutable matrices whether they want to or not.

Unfortunately the name Matrix is already taken and I don't think we can change that for backwards compatibility. I did experiment with making Matrix into a new class that would warn if you try to mutate it (#16790) but I don't think that approach can work.

So we need something like ImMatrix or perhaps just to recommend ImmutableMatrix everywhere in the docs.

[tom-pytel]

You read my mind, I was thinking IMatrix...

[tom-pytel]

Add some deprecation warnings further down the line when instantiating Matrix but not MutableMatrix, MutableDenseMatrix or MutableSparseMatrix forcing people to use those classes if they want a mutable matrix and eventually even later we can reclaim Matrix for immutable?

[oscarbenjamin]

I suppose IMatrix is as close as we can get to Matrix. Then all references to Matrix can be changed either to IMatrix or MutableMatrix:

$ git grep ' Matrix(' | wc -l
    3098

I'm not sure if we'll ever be able to deprecate this fully. It's been too prominent in the docs. Anyone who has written their code in the currently recommended way would have to change it. Look how much of SymPy's own codebase depends on matrices being mutable (#16790) and project that across the userbase.

I think this will just be a long-term soft-deprecation like numpy.matrix which is still used in SymPy in places even though it now gives PendingDeprecationWarning:

sympy/pytest.ini

Lines 7 to 11 in ee2c266

	# np.matrix is deprecated but still used in sympy.physics.quantum:
	# https://github.com/sympy/sympy/issues/15563
	# This prevents ~800 warnings showing up running the tests under pytest.
	filterwarnings =
	ignore:.*the matrix subclass is not the recommended way.*:PendingDeprecationWarning

[tom-pytel]

$ git grep "\bMatrix\s*("

Actually most of those hits are tests, the matrix module itself and .rst files. When those are excluded the count is 767. Not now but eventually if you want to remove mutable matrices from wherever they are not needed in the SymPy codebase you would only need to change all Matrix to IMatrix, see what breaks and then change those to MutableMatrix.

[tom-pytel]

@sylee957

I'm not sure that these things can be easily determined before computing eigenvalues though.

So whatever, I had confused about reals_only and orthogonally-diagonalizable before, and reals_only should also be True for hermitian matrices for some sense,
but as @Pristine-Cat 's idea, I'd consider generalizing it as a new property to test if all eigenvalues are real.

If I understand the reals_only flag from the is_diagonalizable code it means they want only real eigenvals but why is this flag there to begin with? I don't see it used in the SymPy codebase. And why a flag and not a post-diagonalize reality check on D and/or P via something like is_real_matrix like @oscarbenjamin suggested since is_diagonalizable essentially does all the work of diagonalize anyway apart from the normalization?

In fact you could even flip the hierarchy around and have is_diagonalizable call diagonalize(normalize=False, ...) to determine diagonalizability.

[tom-pytel]

@oscarbenjamin - Just sent up a test of what I understand you to mean - #18253