Fix kron ambiguity for triangular/hermitian matrices

[dkarrasch]

This resolves method ambiguities in SparseArrays.jl, that have been introduced in #53186. As usual, we would want the result to be sparse if any one Kronecker factor is sparse. In that case, however, we would need a specialized kernel to begin with. So, this is more like a quick fix.

This has not been released, so it's not a regression.

@araujoms Is this fine with you?

[araujoms]

There are some problems. First is that you removed the restriction to a subtype of Number. This reintroduces the bug that @jishnub had pointed out here: #53186 (comment) . For ease of reference, the problem is when we have a matrix of matrices:

m = UpperTriangular(reshape(fill(ones(2,2),4), 2, 2))
n = Symmetric(reshape(fill(ones(2,2),4), 2, 2))
kron(m,m)
kron(n,n)

The second is that you removed the parametrization of the functions. I don't understand why, we really want the compiler to specialize on conj and real.

The third is that now if we wrap Hermitian around another wrapper, like Diagonal, we won't have the nice code anymore, but the generic fallback. Of course it would be much more efficient to write something directly for Hermitian(Diagonal), but I was relying on the kronecker products staying conveniently Hermitian. I won't write code for every possible double wrapper, but I might do it for Diagonal.

[dkarrasch]

First is that you removed the restriction to a subtype of Number.

Hm, I see. Dispatching on type parameters is just not recommended (especially when we have so many existing methods already), since it's prone to introducing ambiguities and increases latency.

The second is that you removed the parametrization of the functions. I don't understand why, we really want the compiler to specialize on conj and real.

IIRC, then this is what the compiler does: it specializes on function arguments whenever they are called in the function body, and not merely passed on as an argument. I'll need to double-check.

The third is that now if we wrap Hermitian around another wrapper, like Diagonal, we won't have the nice code anymore, but the generic fallback.

Sure, but now we have stdlibs in an ambiguous state. So we can't kron two symmetric/triangular sparse matrices, not even with some generic fallback. The other option then would be to work out the analogue of #53186 for symmetric/hermitian/triangular sparse matrices, before you get to the Hermitian(Diagonal) case.

[araujoms]

Hm, I see. Dispatching on type parameters is just not recommended (especially when we have so many existing methods already), since it's prone to introducing ambiguities and increases latency.

It's better than introducing a bug, though. If you can fix the underlying bugs to make kron work with the recursive matrices, great, we can remove the specialization to <: Number. I can't.

IIRC, then this is what the compiler does: it specializes on function arguments whenever they are called in the function body, and not merely passed on as an argument. I'll need to double-check.

Well they are just passed on as an argument, no? I'm pretty sure it wasn't specializing without it.

Sure, but now we have stdlibs in an ambiguous state. So we can't kron two symmetric/triangular sparse matrices, not even with some generic fallback. The other option then would be to work out the analogue of #53186 for symmetric/hermitian/triangular sparse matrices, before you get to the Hermitian(Diagonal) case.

Sure, I was just complaining, we can't let a method ambiguity stay.

[dkarrasch]

Not to insist, but to learn something: I read

https://docs.julialang.org/en/v1/manual/performance-tips/#Be-aware-of-when-Julia-avoids-specializing

as it will specialize since in our case, the Function arguments are called within _hermkron!, not just passed as an argument to another function.

[araujoms]

Fair enough, it is specializing then. I did a quick test to confirm it:

f(x) = x^2
g(x) = abs2(x)
test(z::Real) = do_power(z,f)
test(z::Complex) = do_power(z,g)
do_power(z,h) = h(z)

It wasn't actually me who added this parametrization, but I can't find in the PR who did, so I could ask what they had in mind. I guess it doesn't really matter.

[dkarrasch]

Could well be me (who suggested it) at the time. I remember looking it up not so long ago.

[dkarrasch]

How about this now? The signature could be relaxed again once somebody writes sparse kron-kernels, but for now I think this is not too bad.

[araujoms]

I think it's fine now. I'd suggest adding the four lines I posted above to the tests, so that this bug doesn't get accidentally reintroduced.

[araujoms]

Thanks! I've opened an issue to track the underlying problem.