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

change order of exponent for neg. matrix powers #20071

Closed
wants to merge 1 commit into from

Conversation

felixrehren
Copy link
Contributor

According to Higham-Lin [1], numerical accuracy is higher for the calculation A^-k, when A is a matrix and k a positive integer, if it is calculated as (A^-1)^k.

[1] A Schur-Pade algorithm for fractional powers of a matrix, SIAM J. Matrix Anal. & Appl., Vol 32/3 1056-1078
Section 6:
"Algorithm 6.1 [A^-k = (A^k)^-1] inverts A^k, which is potentially a much more ill conditioned matrix than A. Intuitively, Algorithm 6.2 [A^-k = (A^-1)^k] should therefore be preferred.
Section 9, Experiment 7: "Algorithms 6.2 ... clearly produce much more accurate results than Algorithm 6.1, as we expected."

According to Higham-Lin [1], numerical accuracy is higher for the calculation `A^-k`, for `A` a matrix and `k` an integer, if it is calculated as `(A^-1)^k`

[1] A Schur-Pade algorithm for fractional powers of a matrix, SIAM J. Matrix Anal. & Appl., Vol 32/3 1056-1078
Section 6:
"Algorithm 6.1 [A^-k = (A^k)^-1] inverts `A^k`, which is potentially a much more ill conditioned matrix than `A`. Intuitively, Algorithm 6.2 [A^-k = (A^-1)^k] should therefore be preferred.
Section 9, Experiment 7: "Algorithms 6.2 ... clearly produce much more accurate results than Algorithm 6.1, as we expected."
@kshyatt kshyatt added the linear algebra Linear algebra label Jan 16, 2017
@tkelman
Copy link
Contributor

tkelman commented Jan 16, 2017

For something like a tridiagonal matrix it could be cheaper to calculate A^(-p) first then invert vs inverting first then taking inv(A)^(-p). Not sure whether this conditioning-vs-cost tradeoff should be made on a type by type basis, or prefer more-accurate fallbacks that are only equivalent in cost for dense matrices though.

@andreasnoack
Copy link
Member

I think we should just add a specialized method for Tridiagonal if this is a concern

@fredrikekre
Copy link
Member

Was included in #21184

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.

5 participants