Better definition of isapprox

[mtfishman]

Followup to the comments on isapprox in #54.

[lkdvos]

copying this here:

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

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 81.73%. Comparing base (9f6111a) to head (3af6e96).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #55      +/-   ##
==========================================
+ Coverage   81.21%   81.73%   +0.52%     
==========================================
  Files           9        9              
  Lines         628      646      +18     
==========================================
+ Hits          510      528      +18     
  Misses        118      118              

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 point, I think I missed a few cross terms (I was thinking they would be subleading but I don't think they are). I was trying to circumvent that issue of adding and subtracting large values by first rewriting it in terms of differences a1 - b1 and a2 - b2.

[lkdvos]

Why not do a similar approach to other parts of the code, and check if arg1 is equal, if so compare arg2, and vice versa, and if nothing else actually instantiate the difference? I think here the precision actually really does matter, and taking a square root looks dangerous to me...

[mtfishman]

I added in cross terms but also dropped some terms that are second order in the differences a1 - b1 and a2 - b2, with this new version I see:

julia> using StableRNGs

julia> rng = StableRNG(123);

julia> a = randn(rng, 100, 100) ⊗ randn(rng, 100, 100);

julia> b = (arg1(a) + 1e-1 * randn(rng, size(arg1(a)))) ⊗ (arg2(a) + 1e-1 * randn(rng, size(arg2(a))));

julia> err = norm(collect(a) - collect(b)) / max(norm(a), norm(b))
0.14243654495933608

julia> rtol = 0.143; isapprox(a, b; rtol), isapprox(collect(a), collect(b); rtol)
(false, true)

julia> rtol = 0.1435; isapprox(a, b; rtol), isapprox(collect(a), collect(b); rtol)
(true, true)

so it looks like the approximation is pretty good.

[mtfishman]

Why not do a similar approach to other parts of the code, and check if arg1 is equal, if so compare arg2, and vice versa, and if nothing else actually instantiate the difference? I think here the precision actually really does matter, and taking a square root looks dangerous to me...

We can definitely do that, some tests rely on the current more general behavior that allows both factors to be different from each other but maybe for those we can just materialize the Kronecker product, I was curious to see if something like this could work that just depends on differences between factors.

[mtfishman]

Why not do a similar approach to other parts of the code, and check if arg1 is equal, if so compare arg2, and vice versa, and if nothing else actually instantiate the difference?

I worked on it a bit more and wrote a more accurate distance function for Kronecker arrays that includes cross terms and second order terms, but even so in the end I agree with you that we should analyze the precision more carefully before using that.

I've included that distance function in the package with a test in case we want to use it at some point, but for now I made isapprox stricter where it only works if the first or second arguments exactly match (in general I've been stricter with materializing where I don't materialize automatically to force users to be aware of when materializing occurs, but we can loosen that if that seems to strict).

I rewrote tests that were relying on the less strict behavior of isapprox from before by manually calling isapprox on each term, which I think is more explicit anyway.