diff --git a/include/ceres/jet.h b/include/ceres/jet.h
index e45caf2a2..38a69f481 100644
--- a/include/ceres/jet.h
+++ b/include/ceres/jet.h
@@ -632,28 +632,34 @@ inline Jet<T, N> hypot(const Jet<T, N>& x, const Jet<T, N>& y) {
 
 template <typename T, int N>
 inline Jet<T, N> fmax(const Jet<T, N>& x, const Jet<T, N>& y) {
-  return x < y ? y : x;
+  using std::isgreater;
+  return isnan(y.a) || isgreater(x.a, y.a) ? x : y;
 }
+
 template <typename T, int N>
 inline Jet<T, N> fmax(const Jet<T, N>& x, const T& y) {
-  return x < y ? Jet<T, N>{y} : x;
+  return fmax(x, Jet<T, N>{y});
 }
+
 template <typename T, int N>
 inline Jet<T, N> fmax(const T& x, const Jet<T, N>& y) {
-  return x < y ? y : Jet<T, N>{x};
+  return fmax(Jet<T, N>{x}, y);
 }
 
 template <typename T, int N>
 inline Jet<T, N> fmin(const Jet<T, N>& x, const Jet<T, N>& y) {
-  return y < x ? y : x;
+  using std::isless;
+  return isnan(x.a) || isless(y.a, x.a) ? y : x;
 }
+
 template <typename T, int N>
 inline Jet<T, N> fmin(const Jet<T, N>& x, const T& y) {
-  return y < x ? Jet<T, N>{y} : x;
+  return fmin(x, Jet<T, N>{y});
 }
+
 template <typename T, int N>
 inline Jet<T, N> fmin(const T& x, const Jet<T, N>& y) {
-  return y < x ? y : Jet<T, N>{x};
+  return fmin(Jet<T, N>{x}, y);
 }
 
 // erf is defined as an integral that cannot be expressed analytically
diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc
index bceb9e480..32b136175 100644
--- a/internal/ceres/jet_test.cc
+++ b/internal/ceres/jet_test.cc
@@ -834,6 +834,16 @@ TEST(Jet, Jet) {
     VL << "z = " << z;
     ExpectJetsClose(x, z);
   }
+  {
+    J z = fmax(std::numeric_limits<double>::quiet_NaN(), x);
+    VL << "z = " << z;
+    ExpectJetsClose(x, z);
+  }
+  {
+    J z = fmax(x, std::numeric_limits<double>::quiet_NaN());
+    VL << "z = " << z;
+    ExpectJetsClose(x, z);
+  }
 
   {
     J z = fmin(x, y);
@@ -865,6 +875,17 @@ TEST(Jet, Jet) {
     VL << "z = " << z;
     ExpectJetsClose(J{y.a}, z);
   }
+  {
+    J z = fmin(x, std::numeric_limits<double>::quiet_NaN());
+    VL << "z = " << z;
+    ExpectJetsClose(x, z);
+  }
+  {
+    J z = fmin(std::numeric_limits<double>::quiet_NaN(), x);
+    VL << "z = " << z;
+    ExpectJetsClose(x, z);
+  }
+
   {  // copysign(x, +1)
     J z = copysign(x, J{+1});
     VL << "z = " << z;