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

WIP: Add specialized methods for inverses of small matrices #12454

Closed
wants to merge 1 commit into from

Conversation

KristofferC
Copy link
Member

Same as #12452 but for inverses

The following benchmark shows the difference between master and PR (using #12452):

function f(N, n)
    A = rand(n, n)
    @time for i = 1:N
        inv(A)
    end
end

Master

f(10^5, 2)
#  0.094840 seconds (800.00 k allocations: 48.828 MB, 5.10% gc time) 
f(10^5, 3)
#  0.124090 seconds (800.00 k allocations: 59.509 MB, 4.77% gc time)

PR

f(10^5, 2);
#0.009085 seconds (100.00 k allocations: 9.155 MB)
f(10^5, 3);
#0.012562 seconds (100.00 k allocations: 13.733 MB)

function small_inv{T}(A::StridedMatrix{T})
n = chksquare(A)
d = det(A)
if d == zero(T)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is exact comparison with zero reasonable here?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably not. I thought about it but I didn't know what would be reasonable due to the generic T. eps(T) maybe?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it is fine to throw for exact zeros only. That is similar to what LAPACK does.

@andreasnoack
Copy link
Member

These functions actually assume commutativity for multiplication of the elements so the elements should probably be restricted to subtypes of Union{Real,Complex}.

@KristofferC
Copy link
Member Author

I restricted the types for when the small_inv method is called now.

@KristofferC
Copy link
Member Author

Could someone restart the AppVeyor build for me?

@tkelman
Copy link
Contributor

tkelman commented Aug 6, 2015

You need to rebase and force push here, you hit a build-number collision appveyor/ci#353

@KristofferC
Copy link
Member Author

To note, the timings here are with the hand coded determinants in #12452. These might be unacceptable due to numerical instability. WIth only #12460 the timings in this PR are:

f(10^5, 2);
#  0.045274 seconds (700.00 k allocations: 41.199 MB, 8.82% gc time)
f(10^5, 3);
#  0.052139 seconds (700.00 k allocations: 51.880 MB, 6.88% gc time)

which is only about a factor 2x faster.

@StefanKarpinski
Copy link
Member

This looks pretty reasonable. @andreasnoack, please review and merge if you think this is good to go.

@andreasnoack
Copy link
Member

There are two issues here.

  1. I think most of the speed up will depend on the faster determinant in WIP: Add specialized methods for determinant of small matrices #12452, and they have a larger error so I'd like to find another solution for them.
  2. Even with more precise determinants, I think these inverses might be imprecise relative to inverses based on the LU. I'll have to run some tests to confirm that. @KristofferC have you looked into the errors?

@mlubin
Copy link
Member

mlubin commented Aug 12, 2015

Couldn't you hand-code a 2x2 LU, for example, to work around the numerical stability issue?
That seems to be the motivation behind the proposal for the modified hand-coded determinant also (#12452 (comment)).

@KristofferC KristofferC changed the title Add specialized methods for inverses of small matrices WIP: Add specialized methods for inverses of small matrices Aug 12, 2015
@andreasnoack
Copy link
Member

@mlubin Sure. I think that is the way to procees, but the LU based inverse is a little verbose. The 2x2 is okay, but the 3x3 case is quite messy and has six branches.

You can try out inv(SymPy.Sym[Sym("a$i$j") for i = 1:3, j = 1:3])

@KristofferC
Copy link
Member Author

Closing this because a user can easily add their own hand coded inverses if they feel the default speed is not fast enough. Better to provide numerically stable functions by default.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants