Checking the symmetry of a matrix

General Usage
, ,
1

I am trying to check the symmetry of a matrix which I was expecting to be symmetric.

I have printed A^T -A to check the deviation from symmetry. It is actually populated with very small numbers, as shown below.

 0.0          0.0           0.0           0.0           0.0           1.13687e-13   9.09495e-13   5.68434e-14   0.0          0.0           1.13687e-13   1.42109e-14
  0.0          0.0           0.0           0.0           0.0          -9.09495e-13   0.0           0.0           0.0          0.0          -9.09495e-13   0.0
  0.0          0.0           0.0           0.0           0.0          -9.09495e-13   2.84217e-14   0.0           0.0          5.68434e-14   0.0           0.0
  0.0          0.0           0.0           0.0           0.0           0.0          -7.10543e-15   2.84217e-14   0.0          0.0          -2.77556e-17   1.42109e-14
  0.0          0.0           0.0           0.0           0.0           0.0           0.0           0.0           0.0          2.22045e-16   0.0           3.55271e-15
 -1.13687e-13  9.09495e-13   9.09495e-13   0.0           0.0           0.0           0.0           0.0           0.0          0.0           3.55271e-15  -1.81899e-12
 -9.09495e-13  0.0          -2.84217e-14   7.10543e-15   0.0           0.0           0.0           0.0           0.0          0.0           0.0           0.0
 -5.68434e-14  0.0           0.0          -2.84217e-14   0.0           0.0           0.0           0.0           0.0          0.0           2.84217e-14   0.0
  0.0          0.0           0.0           0.0           0.0           0.0           0.0           0.0           0.0          1.13687e-13   0.0           0.0
  0.0          0.0          -5.68434e-14   0.0          -2.22045e-16   0.0           0.0           0.0          -1.13687e-13  0.0          -4.44089e-16   0.0
 -1.13687e-13  9.09495e-13   0.0           2.77556e-17   0.0          -3.55271e-15   0.0          -2.84217e-14   0.0          4.44089e-16   0.0           0.0
 -1.42109e-14  0.0           0.0          -1.42109e-14  -3.55271e-15   1.81899e-12   0.0           0.0           0.0          0.0           0.0           0.0

Please let me know if there is a way to give tolerance on checking the equality of terms about the main diagonal.

Thanks.

2

Are you aware of LinearAlgebra.issymmetric()?

afaik issymmetric() doesn’t have a way of specifying a tolerance though.

There is the build-in function isapprox, which allows you to specify absolute and relative tolerances. You could just broadcast this over the full array (A - A') and reduce with all, e.g you could define:.

isapproxsymmetric(A) = all(isapprox.(A - A', 0; rtol=1e-10))
3

Don’t broadcast it! isapprox works directly for matrices. Just do isapprox(A, A', rtol=1e-12) or whatever your tolerance is.

The reason you don’t want to broadcast isapprox elementwise is that when you are asking “is this difference small” you need to know “small compared to what?” And the correct scale to compare to is generally based on the whole matrix.

For example, this matrix A is very non-symmetrical, and in fact is defective:

julia> A = [1e-15 1e-15; 0 1e-15]
2×2 Matrix{Float64}:
 1.0e-15  1.0e-15
 0.0      1.0e-15

whereas this matrix B is nearly Hermitian:

julia> B = [1 1e-15; 0 1]
2×2 Matrix{Float64}:
 1.0  1.0e-15
 0.0  1.0

even though A - A^T = B - B^T!

julia> A - A' == B - B'
true

If you use isapprox on the whole matrix, it will tell you this:

julia> isapprox(A, A', rtol=1e-12)
false

julia> isapprox(B, B', rtol=1e-12)
true

whereas if you try to use isapprox element-wise it will fail for B:

julia> all(isapprox.(B, B', rtol=1e-12))
false

The reason is that, if you look elementwise, the test 0 ≈ 1e-15 must necessarily fail with any relative tolerance, because you don’t have any yardstick to know that 1e-15 is small in an “absolute” sense. See also this discussion and the isapprox documentation on comparison to zero.

9 Likes
4

Nice! Totally makes sense!
I wasn’t aware isapprox is defined for Matrix.

5

It’s defined for any type with - and norm (essentially, the concept of a “relative error” makes sense for any normed vector space). My mistake, we currently define it only for AbstractArray types and Numbers (and a few other odds and ends like Missing).

But we could define for any type with - and norm, since that’s all we really need. (I guess we also need eltype in order to choose the default rtol, or maybe eps.)

2 Likes
6

Also, as explained in links above, comparing to zero with a relative tolerance (rtol) is equivalent to ==, which is almost never what you want.

1 Like