Refactor algorithm selection logic

[mtfishman]

This refactors the algorithm selection logic (select_algorithm) in terms of two layers. In this new design, there is a generic function select_algorithm for any factorization function and matrix type that processes the alg keyword argument, and either uses the algorithm that is specified or if none is specified forwards to a lower level generic function default_algorithm, which then chooses a default algorithm depending on the factorization function and matrix type.

The motivation for this was that I started making a PR to make select_algorithm stricter about how it handles keyword arguments, but it was cumbersome to do since the select_algorithm logic for processing keyword arguments was repeated for each factorization function. In addition, I added logic for allowing an algorithm type (in addition to a Symbol and algorithm object) to be passed as alg, which was also easier to implement in the new design.

@lkdvos and I were discussing the naming of these two layers of functions (currently select_algorithm and default_algorithm), he proposed select_algorithm -> _select_algorithm and default_algorithm -> select_algorithm. I'm definitely open to suggestions, curious what you think @Jutho.

To-do:

• Update the docstrings and add docs for default_algorithm.
• Dedicated tests for select_algorithm/default_algorithm.
• Investigate type stability.
• Decide on the names.

[codecov]

Codecov Report

Attention: Patch coverage is 94.31818% with 5 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/implementations/truncation.jl	70.00%	3 Missing ⚠️
src/algorithms.jl	89.47%	2 Missing ⚠️

Files with missing lines	Coverage Δ
src/MatrixAlgebraKit.jl	100.00% <ø> (ø)
src/implementations/orthnull.jl	96.77% <100.00%> (ø)
src/interface/eig.jl	73.33% <100.00%> (+16.19%)	⬆️
src/interface/eigh.jl	73.33% <100.00%> (+16.19%)	⬆️
src/interface/lq.jl	100.00% <100.00%> (+36.36%)	⬆️
src/interface/polar.jl	63.63% <100.00%> (+7.38%)	⬆️
src/interface/qr.jl	100.00% <100.00%> (+36.36%)	⬆️
src/interface/schur.jl	50.00% <100.00%> (+26.92%)	⬆️
src/interface/svd.jl	73.33% <100.00%> (+16.19%)	⬆️
src/algorithms.jl	82.43% <89.47%> (+0.46%)	⬆️
... and 1 more

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[mtfishman]

What this PR calls _select_algorithm was called algorithm_or_select_algorithm in #19, that is at least a more descriptive name for what it is doing.

EDIT: So to clarify, the proposal would be to change select_algorithm -> algorithm_or_select_algorithm or maybe_select_algorithm and default_algorithm -> select_algorithm.

[lkdvos]

Looks like a reasonable change to me, should definitely offer quite a bit of flexibility.

If we decide on this, can you also look at exporting default_algorithm, since that would be the main point to overload, and check if the docs need updating for this new approach?

[mtfishman]

If we decide on this, can you also look at exporting default_algorithm, since that would be the main point to overload, and check if the docs need updating for this new approach?

What do you think of making it public but not exporting it for now, indicating it is more of a syntax for overloading rather than using externally? (EDIT: I made default_algorithm and select_algorithm public for now but I'm happy to change that.)

[mtfishman]

Currently, if you pass an algorithm type or Symbol as the alg keyword argument, select_algorithm and therefore the factorizations are not type stable. I thought constant propagation would handle that but I guess not. I tried playing around with @inline but it didn't seem to help.

[lkdvos]

I don't have any strong feelings about public or exported in this case. Given that our use-cases are mostly to use this as a low-level library to build on, I don't think it matters too much. As you noticed public doesn't work in 1.10, I think the recommended way of resolving that is using Compat.jl to handle this instead of doing that manually.

I'm a bit more worried about the type stability though.
Was this type-stable before?
There's also Base.@constprop :aggressive as an option, or explicitly Val-ing some arguments in the internal implementations.

[mtfishman]

I don't have any strong feelings about public or exported in this case. Given that our use-cases are mostly to use this as a low-level library to build on, I don't think it matters too much. As you noticed public doesn't work in 1.10, I think the recommended way of resolving that is using Compat.jl to handle this instead of doing that manually.

There are long debates about that: https://discourse.julialang.org/t/is-compat-jl-worth-it-for-the-public-keyword/119041, the approach I use here is the one recommended in the Julia docs: https://docs.julialang.org/en/v1.12-dev/manual/modules/#Export-lists. Some people are weary of Compat.jl because it causes a lot of invalidations and uses type piracy, so if there is a simple way to avoid it it seems better to do so, but I'm fine either way.

I'm a bit more worried about the type stability though. Was this type-stable before? There's also Base.@constprop :aggressive as an option, or explicitly Val-ing some arguments in the internal implementations.

I doubt it was type stable before but I can double check. I also tried Base.@constprop :aggressive but that didn't seem to help, though I can explore that a bit more, and I can also look into if some use of Val can help.

[mtfishman]

I think basically what we are doing is analogous to this minimal example:

julia> f(; x) = Val{x}()
f (generic function with 1 method)

julia> @code_warntype f(; x=:x)
MethodInstance for Core.kwcall(::@NamedTuple{x::Symbol}, ::typeof(f))
  from kwcall(::NamedTuple, ::typeof(f)) @ Main REPL[56]:1
Arguments
  _::Core.Const(Core.kwcall)
  @_2::@NamedTuple{x::Symbol}
  @_3::Core.Const(Main.f)
Locals
  @_4::Symbol
  x::Symbol
Body::Val
1 ─       Core.NewvarNode(:(@_4))
│   %2  = Core.isdefined(@_2, :x)::Core.Const(true)
└──       goto #3 if not %2
2 ─       (@_4 = Core.getfield(@_2, :x))
└──       goto #4
3 ─       Core.Const(:(Core.UndefKeywordError(:x)))
└──       Core.Const(:(@_4 = Core.throw(%6)))
4 ┄ %8  = @_4::Symbol
│         (x = %8)
│   %10 = Base.keys(@_2)::Core.Const((:x,))
│   %11 = (:x,)::Core.Const((:x,))
│   %12 = Base.diff_names(%10, %11)::Core.Const(())
│   %13 = Base.isempty(%12)::Core.Const(true)
└──       goto #6 if not %13
5 ─       goto #7
6 ─       Core.Const(:(Base.kwerr(@_2, @_3)))
7 ┄ %17 = Main.:(var"#f#1")::Core.Const(Main.var"#f#1")
│   %18 = x::Symbol
│   %19 = (%17)(%18, @_3)::Val
└──       return %19

julia> @code_warntype (() -> f(; x=:x))()
MethodInstance for (::var"#5#6")()
  from (::var"#5#6")() @ Main REPL[61]:1
Arguments
  #self#::Core.Const(var"#5#6"())
Body::Val{:x}
1 ─ %1 = (:x,)::Core.Const((:x,))
│   %2 = Core.apply_type(Core.NamedTuple, %1)::Core.Const(NamedTuple{(:x,)})
│   %3 = (:x,)::Core.Const((:x,))
│   %4 = (%2)(%3)::Core.Const((x = :x,))
│   %5 = Core.kwcall(%4, Main.f)::Core.Const(Val{:x}())
└──      return %5

julia> using TestExtras

julia> @constinferred f(; x=:x)
Val{:x}()

julia> @constinferred f(; x=Symbol("x"))
Test Failed at REPL[75]:1
  Expression: @constinferred f(; x = Symbol("x"))
   Evaluated: Val{:x} != Val

ERROR: There was an error during testing

so it seems like at least in principle it should be inferrable.

[Jutho]

I am not sure using multiple dispatch via _select_algorithm to deal with the different cases of alg keywords is the right approach, as it might put more strain on the compiler. I don't think there is need to be able to extend this functionality, so something like

function select_algorithm(f, A; alg=nothing, kwargs...)
	if isnothing(alg)
		return default_algorithm(f, A; kwargs...)
	elseif alg isa AbstractAlgorithm
		isempty(kwargs) ||
			throw(ArgumentError("Additional keyword arguments are not allowed when an algorithm is specified."))
		return alg
	elseif
	...
	end
end

should work. Or does the _select_algorithm lowering actually help with the lowering?

[mtfishman]

I'm totally fine with that, using dispatch was a bit arbitrary. I generally like using dispatch rather than if-statements purely as a style thing, but if you and the compiler prefer if-statements I can switch to that.

[mtfishman]

I'm trying this out locally, an unfortunate side affect of this is that functions like left_orth/right_orth become type unstable, since before they were calling _select_algorithm which is type stable and select_algorithm isn't with the new design (I think what happens is that certain branches of select_algorithm aren't type stable, so then the entire thing becomes type unstable even for cases that used to be type stable).

I'll investigate that more. But maybe it would be better to just change select_algorithm to accept alg as a positional argument and go back to using dispatch. That should solve all type stability problems being discussed above as well.

[Jutho]

If the _select_algorithm lowering is better for type stability, I am certainly fine with that, and there is no need to waste time trying to trick the compiler to make the if elseif else construction inferable as well.

[Jutho]

Is the idea that if we only modify this list for the eig_full tests, this will test the selection/default functionality and should thus also work for the other factorisations? I am fine with this, and definitely think we should not run all tests 3 times over as this would be a large amount of unnecessary compute. Just checking if this is the plan?

[mtfishman]

For now I just added it here to get some more feedback on the design before going through the work of adding it to all of the factorization tests, and also as you say I also wanted to hear your opinion on whether we even want to test all of the factorizations.

Maybe I could add it to one more factorization test, and also add some standalone tests of select_algorithm and default_algorithm.

[Jutho]

Standalone select_algorithm and default_algorithm tests are definitely a good idea. Regarding the names: I like those names 😄 !

[lkdvos]

Regarding the type stability, I just remembered/realized that this change has an additional side-effect that might affect this: before we didn't have two layers (select_algorithm and default_algorithm), and as a result the implementations were written directly in terms of select_algorithm(::typeof(f), ...), while now there are a bunch of definitions for select_algorithm(f, ...). Given that Julia doesn't specialize on functions by default, I'm wondering if this might be affecting the outcomes. In other words, I'm wondering if it's not better to write select_algorithm(f::F, ...) where {F} if we want to have these two layers of functions in there.

[mtfishman]

Ok, I've solved most type instability issues and added more tests. I had to rely both on forcing method specialization and Base.@constprop :aggressive as @lkdvos suggested.

Oddly, the one case that still isn't type stable is specifying alg as a type, i.e. svd_compact!(randn(2, 2); alg=LAPACK_QRIteration). Any suggestions with that would be welcome...

[mtfishman]

This is ready for review. All type stability issues are fixed, except for one issue in Julia 1.10, I would vote for letting that slide (as opposed to making the code more complicated to fix it) but please let me know if you want me to look into that.

In the latest, I was able to avoid all uses of Base.@constprop :aggressive. @lkdvos pointed out to me that I was using @constinferred wrong and should have been interpolating the alg input in the tests.

[mtfishman]

This is the case I mentioned that isn't type stable in older versions of Julia. I'd vote for letting this slide and telling users to use a newer versions of Julia if they care about type stability, let me know what you think.

[lkdvos]

I'm assuming that this is also just some strange interplay with the non-heterogeneous tuple of algs, so this is more of a test artifact than a real problem. Definitely okay with letting this slide.

[lkdvos]

I only looked over this briefly on my phone so it would be good if @Jutho can also give his blessing, but this looks like a nice implementation now, I'm very happy the constprop annotations weren't necessary in the end.

[Jutho]

Very nice PR, only pending comment is a suggestion in an attempt to further improve the clarity of the select_algorithm doc string.

[Jutho]

This looks ready to be merged! Feel free to merge yourself @mtfishman , otherwise without reaction I will merge this tonight.

[mtfishman]

Thanks again for the careful review.

Should we register v0.2 after #25 gets merged?

[Jutho]

Yes that sounds good to me. I went through #25 and do not have any comments.