Consistency checks missing when constructing linear algebra objects

Specific Domains Numerics
1

I have the following question which virtually applies to all defined linear algebra objects in Julia/LinearAlgebra:

Why there are no consistency checks when constructing linear algebra objects?

As an example, let’s consider the Schur data type, which is defined as:

struct Schur{Ty,S<:AbstractMatrix} <: Factorization{Ty}
    T::S
    Z::S
    values::Vector
end

This structure contains two matrices T and Z of the same type and a vector values. It has nothing specific related to the Schur decomposition, as claimed by the documentation. This can be easily checked by the following construction:

julia> F = Schur(rand(3,3),rand(2,1),rand(0))
Schur{Float64, Matrix{Float64}}
T factor:
3×3 Matrix{Float64}:
 0.851595  0.622935  0.416998
 0.493628  0.782394  0.54265
 0.522581  0.995593  0.996112
Z factor:
2×1 Matrix{Float64}:
 0.6150198969878795
 0.08655354595467768
eigenvalues:
Float64[]

To avoid such misconstructions, I am missing the following minimal checks: T and Z should be both square matrices of the same order, say n, and values should be a vector of length n. More involved tests could also check that T is quasi-upper triangular in the general real case or upper triangular in the general complex case.

I would like to add that as far as some linear algebra data types serve only to pack the computed results, this setup could be acceptable. However, if the constructed objects are arguments of other functions, then the situation is different. As an example, this lack of checks is disturbing, especially when in a next step I would like to “reorder” the eigenvalues:

ordschur(F,falses(5))
 ** On entry to DTRSEN parameter number  8 had an illegal value
ERROR: ArgumentError: invalid argument #8 to LAPACK call

Here ordchur accepts the above Schur object without complaints, delivers its elements to the trsen! wrapper to the LAPACK routine DTRSEN, and, finally, the inconsistency is discovered only at the level of this routine with an error message which still requires further inspection of the LAPACK code.

The inconsistency could even remain undetected, leading to completely nonsense results:

ordschur(Schur(rand(3,3),rand(3,1),rand(0)),falses(5))
Schur{Float64, Matrix{Float64}}
T factor:
3×3 Matrix{Float64}:
 0.848517   0.385575  0.183215
 0.0834543  0.938042  0.919729
 0.746262   0.886942  0.0102043
Z factor:
3×1 Matrix{Float64}:
 0.3533799437577295
 0.7928718694559094
 0.01924732093261672
eigenvalues:
3-element Vector{ComplexF64}:
  0.8485167187308525 + 0.1793818626275508im
  0.9380416658902896 + 0.9031867257980865im
 0.01020433757434469 - 0.9031867257980865im

I would like to hear other opinions on how much checking is necessary to build only meaningful linear algebra objects.

2 Likes
2

Opened an issue:
Consistency checks missing when constructing linear algebra objects · Issue #44191 · JuliaLang/julia (github.com)

2 Likes
3

LinAlg experts waiting for input:

@stevengj: If we want to have a catch-all issue like this, it really needs to come with a checklist of constructors that need checks, i.e. someone needs to do a survey.

2 Likes
4

You generally shouldn’t construct linear algebra objects by calling the constructors directly and if you do then it’s up to you to ensure that you do it correctly. E.g. by writing your own tests for the functions that end up calling the constructors directly.

However, we do generally try to catch all illegal values in Julia before we call LAPACK so it looks like the checks in our wrapper for xtrsen doesn’t have all the required checks.