CG singular failure seems to be tested on non psd/nsd matrices

[adconner]

I'm writing tests for a pivoted QR implementation, and I noticed that CG is tested in test_iterative_singular in test_singular.py. The issue is that construct_singular_matrix does not respect tags, in particular it doesn't produce psd and nsd matrices (not even symmetric, eg). One could consider this a bug in construct_singular_matrix, but this testing of CG seems to be the only location where singular matrices respecting tags is requested.

As a secondary matter, it seems we are testing that the iterative solvers surely fail on singular inputs. Is this really a guarantee these solvers provide? Also why do we test on special case constructed matrices only in construct_singular_matrix?

[patrick-kidger]

Your first point – fixed in #214! Thanks for spotting and reporting this.

As a secondary matter, it seems we are testing that the iterative solvers surely fail on singular inputs. Is this really a guarantee these solvers provide?

Interesting point. The original plan for Lineax was that a solver would either return the correct answer or error, although as you've highlighted in a previous thread that this is perhaps overly ambitious in general. So this test originated with the desire to check that we would in fact get the desired behaviour of an error on singular matrices. I no longer recall whether that's considered a formal guarantee for the solvers though... it might just be a test we stood up and noticed worked empirically.

Also why do we test on special case constructed matrices only in construct_singular_matrix?

I'm not sure I understand the question – do why not generate a 'general' singular matrix? (However 'general' is defined?)

[adconner]

Also why do we test on special case constructed matrices only in construct_singular_matrix?

I'm not sure I understand the question – do why not generate a 'general' singular matrix? (However 'general' is defined?)

I meant why did we have the line

    if isinstance(solver, (lx.Diagonal, lx.CG, lx.BiCGStab, lx.GMRES)):
        return tuple(matrix.at[0, :].set(0) for matrix in matrices)

I wondered in particular if it was a hack to get the tests to pass.

[patrick-kidger]

Ah, that is because those solvers require square matrices, and the other methods of producing singularity give rectangular matrices.

[adconner]

Oh I see, and in the other case sometimes we emit full rank nonsquare matrices. So when you write "singular" you don't mean "rank deficient"? That's a bit surprising of a usage to me. Would it be worthwhile to be more explicit about our case classes? There are just 3: full rank square, full rank nonsquare, rank deficient (corresponding to the values of the well_posed argument)

[patrick-kidger]

Ah, the terminology in the tests dates from a prerelease version of Lineax where we were categorizing operators into 'nonsingular' and 'singular'. Which of course isn't really correct, and is the reason why in most of the Lineax documentation and source code, we swapped before release to instead saying things like 'well posed' or 'underdetermined' instead. ('Singular' still creeps in in a few places when we're talking about specifically square matrices, and there is no ambiguity.)

I don't think we ever updated the terminology in the tests to reflect that though, and I'm accidentally letting that creep into my language here; you're right to call out the abuse of terminology. I agree with your trifurcation as being the better way to describe things. (We could probably tweak our use of this language in some places in the code.)