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: General 1-norm estimation for some functions of sparse matrices #12609

Closed
wants to merge 2 commits into from

Conversation

marcusps
Copy link
Contributor

@mfasi @tkelman

Here is what I had in mind for a more general normest routine. At the moment there is something wrong with Ac_ldiv_B! that I don't quite understand (something about the types keeps this from compiling).

This attempt of a general normest is also not type stable (the pre routine in particular mucks that up) but that is probably fixable by replacing it with an object instead of a function.

I also put normest in a separate file just to make it easier to test it while I play with it. This pull request is really just meant to bring attention to the idea and discuss the proper way to do it.

EDIT:

This is motivated by the usefulness of such a generic function for the fast and accurate computation of f(A)*v without explicitly computing f(A), which is especially attractive when A is sparse but f(A) may not be. The one interesting case where this can be applied is expm, but potentially others as well, such as logm, sqrtm (see below).

Computing the Action of the Matrix Exponential, with an Application to Exponential Integrators, Awad H. Al-Mohy and Nicholas J. Higham, SIAM Journal on Scientific Computing 2011 33:2, 488-511. (preprint)

Computing f(A)b for Matrix Functions f, Philip I. Davies and Nicholas J. Higham. Lecture Notes in Computational Science and Engineering Volume 47, 2005, pp 15-24 (preprint)

Estimating the condition number of f(A)b, Edvin Deadman. Numerical Algorithms December 2014 (article)

An Improved Schur--Padé Algorithm for Fractional Powers of a Matrix and their Fréchet Derivatives, Nicholas J Higham, Lijing Lin, SIAM. J. Matrix Anal. & Appl. 32 1056 2013. (preprint)

@kshyatt kshyatt added the sparse Sparse arrays label Aug 13, 2015
+ Add tests
+ Remove check for parallel columns in complex matrices
+ Remove redundant capping of number of blocks
+ Add more robust default value for number of blocks
+ Minor change in complex "sign" for type stability
@marcusps
Copy link
Contributor Author

I have changed my approach slightly, and the generalization seems to work (this branch includes some minor changes to @mfasi's implementation of Higham and Tisseur's algorithm for norm estimation, but this first commit with the minor changes has been submited as a separate pull request).

I had to change how cond() and normestinv() are tested, because the algorithm implemented is probabilistic and only provides upper and lower bounds to these two quantities. Other than that, it should be easy to do 1 norm estimation for matrix powers (needed in ExpmV.jl) as well as other functions of sparse matrices.

(a-b)/b
end

test_normestinv_Ac = Float64[begin
Copy link
Contributor

Choose a reason for hiding this comment

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

trailing whitespace on these - make sure they pass make check-whitespace locally

@marcusps marcusps force-pushed the normest-general branch 2 times, most recently from 3714a9c to 6ec4d0a Compare August 19, 2015 01:52
@marcusps marcusps changed the title WIP: General 1-norm estimation for sparse matrix function RFC: General 1-norm estimation for sparse matrix function Aug 19, 2015
@marcusps marcusps changed the title RFC: General 1-norm estimation for sparse matrix function RFC: General 1-norm estimation for some functions of sparse matrices Aug 19, 2015
@marcusps
Copy link
Contributor Author

@mfasi I have a fully working version now. I am curious to hear what you think. There maybe a slight performance hit with normestinv() written this way, but I think it should be minor (although I have not measured).

@KristofferC
Copy link
Member

Before you complicate the code with @eval I would suggest benchmarking to see if a significant overhead is induced. Maybe for smallish matrices it is relevant.

@tkelman
Copy link
Contributor

tkelman commented Oct 14, 2015

Any updates here? Did other later PR's cover this functionality?

@marcusps
Copy link
Contributor Author

No, not yet. I haven't had time to do careful benchmarking of this code vs what @mfasi suggested.

@marcusps
Copy link
Contributor Author

I have some code elsewhere that does, but I have not delved into detailed
benchmarking yet — I just haven't had time.

On Wed, Oct 14, 2015 at 4:24 PM Tony Kelman notifications@github.com
wrote:

Any updates here? Did other later PR's cover this functionality?


Reply to this email directly or view it on GitHub
#12609 (comment).

@GiggleLiu
Copy link
Contributor

why this pr is not merged yet?

@ViralBShah
Copy link
Member

Would be nice to revive this, but it is also really old and may just need a fresh start.

@marcusps
Copy link
Contributor Author

marcusps commented Dec 20, 2018 via email

@ViralBShah
Copy link
Member

Enough has changed since this was created that we'll need to substantially redo this if we want it.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
sparse Sparse arrays
Projects
None yet
Development

Successfully merging this pull request may close these issues.

7 participants