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

[RFC] Add estimates of operator norms #39058

Closed
wants to merge 77 commits into from

Conversation

sethaxen
Copy link
Contributor

@sethaxen sethaxen commented Dec 31, 2020

This PR adds three methods:

  • opnormest1(A): lower-bound estimate of opnorm(A, 1)
  • opnormestprod1(As...)/opnormestprod1(A,p::Integer): lower-bound estimate of opnorm(prod(As), 1)
  • opnormestinv1(A): lower-bound estimate of opnorm(inv(A), 1)

The implementation of opnormest1 comes from the paper by Higham and Tisseur. It is especially useful when computing the actions of matrix functions (e.g. exp(A)*b) for large and sparse A. Edit: It also works on any object that implements a minimal interface of a linear operator.

There was already an implementation of opnormestinv1 opnormest(inv, A, 1) (called opnormestinv(A)) in SparseArrays which served as the starting point for this implementation and is now superseded. From a quick scan, future PRs could make use of this estimate for computing condition numbers and for speeding up the matrix logarithm, which currently calls opnorm(A^p, 1) for powers of 1 through 5. Even for a dense A, the estimate is faster and on my machine gives a 20% speed-up for a large matrix. Edit: this implementation also minimizes the number of allocations, which causes a substantial speed-up

I've marked this as RFC because I'm not very satisfied with the interface thus far, in particular the interface to return vectors corresponding to the estimate (see docstrings).

Edit: I've added the opnormest function which is exported. The public interface is now:

  • opnormest(A, p=2): estimate opnorm(A, p)
  • opnormest(f, A, p=2): estimate opnorm(f(A), p) for certain f (pinv, inv, and prod already supported). While clean, this interface is not very extensible. e.g., if we want to support f(A) = A^n or f(A) = A^(-n), there's not an obvious existing f we can use (except maybe Base.Fix2, but we would need to add ^(n::Integer) = Base.Fix2(^, n) to make it user-friendly).

I also added a power iteration method for estimating the 2-norm, which I believe is the approach matlab's normest function uses.

Relates #12609, JuliaLang/LinearAlgebra.jl#232

@sethaxen sethaxen changed the title [RFC] Add estimates of operator 1-norm [RFC] Add estimates of operator norms Jan 3, 2021
end
@test_throws DimensionMismatch SparseArrays.opnormestinv(sprand(3,5,.9))
@test_throws DimensionMismatch (@test_deprecated SparseArrays.opnormestinv(sprand(3,5,.9)))
Copy link
Contributor Author

Choose a reason for hiding this comment

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

These currently error. It doesn't seem that chaining @test_deprecated for these tests actually works.

@DilumAluthge
Copy link
Member

DilumAluthge commented Jan 14, 2022

@sethaxen It looks like this PR touches five files. Three of those files are in the LinearAlgebra stdlib, and the remaining two files are in the SparseArrays stdlib.

We have moved the SparseArrays stdlib to an external repository. Therefore, for the portion of this PR that modifies the SparseArrays stdlib, can you please:

  1. Remove the SparseArrays changes from this PR.
  2. Open a new PR against https://github.com/JuliaLang/SparseArrays.jl with the SparseArrays changes.

Thank you!

@DilumAluthge DilumAluthge removed the sparse Sparse arrays label Jan 14, 2022
@DilumAluthge
Copy link
Member

We have moved the LinearAlgebra stdlib to an external repo: https://github.com/JuliaLang/LinearAlgebra.jl

@sethaxen If you think that this PR is still relevant, please open a new PR on the LinearAlgebra.jl repo.

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

Successfully merging this pull request may close these issues.

4 participants