Introduce AbstractKroneckerArray
#54
Conversation
Codecov Report❌ Patch coverage is
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
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
| # TODO: this definition doesn't fully retain the original meaning: | ||
| # ‖a - b‖ < atol could be true even if the following check isn't |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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)))
endwhich would work for SectorArrays since a1 == b1.
There was a problem hiding this comment.
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:
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
| @doc """ | ||
| arg1(AB::AbstractKroneckerArray{T, N}) | ||
| Extract the first factor (`A`) of the Kronecker array `AB = A ⊗ B`. | ||
| """ arg1 | ||
|
|
||
| @doc """ | ||
| arg2(AB::AbstractKroneckerArray{T, N}) | ||
| Extract the second factor (`B`) of the Kronecker array `AB = A ⊗ B`. | ||
| """ arg2 |
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
But presumably better to change that in a separate PR, this is currently non-breaking
There was a problem hiding this comment.
Agreed that kind of change should be done in a separate PR. factors(x, [i]) is probably a better interface.
| arguments(a::AbstractKroneckerArray) = (arg1(a), arg2(a)) | ||
| arguments(a::AbstractKroneckerArray, n::Int) = arguments(a)[n] | ||
| argument_types(a::AbstractKroneckerArray) = argument_types(typeof(a)) |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
Your suggestion of factors(x, [i]) would address this issue.
|
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 |
What it says in the title :)