Skip to content

fix(matrices): Use QQ[x,y] for Matrix.inv over RR[x,y] #26833

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.

Sign up for GitHub

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

Jump to bottom
Merged
merged 1 commit into from
Jul 19, 2024
Merged

fix(matrices): Use QQ[x,y] for Matrix.inv over RR[x,y] #26833

merged 1 commit into from
Jul 19, 2024

Conversation

oscarbenjamin
Copy link
Collaborator

Use exact domains for matrices of polynomials because linear algebra with matrices of polynomials with float coefficients behaves poorly. Usual techniques for controlling floating point error cannot be used with polynomials and polynomial division fails with floating point coefficients.

This fixes gh-26821 in which both matrix inverse and matrix exponential over QQ[x,y,...] were found to produce expressions that are correct but have unnecessarily large floating point coefficients like 1e200. This also likely fixes other cases in which Matrix.inv would return incorrect results because of explosive rounding errors.

References to other Issues or PRs

Backport of gh-26822 for 1.13 branch

Brief description of what is fixed or changed

Other comments

Release Notes

  • matrices
    • Matrix.inv() now uses exact rational coefficients internally for matrices of polynomials or rational functions if the matrix has symbolic expressions with float coefficients. This fixes an issue first seen in SymPy 1.13.0 where unnecessarily large floats were present in the expressions for the matrix inverse (Matrix exponential badly scaled coefficients for RR[x,y,...]. #26821) and likely also prevents Matrix.inv from returning some inaccurate results for some matrices containing floats.

@oscarbenjamin
fix(matrices): Use QQ[x,y] for Matrix.inv over RR[x,y]
da3a790 
Use exact domains for matrices of polynomials because linear algebra
with matrices of polynomials with float coefficients behaves poorly.
Usual techniques for controlling floating point error cannot be used
with polynomials and polynomial division fails with floating point
coefficients.

This fixes sympygh-26821 in which both matrix inverse and matrix exponential
over QQ[x,y,...] were found to produce expressions that are correct but
have unnecessarily large floating point coefficients like 1e200. This
also likely fixes other cases in which Matrix.inv would return incorrect
results because of explosive rounding errors.
@sympy-bot
Copy link

sympy-bot commented Jul 19, 2024 •
edited
Loading

✅

Hi, I am the SymPy bot. 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
    • Matrix.inv() now uses exact rational coefficients internally for matrices of polynomials or rational functions if the matrix has symbolic expressions with float coefficients. This fixes an issue first seen in SymPy 1.13.0 where unnecessarily large floats were present in the expressions for the matrix inverse (Matrix exponential badly scaled coefficients for RR[x,y,...]. #26821) and likely also prevents Matrix.inv from returning some inaccurate results for some matrices containing floats. (#26833 by @oscarbenjamin)

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

Click here to see the pull request description that was parsed. 
Use exact domains for matrices of polynomials because linear algebra with matrices of polynomials with float coefficients behaves poorly. Usual techniques for controlling floating point error cannot be used with polynomials and polynomial division fails with floating point coefficients.

This fixes gh-26821 in which both matrix inverse and matrix exponential over QQ[x,y,...] were found to produce expressions that are correct but have unnecessarily large floating point coefficients like 1e200. This also likely fixes other cases in which Matrix.inv would return incorrect results because of explosive rounding errors.

<!-- 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. -->

Backport of gh-26822 for 1.13 branch

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


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

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
   * Matrix.inv() now uses exact rational coefficients internally for matrices of polynomials or rational functions if the matrix has symbolic expressions with float coefficients. This fixes an issue first seen in SymPy 1.13.0 where unnecessarily large floats were present in the expressions for the matrix inverse (https://github.com/sympy/sympy/issues/26821) and likely also prevents Matrix.inv from returning some inaccurate results for some matrices containing floats.
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

@oscarbenjamin oscarbenjamin added this to the SymPy 1.13 milestone Jul 19, 2024
@github-actions GitHub Actions
Copy link

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)



Significantly changed benchmark results (master vs previous release)


| Change   | Before [a36a8b23] <sympy-1.12.1^0>   | After [c1b1d2d7]    |   Ratio | Benchmark (Parameter)                                                |
|----------|--------------------------------------|---------------------|---------|----------------------------------------------------------------------|
| -        | 71.5±0.7ms                           | 45.8±0.6ms          |    0.64 | integrate.TimeIntegrationRisch02.time_doit(10)                       |
| -        | 70.0±1ms                             | 43.7±0.3ms          |    0.62 | integrate.TimeIntegrationRisch02.time_doit_risch(10)                 |
| +        | 18.0±0.9μs                           | 29.9±0.1μs          |    1.67 | integrate.TimeIntegrationRisch03.time_doit(1)                        |
| -        | 5.39±0.05ms                          | 2.92±0.01ms         |    0.54 | logic.LogicSuite.time_load_file                                      |
| -        | 71.7±0.4ms                           | 28.3±0.3ms          |    0.4  | polys.TimeGCD_GaussInt.time_op(1, 'dense')                           |
| -        | 25.9±0.3ms                           | 16.9±0.06ms         |    0.65 | polys.TimeGCD_GaussInt.time_op(1, 'expr')                            |
| -        | 73.1±0.2ms                           | 28.6±0.1ms          |    0.39 | polys.TimeGCD_GaussInt.time_op(1, 'sparse')                          |
| -        | 256±2ms                              | 123±0.6ms           |    0.48 | polys.TimeGCD_GaussInt.time_op(2, 'dense')                           |
| -        | 257±2ms                              | 124±0.4ms           |    0.48 | polys.TimeGCD_GaussInt.time_op(2, 'sparse')                          |
| -        | 650±4ms                              | 367±2ms             |    0.56 | polys.TimeGCD_GaussInt.time_op(3, 'dense')                           |
| -        | 655±3ms                              | 367±1ms             |    0.56 | polys.TimeGCD_GaussInt.time_op(3, 'sparse')                          |
| -        | 506±2μs                              | 292±3μs             |    0.58 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(1, 'dense')            |
| -        | 1.77±0.02ms                          | 1.04±0ms            |    0.59 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(2, 'dense')            |
| -        | 5.79±0.03ms                          | 3.07±0.02ms         |    0.53 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(3, 'dense')            |
| -        | 451±3μs                              | 232±0.9μs           |    0.51 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(1, 'dense')               |
| -        | 1.46±0.01ms                          | 671±3μs             |    0.46 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(2, 'dense')               |
| -        | 4.85±0.02ms                          | 1.63±0.01ms         |    0.34 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(3, 'dense')               |
| -        | 379±3μs                              | 207±1μs             |    0.55 | polys.TimeGCD_SparseGCDHighDegree.time_op(1, 'dense')                |
| -        | 2.42±0.05ms                          | 1.23±0.01ms         |    0.51 | polys.TimeGCD_SparseGCDHighDegree.time_op(3, 'dense')                |
| -        | 9.83±0.05ms                          | 4.34±0.03ms         |    0.44 | polys.TimeGCD_SparseGCDHighDegree.time_op(5, 'dense')                |
| -        | 363±3μs                              | 170±0.9μs           |    0.47 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(1, 'dense')            |
| -        | 2.48±0.01ms                          | 903±5μs             |    0.36 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(3, 'dense')            |
| -        | 9.39±0.02ms                          | 2.61±0.01ms         |    0.28 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(5, 'dense')            |
| -        | 1.03±0ms                             | 416±3μs             |    0.4  | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'dense')           |
| -        | 1.77±0.02ms                          | 515±5μs             |    0.29 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'sparse')          |
| -        | 5.92±0.04ms                          | 1.77±0.01ms         |    0.3  | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'dense')           |
| -        | 8.43±0.2ms                           | 1.51±0ms            |    0.18 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'sparse')          |
| -        | 293±0.9μs                            | 65.0±0.8μs          |    0.22 | polys.TimePREM_QuadraticNonMonicGCD.time_op(1, 'sparse')             |
| -        | 3.41±0.04ms                          | 386±2μs             |    0.11 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'dense')              |
| -        | 3.96±0.07ms                          | 282±2μs             |    0.07 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'sparse')             |
| -        | 7.06±0.08ms                          | 1.26±0.01ms         |    0.18 | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'dense')              |
| -        | 8.74±0.04ms                          | 856±6μs             |    0.1  | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'sparse')             |
| -        | 5.09±0.03ms                          | 3.01±0.01ms         |    0.59 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(2, 'sparse') |
| -        | 12.0±0.1ms                           | 6.48±0.04ms         |    0.54 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'dense')  |
| -        | 22.5±0.08ms                          | 9.20±0.07ms         |    0.41 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'sparse') |
| -        | 5.32±0.02ms                          | 882±2μs             |    0.17 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(1, 'sparse')    |
| -        | 12.7±0.03ms                          | 7.01±0.09ms         |    0.55 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(2, 'sparse')    |
| -        | 100±0.7ms                            | 25.6±0.2ms          |    0.26 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'dense')     |
| -        | 168±0.9ms                            | 54.2±0.2ms          |    0.32 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'sparse')    |
| -        | 175±0.8μs                            | 114±1μs             |    0.65 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(1, 'dense')      |
| -        | 368±3μs                              | 218±0.8μs           |    0.59 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(1, 'sparse')     |
| -        | 4.24±0.03ms                          | 826±4μs             |    0.19 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'dense')      |
| -        | 5.19±0.03ms                          | 391±6μs             |    0.08 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'sparse')     |
| -        | 19.9±0.4ms                           | 2.75±0.03ms         |    0.14 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'dense')      |
| -        | 22.6±0.2ms                           | 643±3μs             |    0.03 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'sparse')     |
| -        | 485±7μs                              | 135±2μs             |    0.28 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(1, 'sparse') |
| -        | 4.62±0.02ms                          | 603±1μs             |    0.13 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'dense')  |
| -        | 5.22±0.05ms                          | 138±0.8μs           |    0.03 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'sparse') |
| -        | 13.0±0.08ms                          | 1.29±0.01ms         |    0.1  | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'dense')  |
| -        | 13.9±0.2ms                           | 144±1μs             |    0.01 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'sparse') |
| -        | 134±0.3μs                            | 76.4±0.9μs          |    0.57 | solve.TimeMatrixOperations.time_rref(3, 0)                           |
| -        | 252±0.9μs                            | 88.6±1μs            |    0.35 | solve.TimeMatrixOperations.time_rref(4, 0)                           |
| -        | 24.1±0.2ms                           | 10.3±0.04ms         |    0.43 | solve.TimeSolveLinSys189x49.time_solve_lin_sys                       |
| -        | 29.2±0.4ms                           | 15.4±0.1ms          |    0.53 | solve.TimeSparseSystem.time_linsolve_Aaug(20)                        |
| -        | 55.6±0.4ms                           | 25.0±0.3ms          |    0.45 | solve.TimeSparseSystem.time_linsolve_Aaug(30)                        |
| -        | 28.9±0.1ms                           | 15.4±0.2ms          |    0.53 | solve.TimeSparseSystem.time_linsolve_Ab(20)                          |
| -        | 55.6±0.5ms                           | 24.8±0.1ms          |    0.45 | solve.TimeSparseSystem.time_linsolve_Ab(30)                          |



Full benchmark results can be found as artifacts in GitHub Actions

(click on checks at the top of the PR).

      		
    

  





        


            

            
            

              
            

        

      


      

            
  
  

    
  
    
    

  

  

  


    












      

  
          
 
  

    
  

  



        @oscarbenjamin
  oscarbenjamin



        merged commit 2e489cf
      into

      
  
    sympy:1.13
  


    Jul 19, 2024

  







  

  
  

  
      @oscarbenjamin
  oscarbenjamin


    
    deleted the
    
      
        pr_matrix_inverse_exact_113
    
    branch

    July 19, 2024 08:31    
    




  

  

    
  


  

      


    @jfbu
jfbu



    mentioned this pull request
    
      Aug 4, 2024
    





  


  

      
     Merged

  







  













  
  

    

      

  




      


    


    

      

        
      

  
    Sign up for free
    to join this conversation on GitHub.
    Already have an account?
    Sign in to comment


  



      

    

  




  
          



      

  

    
    

    

  

    Reviewers
  


    

    No reviews






    

  


      
  

    Assignees
  



        


    No one assigned







      


  


  

    Labels
  



    

      

    domainmatrix

  

    floats

  

    matrices

  

    regression









      

  

    

        

    Projects
  


        




    None yet




  



      


  

    
  

    Milestone
  


      
    
  
    SymPy 1.13
  





    
      
          


  

    

      
        

          
  

    Development
  



            


Successfully merging this pull request may close these issues.





  
  



      
    

  




      

    

  
      

  

    

      2 participants
    

    

        
          @oscarbenjamin 
        
          @sympy-bot