What causes `AssertionError` when asserting a diagonalized vector in Julia?

Viewed 261

I've been doing some exercises in julia, and I am currently trying to @assert a vector that has been diagonalized into a matrix against a "solution matrix" given in the exercise notebook. However, I get an AssertionError when asserting my code against the provided solution. Example of my code:

julia> using LinearAlgebra

julia> A =
[
 140   97   74  168  131
  97  106   89  131   36
  74   89  152  144   71
 168  131  144   54  142
 131   36   71  142   36
]

5×5 Matrix{Int64}:
 140   97   74  168  131
  97  106   89  131   36
  74   89  152  144   71
 168  131  144   54  142
 131   36   71  142   36

julia> A_eigv = eigen(A).values

5-element Vector{Float64}:
 -128.49322764802145
  -55.887784553057
   42.752167279318854
   87.16111477514494
  542.4677301466137

julia> A_diag = Diagonal(A_eigv)

5×5 Diagonal{Float64, Vector{Float64}}:
 -128.493     ⋅        ⋅        ⋅         ⋅ 
     ⋅     -55.8878    ⋅        ⋅         ⋅ 
     ⋅        ⋅      42.7522    ⋅         ⋅ 
     ⋅        ⋅        ⋅      87.1611     ⋅ 
     ⋅        ⋅        ⋅        ⋅      542.468

julia> @assert A_diag ==  [-128.493    0.0      0.0      0.0       0.0;
    0.0    -55.8878   0.0      0.0       0.0;
    0.0      0.0     42.7522   0.0       0.0;
    0.0      0.0      0.0     87.1611    0.0;
    0.0 0.0      0.0      0.0     542.468]

AssertionError: A_diag == [-128.493 0.0 0.0 0.0 0.0; 0.0 -55.8878 0.0 0.0 0.0; 0.0 0.0 42.7522 0.0 0.0; 0.0 0.0 0.0 87.1611 0.0; 0.0 0.0 0.0 0.0 542.468]

Stacktrace:
 [1] top-level scope
   @ In[90]:1
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1094

My inital assumption was that a difference in number of decimals was the cause of the error. I therefore replaced == with ≈ (\approx). However, as the code example below shows, the error persists:

julia> @assert A_diag ≈ #\approx
[-128.493    0.0      0.0      0.0       0.0;
    0.0    -55.8878   0.0      0.0       0.0;
    0.0      0.0     42.7522   0.0       0.0;
    0.0      0.0      0.0     87.1611    0.0;
    0.0 0.0      0.0      0.0     542.468]

AssertionError: A_diag ≈ [-128.493 0.0 0.0 0.0 0.0; 0.0 -55.8878 0.0 0.0 0.0; 0.0 0.0 42.7522 0.0 0.0; 0.0 0.0 0.0 87.1611 0.0; 0.0 0.0 0.0 0.0 542.468]

Stacktrace:
 [1] top-level scope
   @ In[97]:1
 [2] eval
   @ ./boot.jl:360 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1094

I've been reading through my code multiple times now, and I am at a loss. The values in my diagonal matrix (A_diag) are seemingly identical to the solution matrix. Furthermore, setting the statement to approximately equal (\approx) renders the same error, so I assume I can count out decimal error.

My main question is: What causes the AssertionError?

1 Answers

No, the dots are treated as 0 for the purpose of testing equality.

julia> Diagonal(1:2) == [1 0; 0 2]
true

Your problem is that your arrays are actually not equal; -128.49322764802145 is not the same as -128.493. (The pretty-printed version of the array truncates floats for display, but that's not the true underlying value!).

[Edit:]

Using ≈ (\approx) will also fail in this case. The reason for this is explained in the documentation for isapprox()

The binary operator ≈ is equivalent to isapprox with the default arguments,

if an atol > 0 is not specified, rtol defaults to the square root of eps of the type of x or y, whichever is bigger (least precise).

Essentially, what this means is that ≈ will test for approximate equality with a relative tolerance of √eps() which is approximately equal to 1.5e-8, or 0.0000015%. This tolerance is waaay too low, and increasing the tolerance will resolve the issue. E.g.:

# Option 1: Absolute tolerance. Set to a reasonable max deviation:
julia> isapprox(A_diag, sol_mat, atol = 1e-3)
true

# Option 2: Relative tolerance. Setting rtol = 1e-n, n is the number of significant digits in either matrices will work in most cases.
julia> isapprox(A_diag, sol_mat, rtol = 1e-6)
true

Since the solution matrix provides the values in six significant digits, another alternative is to round the values in A_diag to this number og digits and test for equality. E.g.:

julia> round.(A_diag, RoundNearestTiesUp, sigdigits=6) == 
[-128.493    0.0      0.0      0.0       0.0;
    0.0    -55.8878   0.0      0.0       0.0;
    0.0      0.0     42.7522   0.0       0.0;
    0.0      0.0      0.0     87.1611    0.0;
    0.0      0.0      0.0      0.0     542.468]
true
Related