issymmetric with rtol, atol

[jamblejoe]

Hi,
would be great if one could pass arguments to issymmetric, which specifies the absolut and/or relative tolerance between the upper and lower triangle, similar to isapprox(x,y; atol=..., rtol=...).
Thanks!

[StefanKarpinski]

Would this check isapprox(X[i,j], X[j,i], atol=atol, rtol=rtol) for all indices or check that isapprox(X, transpose(X), atol=atol, rtol=rtol). I'm guessing the latter is more correct.

[StefanKarpinski]

cc @stevengj

[stevengj]

I think that isapprox(X, transpose(X), atol=atol, rtol=rtol) is definitely more reasonable here than an element-wise test. Probably should also pass through a norm argument.

(Using an approximate-equality test by default in things like Cholesky routines is a tricker question, which seems like it would require some careful analysis. See also the discussion in #182.)

[jlapeyre]

I wrote a package to try to approach this in a systematic way. You can have a default notion of closeness for a particular method. But, there is a uniform interface allowing they user to specify what they want.
IsApprox.jl

In particular, for issymmetric (as for many other functions), you can choose exact equality, a matrix norm, or element-wise equality (for a quick exit).

https://github.com/jlapeyre/IsApprox.jl/blob/11af339b500ef7ef37735d689f1d95a441360813/src/base_applications.jl#L53-L92

[stevengj]

What's the advantage of calling issymmetric(A, Approx(...)) over isapprox(A, A', ...)?

[oscardssmith]

It should be about 2x faster.

[jlapeyre]

What's the advantage of calling issymmetric(A, Approx(...)) over isapprox(A, A', ...)?

@stevengjl, In this case, the former reduces to the latter. The sequence of calls is

issymmetric(A, Approx())
ishermitian(A, Approx())
isapprox(Approx(), A, adjoint(A))
# and the final call above goes to the following method.
Base.isapprox(a::Union{EachApprox, Approx}, x, y) = isapprox(x, y; pairs(a.kw)...)

(isapprox takes AbstractApprox as the first argument so that it may be used as a replacement in Base without breaking code)

An advantage is that you can drop in EachApprox, instead and get element-wise comparison:

m = rand(1000,1000);
A = m + m' .* (1 + rand() * 1e-12);

julia> @btime IsApprox.issymmetric($A, EachApprox(atol=1e-15))
  8.555 ns (0 allocations: 0 bytes)
false

julia> @btime IsApprox.issymmetric($A, Approx(atol=1e-15))
  3.638 ms (2 allocations: 7.63 MiB)
false

(EachApprox is not a great name, maybe ElementWise is better.)
The code paths diverge above after the call to ishermitian.

The idea is to do this uniformly for everything. For example iszero behaves the same way (that is, as you would expect):

    m = [9.914830625828892e-11 4.473996666408231e-11 3.066184129948095e-11;
         9.914231937308404e-12 7.07270051593276e-11 6.487445728681865e-11;
         9.587109787642197e-11 6.520784660158512e-11 1.2903560723094866e-11]

    @test iszero(m, Approx(atol=2e-10)) # in norm
    @test ! iszero(m, Approx(atol=1e-10))
    @test iszero(m, EachApprox(atol=2e-10)) # element-wise
    @test iszero(m, EachApprox(atol=1e-10))
    @test ! iszero(m) # exact
    @test ! iszero(m, Equal()) # exact

What I really want to avoid is that each of the many predicates is treated separately to add support for closeness in norm and or elementwise, (and as I recently added, up to a global phase) in a different, and ad hoc, way. I wrote this package after seeing exactly this happen in a number of quantum information packages.

[fredrikekre]

Duplicate of #558.

(No need to discuss this in two places as is happening now).