diff --git a/base/floatfuncs.jl b/base/floatfuncs.jl
index f34f6875c7b6e..adb7bb20d1025 100644
--- a/base/floatfuncs.jl
+++ b/base/floatfuncs.jl
@@ -170,10 +170,11 @@ end
 
 # isapprox: approximate equality of numbers
 """
-    isapprox(x, y; rtol::Real=sqrt(eps), atol::Real=0)
+    isapprox(x, y; rtol::Real=sqrt(eps), atol::Real=0, nans::Bool=false)
 
 Inexact equality comparison: `true` if `norm(x-y) <= atol + rtol*max(norm(x), norm(y))`. The
-default `atol` is zero and the default `rtol` depends on the types of `x` and `y`.
+default `atol` is zero and the default `rtol` depends on the types of `x` and `y`. The keyword
+argument `nans` determines whether or not NaN values are considered equal (defaults to false).
 
 For real or complex floating-point values, `rtol` defaults to
 `sqrt(eps(typeof(real(x-y))))`. This corresponds to requiring equality of about half of the
@@ -188,8 +189,8 @@ approximately equal component-wise.
 The binary operator `≈` is equivalent to `isapprox` with the default arguments, and `x ≉ y`
 is equivalent to `!isapprox(x,y)`.
 """
-function isapprox(x::Number, y::Number; rtol::Real=rtoldefault(x,y), atol::Real=0)
-    x == y || (isfinite(x) && isfinite(y) && abs(x-y) <= atol + rtol*max(abs(x), abs(y)))
+function isapprox(x::Number, y::Number; rtol::Real=rtoldefault(x,y), atol::Real=0, nans::Bool=false)
+    x == y || (isfinite(x) && isfinite(y) && abs(x-y) <= atol + rtol*max(abs(x), abs(y))) || (nans && isnan(x) && isnan(y))
 end
 
 const ≈ = isapprox
diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl
index dff7113069929..fe49a356ab7bd 100644
--- a/base/linalg/generic.jl
+++ b/base/linalg/generic.jl
@@ -1216,13 +1216,13 @@ promote_leaf_eltypes(x::Union{AbstractArray,Tuple}) = mapreduce(promote_leaf_elt
 # `a ≈ a` is `true`.
 function isapprox(x::AbstractArray, y::AbstractArray;
     rtol::Real=Base.rtoldefault(promote_leaf_eltypes(x),promote_leaf_eltypes(y)),
-    atol::Real=0, norm::Function=vecnorm)
+    atol::Real=0, nans::Bool=false, norm::Function=vecnorm)
     d = norm(x - y)
     if isfinite(d)
         return d <= atol + rtol*max(norm(x), norm(y))
     else
         # Fall back to a component-wise approximate comparison
-        return all(ab -> isapprox(ab[1], ab[2]; rtol=rtol, atol=atol), zip(x, y))
+        return all(ab -> isapprox(ab[1], ab[2]; rtol=rtol, atol=atol, nans=nans), zip(x, y))
     end
 end
 
diff --git a/base/linalg/uniformscaling.jl b/base/linalg/uniformscaling.jl
index 546b4fd6ae5bd..583ec9435a6f5 100644
--- a/base/linalg/uniformscaling.jl
+++ b/base/linalg/uniformscaling.jl
@@ -165,8 +165,10 @@ broadcast(::typeof(/), J::UniformScaling,x::Number) = UniformScaling(J.λ/x)
 
 ==(J1::UniformScaling,J2::UniformScaling) = (J1.λ == J2.λ)
 
-isapprox{T<:Number,S<:Number}(J1::UniformScaling{T}, J2::UniformScaling{S};
-                              rtol::Real=Base.rtoldefault(T,S), atol::Real=0) = isapprox(J1.λ, J2.λ, rtol=rtol, atol=atol)
+function isapprox{T<:Number,S<:Number}(J1::UniformScaling{T}, J2::UniformScaling{S};
+                              rtol::Real=Base.rtoldefault(T,S), atol::Real=0, nans::Bool=false)
+    isapprox(J1.λ, J2.λ, rtol=rtol, atol=atol, nans=nans)
+end
 
 function copy!(A::AbstractMatrix, J::UniformScaling)
     size(A,1)==size(A,2) || throw(DimensionMismatch("a UniformScaling can only be copied to a square matrix"))
diff --git a/test/floatapprox.jl b/test/floatapprox.jl
index fe76153082173..ab432de6f71d8 100644
--- a/test/floatapprox.jl
+++ b/test/floatapprox.jl
@@ -59,3 +59,13 @@
 @test [0,1] ≈ [1e-9, 1]
 @test [0,Inf] ≈ [0,Inf]
 @test [0,Inf] ≉ [0,-Inf]
+
+# issue #19936
+for elty in (Float16,Float32,Float64)
+    nan  = elty(NaN)
+    half = elty(0.5)
+    @test !isapprox(nan, nan)
+    @test isapprox(nan, nan, nans=true)
+    @test !isapprox([half, nan, half], [half, nan, half])
+    @test isapprox([half, nan, half], [half, nan, half], nans=true)
+end