Introduce AbstractKroneckerArray

[lkdvos]

What it says in the title :)

[codecov]

Codecov Report

❌ Patch coverage is 79.45205% with 15 lines in your changes missing coverage. Please review.
✅ Project coverage is 81.21%. Comparing base (06d9939) to head (f6a302a).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/kroneckerarray.jl	78.72%	10 Missing ⚠️
src/linearalgebra.jl	81.25%	3 Missing ⚠️
src/matrixalgebrakit.jl	77.77%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #54      +/-   ##
==========================================
- Coverage   81.87%   81.21%   -0.67%     
==========================================
  Files           9        9              
  Lines         618      628      +10     
==========================================
+ Hits          506      510       +4     
- Misses        112      118       +6     

Flag	Coverage Δ
docs	0.00% <0.00%> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

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

[mtfishman]

Good catch, I was a bit lazy with this definition but as you say we would need to be more careful about the tolerances. This is a case where it would be easier to special case for SectorArray.

[mtfishman]

The best I could come up with is:

using LinearAlgebra: promote_leaf_eltypes
function Base.isapprox(
        a::AbstractKroneckerArray, b::AbstractKroneckerArray;
        atol::Real = 0,
        rtol::Real = Base.rtoldefault(promote_leaf_eltypes(a), promote_leaf_eltypes(b), atol),
        norm::Function = norm
    )
    a1, a2 = arg1(a), arg2(a)
    b1, b2 = arg1(b), arg2(b)
    # Approximation of:
    # norm(a - b) = norm(a1 ⊗ a2 - b1 ⊗ b2)
    #             = norm((a1 - b1) ⊗ a2 + b1 ⊗ (a2 - b2) + (a1 - b1) ⊗ (a2 - b2))
    diff1 = norm(a1 - b1)
    diff2 = norm(a2 - b2)
    d = diff1 * norm(a2) + norm(b1) * diff2 + diff1 * diff2
    return iszero(rtol) ? d <= atol : d <= max(atol, rtol * max(norm(a), norm(b)))
end

which would work for SectorArrays since a1 == b1.

[lkdvos]

I'm not sure I'm fully following your derivation, I don't see how you can split the norms like that?
In principle there is the following formula:

$$|| a1 \otimes a2 - b1 \otimes b2 ||^2 = || a1 \otimes a2 ||^2 + || b1 \otimes b2 ||^2 - 2 Re \langle a1 \otimes a2, b1 \otimes b2 \rangle$$

But this doesn't behave well with finite precision, since we are basically subtracting equal magnitude numbers to attempt to find something of much smaller magnitude

[mtfishman]

I was realizing that if we are using these functions outside of KroneckerArrays.jl, these names are a bit too general. We can of course call them as KroneckerArrays.arg1(a), etc. explicitly, or rename them to something more explicit like kron_arg1(a), etc. Any opinion on that?

[mtfishman]

I guess arg1 and arg2 are also used for non-Kronecker data structures such as CartesianProduct so maybe kron_arg* isn't the best name.

[lkdvos]

factors(x, [i])?

[lkdvos]

But presumably better to change that in a separate PR, this is currently non-breaking

[mtfishman]

Agreed that kind of change should be done in a separate PR. factors(x, [i]) is probably a better interface.

[mtfishman]

I think in practice these aren't used so maybe they can be removed for now, I think I was undecided whether I would use arguments vs. arg1/arg2 but in practice I found arg1/arg2 were easier to use. That could be a separate PR of course, just bringing it up here since this PR reminded me of this issue.

[mtfishman]

Your suggestion of factors(x, [i]) would address this issue.

[mtfishman]

Looks good to me, ready to merge?

I see fillarrays.jl should also be made abstract but that is a deeper change that would require changing that code to catch cases where there is a Kronecker product involving Zeros in the broadcast style (which would be made easier by JuliaArrays/FillArrays.jl#385, but I think we could do something ourselves in the meantime). Those definitions are useful for us since the blockwise broadcasting/mapping code in BlockSparseArrays.jl creates Zeros blocks when there are unstored blocks, and those definitions help with removing those from the broadcast expression.