abs(Dual(-0.0,1.0))

[goretkin]

Perhaps we should have abs(Dual(-0.0,1.0)) == Dual(0.0,-1.0)?

julia> abs(Dual(0.0,1.0))
0.0 + 1.0du

julia> abs(Dual(-0.0,1.0))
-0.0 + 1.0du

julia> abs(Dual(-eps(),1.0))
2.220446049250313e-16 - 1.0du

[mlubin]

The current behavior is technically correct but this seems reasonable also. Is there an efficient way to check this?

[goretkin]

I don't know enough about -0.0. Could you explain what makes the current behavior technically correct?

[mlubin]

Absolute value isn't differentiable at zero, so all we can expect is a valid subgradient. Both 1.0 and -1.0 are valid subgradients at that point, because they're both tangent to the graph of absolute value.

[goretkin]

What about the sign on the real part? On Dec 9, 2014 5:59 PM, "Miles Lubin" notifications@github.com wrote:

Absolute value isn't differentiable at zero, so all we can expect is a valid subgradient. Both 1.0 and -1.0 are valid subgradients at that point, because they're both tangent to the graph of absolute value.

— Reply to this email directly or view it on GitHub https://github.com/JuliaDiff/DualNumbers.jl/issues/14#issuecomment-66374791 .

[mlubin]

-0.0 and 0.0 is a weird floating-point curiosity:

julia> 0.0 == -0.0
true
julia> -0.0 >= 0.0
true

[goretkin]

Yeah, I'm with you... On Dec 9, 2014 6:13 PM, "Miles Lubin" notifications@github.com wrote:

-0.0 and 0.0 is a weird floating-point curiosity:

julia> 0.0 == -0.0 true julia> -0.0 >= 0.0 true

— Reply to this email directly or view it on GitHub https://github.com/JuliaDiff/DualNumbers.jl/issues/14#issuecomment-66376648 .

[goretkin]

I was thinking about this again.

Currently we have this situation: julia> z = Dual(-0.0,1.0) -0.0 + 1.0du

julia> abs(z) -0.0 + 1.0du

julia> sqrt(abs2(z)) dual(0.0,NaN)

I understand why Calculus.jl defines the derivatives of square root as it does, however I think it is possible to have a definition for sqrt(::Dual) that does work as expected in the above case.

#https://en.wikipedia.org/wiki/Dual_number#Linear_representation
function mat(z::Dual)
       r = eye(2)*z.re
       r[1,2] = z.du
       return r
end

julia> mat(z)
2x2 Array{Float64,2}:
 -0.0   1.0
 -0.0  -0.0

julia> sqrtm(mat(abs2(z)))
2x2 Array{Complex{Float64},2}:
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im

the final matrix corresponds to Dual(0.0,0.0)

Note that the current definition and the sqrtm(mat(::Dual)) operation above agree in a lot of weird cases already:

julia> sqrtm(mat(Dual(0.0,1.0)) )
2x2 Array{Float64,2}:
 0.0  Inf  
 0.0    0.0

julia> sqrt(Dual(0.0,1.0))
dual(0.0,Inf)

julia> sqrtm(mat(Dual(-0.0,1.0)) )
2x2 Array{Float64,2}:
 -0.0  -Inf  
  0.0    -0.0

julia> sqrt(Dual(-0.0,1.0))
dual(-0.0,-Inf)

just not in the following case

julia> sqrtm(mat(Dual(0.0,-0.0)))
2x2 Array{Complex{Float64},2}:
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im

julia> sqrt(Dual(0.0, -0.0))
dual(0.0,NaN)

[dlfivefifty]

How is sqrtm defined at [0 1; 0 0]? The usual definitions of matrix functions (Taylor series, Cauchy integral formula) don't seem to apply when the spectrum touches a singularity. Or is it entry wise anyaltic continuation?

[goretkin]

@dlfivefifty I don't know, but I can say that it is doing what matches my intuition based on the auto-diff properties of Dual numbers

#repeated from above
julia> sqrtm(mat(Dual(0.0,1.0)) )
2x2 Array{Float64,2}:
 0.0  Inf  
 0.0    0.0

julia> a = [0 1; 0 0]
2x2 Array{Int64,2}:
 0  1
 0  0

julia> sqrtm(a+eye(2)*nextfloat(0.0))
2x2 Array{Float64,2}:
 2.22276e-162  2.24946e161 
 0.0           2.22276e-162

julia> sqrt(Dual(nextfloat(0.0), 1.0))
2.2227587494850775e-162 + 2.2494568972715982e161du

[dlfivefifty]

Ok I guess that's the analytic continuation version.