Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix computation of norms of AbstractArray{<:Measurement} #120

Merged
merged 1 commit into from
May 19, 2022

Conversation

giordano
Copy link
Member

In Julia v1.9 we currently have

julia> A = [(14 ± 0.1) (23 ± 0.2); (-12 ± 0.3) (24 ± 0.4)];

julia> norm(A * inv(A) - I, Inf)
0.0 ± 3.5e-19

julia> norm(A * inv(A) - I)
NaN ± NaN

This is due to the fact that the Inf-norm has now zero value with non-zero
uncertainty, but some checks in the computation of the p-norms assume that if a
number is non-zero then dividing by it wouldn't result in an overflow (or NaN).
However this is not true for numbers of the type 0 ± u with u != 0, which
aren't the additive identity element but nevertheless would result in
overflow (or NaN) when dividing by them.

While a proper solution would be to have in the language two distinct traits for
"additive identity element" (which is already represented by iszero) and
"non-overflowing denominator" (which we currently miss, but I also don't know
any language implementing this trait), Daniel Wennberg pointed out that the
Inf-norm of an array should always have zero-uncertainty if the value is also
zero, thus resulting in a proper "additive identity element" zero which would
satisfy the assumptions in the computations of the other p-norms. This change
defines the method LinearAlgebra.normInf(::AbstractArray{<:Measurement}) to
always force zero uncertainty in the Inf-norm when the value is also zero.


Thanks a lot @danielwe for the help and the in-depth discussion about how to address the problem.

@codecov-commenter
Copy link

codecov-commenter commented May 18, 2022

Codecov Report

Merging #120 (f8cf429) into master (6bd561c) will decrease coverage by 2.68%.
The diff coverage is 100.00%.

❗ Current head f8cf429 differs from pull request most recent head d2ada3d. Consider uploading reports for the commit d2ada3d to get more accurate results

@@            Coverage Diff             @@
##           master     #120      +/-   ##
==========================================
- Coverage   95.54%   92.85%   -2.69%     
==========================================
  Files          12       13       +1     
  Lines         741      546     -195     
==========================================
- Hits          708      507     -201     
- Misses         33       39       +6     
Impacted Files Coverage Δ
src/Measurements.jl 85.18% <ø> (-2.70%) ⬇️
src/linear_algebra.jl 100.00% <100.00%> (ø)
src/comparisons-tests.jl 70.00% <0.00%> (-10.00%) ⬇️
src/conversions.jl 90.47% <0.00%> (-9.53%) ⬇️
src/show.jl 76.47% <0.00%> (-6.87%) ⬇️
src/derivatives-type.jl 94.73% <0.00%> (-5.27%) ⬇️
src/math.jl 94.64% <0.00%> (-3.29%) ⬇️
src/special-functions.jl 95.58% <0.00%> (-0.61%) ⬇️
src/utils.jl 100.00% <0.00%> (ø)
... and 4 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 6bd561c...d2ada3d. Read the comment docs.

In Julia v1.9 we currently have

```julia
julia> A = [(14 ± 0.1) (23 ± 0.2); (-12 ± 0.3) (24 ± 0.4)];

julia> norm(A * inv(A) - I, Inf)
0.0 ± 3.5e-19

julia> norm(A * inv(A) - I)
NaN ± NaN
```

This is due to the fact that the Inf-norm has now zero value with non-zero
uncertainty, but some checks in the computation of the p-norms assume that if a
number is non-zero then dividing by it wouldn't result in an overflow (or NaN).
However this is not true for numbers of the type `0 ± u` with `u != 0`, which
aren't the additive identity element but nevertheless would result in
overflow (or NaN) when dividing by them.

While a proper solution would be to have in the language two distinct traits for
"additive identity element" (which is already represented by `iszero`) and
"non-overflowing denominator" (which we currently miss, but I also don't know
any language implementing this trait), Daniel Wennberg pointed out that the
Inf-norm of an array should always have zero-uncertainty if the value is also
zero, thus resulting in a proper "additive identity element" zero which would
satisfy the assumptions in the computations of the other p-norms.  This change
defines the method `LinearAlgebra.normInf(::AbstractArray{<:Measurement})` to
always force zero uncertainty in the Inf-norm when the value is also zero.
@giordano giordano merged commit a2060ea into master May 19, 2022
@giordano giordano deleted the mg/linear-algebra branch May 19, 2022 00:17
@danielwe
Copy link

danielwe commented May 19, 2022

Let me just add here that there may be an order of limits issue in the argument. My understanding is that linear uncertainty propagation is consistent in the limit that the uncertainty is small compared to the scale over which the function you propagate through changes slope. In the argument, we've taken this limit before taking the limit p -> Inf. I haven't thought about the consequence of reordering these limits, I suspect you'd end up picking out the zero with the largest uncertainty, which I guess is a pretty logical way to define the largest absolute value of two uncertain zeros.

Then again, the next logical step would be to apply the same reasoning to function composition, which would mean that sqrt((0.0 ± 1.0)^2) == 0.0 ± 0.0 is wrong since the composed function is actually equivalent to abs and abs(0.0 ± 1.0) = 0.0 ± 1.0. Consistent treatment of composition along these lines would probably result in finite-p norms of [0.0 ± 1.0, 0.0 ± 0.5] having finite uncertainty too. However, I don't know what the strategy would be for always deciding on the correct value of 0 * Inf in every case. The current strategy is that once uncertainty has vanished it never comes back, which seems like a very pragmatic and reasonable decision.

In other words, as long as this package is happy with (abs(0.0 ± 1.0)^p)^(1/p) == 0.0 ± 0.0 for finite p I think it makes good sense to extend that to p -> Inf.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants