Use gaussian elimination for positive-definite check

[sylee957]

References to other Issues or PRs

Fixes #18955

Brief description of what is fixed or changed

I have changed positive definite check to use division free gaussian elimination, which hopefully gives less bloated expressions, and used berkowitz determinant check for positive semidefinite check which can suffer less from zero testing.

Other comments

Release Notes

• matrices
  • Added is_strongly_diagonally_dominant and is_weakly_diagonally_dominant properties for Matrix.

[sympy-bot]

✅

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:

• matrices
  • Added is_strongly_diagonally_dominant and is_weakly_diagonally_dominant properties for Matrix. (#19205 by @sylee957)

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

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.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

Fixes #18955

#### Brief description of what is fixed or changed

I have changed positive definite check to use division free gaussian elimination, which hopefully gives less bloated expressions, and used berkowitz determinant check for positive semidefinite check which can suffer less from zero testing.

#### Other comments

![Screen Shot 2020-04-28 at 01 02 29](https://user-images.githubusercontent.com/34944973/80393881-112c5980-88ec-11ea-9f50-65cd3613eb6a.png)
![Screen Shot 2020-04-28 at 01 02 47](https://user-images.githubusercontent.com/34944973/80393892-14274a00-88ec-11ea-87d5-d33aaa7164f9.png)


#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
- matrices
  - Added `is_strongly_diagonally_dominant` and `is_weakly_diagonally_dominant` properties for `Matrix`.
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[codecov]

Codecov Report

Merging #19205 into master will decrease coverage by 0.015%.
The diff coverage is 100.000%.

@@              Coverage Diff              @@
##            master    #19205       +/-   ##
=============================================
- Coverage   75.666%   75.651%   -0.016%     
=============================================
  Files          651       651               
  Lines       169336    169383       +47     
  Branches     39970     39981       +11     
=============================================
+ Hits        128131    128140        +9     
- Misses       35605     35637       +32     
- Partials      5600      5606        +6     

[oscarbenjamin]

Are these tests not needed any more?

[sylee957]

I have removed those methods so they cannot be tested

[oscarbenjamin]

If the matrix is not square should we say that is_positive_definite -> False? What about if it is not Hermitian?

[sylee957]

I have left them undefined because I found nobody had asserted that any nonsquare or nonhermitian matrix is not positive definite.

But I have found a generalization that can be extended for nonhermitian matrices in https://www.jstor.org/stable/2317709, which also notes about some properties like the existance of the LR factorization or the positive eigenvalues.

Though False can be asserted for nonsquare matrices, however, I don't think that it loses such triviality because checking whether a matrix is hermitian or square can be done parallel with is_hermitian or is_square.

[oscarbenjamin]

That paper has other useful properties for example that if a matrix positive definite then the diagonal entries are positive (not sure if this only applies for real matrices).

The contrapositive leads to a much cheaper check of positive definiteness: if any diagonal entry is not positive then the matrix is not positive definite.

[oscarbenjamin]

It's probably worth having a shortcut for the case where the answer is probably unknowable. The original issue here came from a symmetric matrix where every entry was a symbol with no assumptions.

I think that if we literally don't know whether any element on the diagonal is positive or not then it's very unlikely that we will be able to get determine if the matrix as a whole is positive definite e.g.:

if all(self[i, i].is_positive is None for i in range(N)):
    return None

Can you think of any examples where we don't know if any of the diagonals is positive or not but we do know that the matrix is or is not positive definite?

[oscarbenjamin]

Actually I guess this would be an example:

In [9]: x = Symbol('x', real=True)                                                                                                             

In [10]: M = Matrix([[sin(x), 2], [2, sin(x)]])                                                                                                

In [11]: M                                                                                                                                     
Out[11]: 
⎡sin(x)    2   ⎤
⎢              ⎥
⎣  2     sin(x)⎦

And in this example we don't know whether anything is positive or negative but we know at least one eigenvalue is negative either way:

In [21]: M = Matrix([[sin(x), 2*sin(x)], [2*sin(x), sin(x)]])                                                                                  

In [22]: M                                                                                                                                     
Out[22]: 
⎡ sin(x)   2⋅sin(x)⎤
⎢                  ⎥
⎣2⋅sin(x)   sin(x) ⎦

In [23]: M.eigenvals()                                                                                                                         
Out[23]: {-sin(x): 1, 3⋅sin(x): 1}

[oscarbenjamin]

Looks good. Merging

[oscarbenjamin]

What exactly is the relationship between row and column diagonal dominance? Are they equivalent?

[sylee957]

They are trivially equivalent for symmetric or hermitian matrices which is sufficient for this test. For other cases I'm not familiar enough to assert.