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

Added initial docs for SingularException #24192

Closed
wants to merge 9 commits into from
Closed

Added initial docs for SingularException #24192

wants to merge 9 commits into from

Conversation

Vaibhavdixit02
Copy link
Contributor

I have added elementary docs for the SingularException.

@Vaibhavdixit02
Copy link
Contributor Author

@mschauer I was unable to understand the different cases of exception thrown, example - Base.LinAlg.SingularException(1), Base.LinAlg.SingularException(2) and so on and couldn't find any reference to the reasoning behind it. And sorry for taking so much time in doing this.

@mschauer
Copy link
Contributor

Ref #24043.

Don't worry, that is why I wrote that writing a good docstring for these exceptions is perhaps not a good first task, #24043 (comment) . On the technical side: A documented exception (example: DimensionMismatch) has

@Vaibhavdixit02
Copy link
Contributor Author

Okay, I had followed the syntax of

struct LAPACKException <: Exception

The main problem I am facing is understanding the different versions of a specific Exception example Base.LinAlg.SingularException(1) and Base.LinAlg.SingularException(2). How how to locate all the places that Exception is thrown so that I could figure out its meaning?

@fredrikekre
Copy link
Member

That is the INFO code returned by http://www.netlib.org/lapack/explore-html/d3/d6a/dgetrf_8f_source.html. Specifically:

INFO is INTEGER
= 0:  successful exit
< 0:  if INFO = -i, the i-th argument had an illegal value
> 0:  if INFO = i, U(i,i) is exactly zero. The factorization
      has been completed, but the factor U is exactly
      singular, and division by zero will occur if it is used
      to solve a system of equations.

@mschauer
Copy link
Contributor

To expand: also Julia routines use this convention, see https://github.com/JuliaLang/julia/blob/master/base/linalg/diagonal.jl#L361

@Vaibhavdixit02
Copy link
Contributor Author

@mschauer I went through julia/doc/src/stdlib/linalg.md and none of the Exceptions have been mentioned here so I am not sure if I should add it from my side I have made changes in accordance with the Info parameter and following the syntax of

struct LAPACKException <: Exception

@mschauer
Copy link
Contributor

Exceptions imported into base can i my opinion be listed in base.md.

I don't now if this is clear, just in case: The exceptions given in the issue are not listed in any .md file because they do not have a docstring yet. (The docstring is what you read when entering

julia> ?DimensionMismatch

  DimensionMismatch([msg])

  The objects called do not have matching dimensionality. Optional argument
  msg is a descriptive error string.

defined here

(that was linked above).

@Vaibhavdixit02
Copy link
Contributor Author

@mschauer this is what you are looking for?

if ex.info > 0
print(io,"matrix is singular")
elseif ex.info <0
print(io,"argument $(ex.info) had illegal value")
Copy link
Member

Choose a reason for hiding this comment

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

This will never happen; we throw a LAPACK exception before we reach this.

Copy link
Member

Choose a reason for hiding this comment

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

Might as well test for it anyway – defensive programming and all that.

Copy link
Contributor

Choose a reason for hiding this comment

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

That test should go into the constructor.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

If it will always be >0 so I think there is no need for checking at all?

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 we should leave it out here. We check all values before calling LAPACK so we should not even get negative values back from LAPACK and as mentioned by @fredrikekre, we check the LAPACK return values so this is already asserted.

The matrix passed is a singular matrix. Optional argument INFO is an INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: if INFO = i, U(i,i) is exactly zero. The factorization
Copy link
Contributor

Choose a reason for hiding this comment

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

Only this case (>0) should occur, see below.

= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: if INFO = i, U(i,i) is exactly zero. The factorization
has been completed, but the factor U is exactly
Copy link
Contributor

Choose a reason for hiding this comment

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

U is specific to the decomposition, this exception occurs in cases where there is no U.

function Base.showerror(io::IO, ex::SingularException)
print(io, "SingularException: ")
if ex.info > 0
print(io,"matrix is singular")
Copy link
Member

Choose a reason for hiding this comment

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

I'm wondering if we really need a custom show method. In particular when the exception gets a doc string. If we really want a custom show method then we should at least print the value of info since it tells you where the factorization failed.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I had put in the custom show method because I found one for ARPACKException. SingularException already does print info if that is what you mean, but adding clarity regarding the significance of info might be helpful I think.

Copy link
Member

Choose a reason for hiding this comment

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

With this show method, the info field is not shown. The positive value can mean slightly different things depending on the factorization so unless we want a fairly long error message it might be better to refer to the doc string for details.

@mschauer
Copy link
Contributor

@andreasnoack Is it true or would it be a useful convention that info = n in translates into A[1:i,1:i] is singular for i >= n?

@andreasnoack
Copy link
Member

andreasnoack commented Oct 18, 2017

Almost but no. The problem is that factorizations might fail when pivoting is disabled and we still use the same exception, e.g.

julia> inv(lufact([0 1; 1 1], Val{false}))
ERROR: Base.LinAlg.SingularException(1)
Stacktrace:
 [1] inv! at ./linalg/lu.jl:308 [inlined]
 [2] inv(::Base.LinAlg.LU{Float64,Array{Float64,2}}) at ./linalg/lu.jl:310

even though the matrix is non-singular so SingularException is slightly misleading.

Update: Now the sentence might make sense.

@kshyatt kshyatt added docs This change adds or pertains to documentation linear algebra Linear algebra labels Oct 18, 2017
Else if info == -9: There was an error return from calculation of eigenvectors.
Else if info ==-14: Did not find any eigenvalues to sufficient accuracy.
Try with a different starting vector or more Lanczos vectors by increasing the value of ncv.
Else: Uspecified ARPACK error.
Copy link
Member

Choose a reason for hiding this comment

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

I don't think this information is accurate in all cases. It's been awhile since I looked at what we're doing in the wrappers, but at some point all we were doing was wrapping the raw return codes from ARPACK, and the meanings of those change across the different low level routines being wrapped underneath.

Consulting the User Guide, DSAUPD (used for real positive semi-definite symmetric matrices) has these error codes

c  INFO    Integer.  (INPUT/OUTPUT)
c          If INFO .EQ. 0, a randomly initial residual vector is used.
c          If INFO .NE. 0, RESID contains the initial residual vector,
c                          possibly from a previous run.
c          Error flag on output.
c          =  0: Normal exit.
c          =  1: Maximum number of iterations taken.
c                All possible eigenvalues of OP has been found. IPARAM(5)  
c                returns the number of wanted converged Ritz values.
c          =  2: No longer an informational error. Deprecated starting
c                with release 2 of ARPACK.
c          =  3: No shifts could be applied during a cycle of the 
c                Implicitly restarted Arnoldi iteration. One possibility 
c                is to increase the size of NCV relative to NEV. 
c                See remark 4 below.
c          = -1: N must be positive.
c          = -2: NEV must be positive.
c          = -3: NCV must be greater than NEV and less than or equal to N.
c          = -4: The maximum number of Arnoldi update iterations allowed
c                must be greater than zero.
c          = -5: WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.
c          = -6: BMAT must be one of 'I' or 'G'.
c          = -7: Length of private work array WORKL is not sufficient.
c          = -8: Error return from trid. eigenvalue calculation;
c                Informational error from LAPACK routine dsteqr.
c          = -9: Starting vector is zero.
c          = -10: IPARAM(7) must be 1,2,3,4,5.
c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
c          = -12: IPARAM(1) must be equal to 0 or 1.
c          = -13: NEV and WHICH = 'BE' are incompatable.
c          = -9999: Could not build an Arnoldi factorization.
c                   IPARAM(5) returns the size of the current Arnoldi
c                   factorization. The user is advised to check that
c                   enough workspace and array storage has been allocated.

whereas ZNAUPD (used for the complex generalized eigenproblem) has these info codes:

c  INFO    Integer.  (INPUT/OUTPUT)
c          If INFO .EQ. 0, a randomly initial residual vector is used.
c          If INFO .NE. 0, RESID contains the initial residual vector,
c                          possibly from a previous run.
c          Error flag on output.
c          =  0: Normal exit.
c          =  1: Maximum number of iterations taken.
c                All possible eigenvalues of OP has been found. IPARAM(5)  
c                returns the number of wanted converged Ritz values.
c          =  2: No longer an informational error. Deprecated starting
c                with release 2 of ARPACK.
c          =  3: No shifts could be applied during a cycle of the 
c                Implicitly restarted Arnoldi iteration. One possibility 
c                is to increase the size of NCV relative to NEV. 
c                See remark 4 below.
c          = -1: N must be positive.
c          = -2: NEV must be positive.
c          = -3: NCV-NEV >= 2 and less than or equal to N.
c          = -4: The maximum number of Arnoldi update iteration 
c                must be greater than zero.
c          = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'
c          = -6: BMAT must be one of 'I' or 'G'.
c          = -7: Length of private work array is not sufficient.
c          = -8: Error return from LAPACK eigenvalue calculation;
c          = -9: Starting vector is zero.
c          = -10: IPARAM(7) must be 1,2,3.
c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
c          = -12: IPARAM(1) must be equal to 0 or 1.
c          = -9999: Could not build an Arnoldi factorization.
c                   User input error highly likely.  Please
c                   check actual array dimensions and layout.
c                   IPARAM(5) returns the size of the current Arnoldi
c                   factorization.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Oh! The issue #24043 was regarding inline documentation of the errors, adding these many error codes to it might make it too bulky I think. Some kind of work around by adding a link to http://www.caam.rice.edu/software/ARPACK/UG/ug.html might be better?

Copy link
Member

Choose a reason for hiding this comment

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

See also #24292

Copy link
Contributor Author

Choose a reason for hiding this comment

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

A decision for #24292 would provide a clearer way then!


Optional Argument info is an Integer.
If info == -1: The matrix is not Hermitian.
Else: The matrix is not positive definite.
Copy link
Member

Choose a reason for hiding this comment

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

This is close but not quite correct since a zero value would mean that the matrix is PD. It would probably also be good to mention that this was tested by a Cholesky factorization, i.e. something like "if info>0 then the Cholesky factorization failed at the $(info) diagonal element which indicates that the matrix is not positive definite."


The matrix passed is a singular matrix. Optional argument INFO is an INTEGER > 0
if INFO = i, The matrix is exactly singular, and division by zero
will occur if it is used to solve a system of equations.
Copy link
Member

Choose a reason for hiding this comment

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

We'll need to add here that when pivoting is disabled for a factorization which exception might not indicate that matrix is singular. E.g. the matrix [0 1l 1 1] will raise this exception for LU without pivoting. Also, please avoid copying the Fortran variable and type names from the LAPACK documentation. INTEGER is not really a thing in Julia and the variable in Julia is info not INFO.

Copy link
Contributor Author

@Vaibhavdixit02 Vaibhavdixit02 Nov 7, 2017

Choose a reason for hiding this comment

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

I want to cross check if I understood it correctly, "if pivoting is disabled the value of info is irrelevant to the singularity of the matrix" is this what you meant? Sorry for not paying attention to INTEGER !

struct SingularException <: Exception
info::BlasInt
end

function Base.showerror(io::IO, ex::SingularException)
Copy link
Member

Choose a reason for hiding this comment

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

We'd need to add similar information here about the case where pivoting is disabled.

@@ -0,0 +1,284 @@
{
Copy link
Member

Choose a reason for hiding this comment

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

This file seems like it got accidentally added.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yeah! 😅

@dkarrasch
Copy link
Member

Seems to be fixed by #31251.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
docs This change adds or pertains to documentation linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

9 participants