[ML] Reduce variability in regression and classification results across our target platforms

docs/CHANGELOG.asciidoc

@@ -49,6 +49,8 @@
 * Improve robustness of anomaly detection to bad input data. (See {ml-pull}1114[#1114].)
 * Adds new `num_matches` and `preferred_to_categories` fields to category output.
   (See {ml-pull}1062[#1062].)
+* Reduce variability of classification and regression results across our target operating systems.
+  (See {ml-pull}1127[#1127].)
 * Switched data frame analytics model memory estimates from kilobytes to megabytes.
   (See {ml-pull}1126[#1126], issue: {issue}54506[#54506].)
 

include/core/CIEEE754.h

@@ -9,8 +9,8 @@
 
 #include <core/ImportExport.h>
 
-#include <stdint.h>
-#include <string.h>
+#include <cstdint>
+#include <cstring>
 
 namespace ml {
 namespace core {
@@ -40,14 +40,14 @@ class CORE_EXPORT CIEEE754 {
     //! as an integer.
     //! \note The actual "exponent" is "exponent - 1022" in two's complement.
     struct SDoubleRep {
-#ifdef __sparc                    // Add any other big endian architectures
-        uint64_t s_Sign : 1;      // sign bit
-        uint64_t s_Exponent : 11; // exponent
-        uint64_t s_Mantissa : 52; // mantissa
+#ifdef __sparc                         // Add any other big endian architectures
+        std::uint64_t s_Sign : 1;      // sign bit
+        std::uint64_t s_Exponent : 11; // exponent
+        std::uint64_t s_Mantissa : 52; // mantissa
 #else
-        uint64_t s_Mantissa : 52; // mantissa
-        uint64_t s_Exponent : 11; // exponent
-        uint64_t s_Sign : 1;      // sign bit
+        std::uint64_t s_Mantissa : 52; // mantissa
+        std::uint64_t s_Exponent : 11; // exponent
+        std::uint64_t s_Sign : 1;      // sign bit
 #endif
     };
 
@@ -57,15 +57,19 @@ class CORE_EXPORT CIEEE754 {
     //!
     //! \note This is closely related to std::frexp for double but returns
     //! the mantissa interpreted as an integer.
-    static void decompose(double value, uint64_t& mantissa, int& exponent) {
+    static void decompose(double value, std::uint64_t& mantissa, int& exponent) {
         SDoubleRep parsed;
         static_assert(sizeof(double) == sizeof(SDoubleRep),
                       "SDoubleRep definition unsuitable for memcpy to double");
         // Use memcpy() rather than union to adhere to strict aliasing rules
-        ::memcpy(&parsed, &value, sizeof(double));
+        std::memcpy(&parsed, &value, sizeof(double));
         exponent = static_cast<int>(parsed.s_Exponent) - 1022;
         mantissa = parsed.s_Mantissa;
     }
+
+    //! Drop \p bits trailing bits from the mantissa of \p value in a stable way
+    //! for different operating systems.
+    static double dropbits(double value, int bits);
 };
 }
 }

include/maths/CBoostedTreeLeafNodeStatistics.h

@@ -459,7 +459,8 @@ class MATHS_EXPORT CBoostedTreeLeafNodeStatistics final {
 
         bool operator<(const SSplitStatistics& rhs) const {
             return COrderings::lexicographical_compare(
-                s_Gain, s_Curvature, s_Feature, rhs.s_Gain, rhs.s_Curvature, rhs.s_Feature);
+                s_Gain, s_Curvature, s_Feature, s_SplitAt, // <- lhs
+                rhs.s_Gain, rhs.s_Curvature, rhs.s_Feature, rhs.s_SplitAt);
         }
 
         std::string print() const {

include/maths/CLbfgs.h

@@ -8,6 +8,7 @@
 #define INCLUDED_ml_maths_CLbfgs_h
 
 #include <maths/CLinearAlgebraShims.h>
+#include <maths/CTools.h>
 
 #include <boost/circular_buffer.hpp>
 
@@ -302,8 +303,9 @@ class CLbfgs {
         return m_BacktrackingMinDecrease * s * las::inner(m_Gx, m_P) / las::norm(m_P);
     }
 
-    constexpr double minimumStepSize() const {
-        return std::pow(m_StepScale, static_cast<double>(MAXIMUM_BACK_TRACKING_ITERATIONS));
+    double minimumStepSize() const {
+        return CTools::stable(std::pow(
+            m_StepScale, static_cast<double>(MAXIMUM_BACK_TRACKING_ITERATIONS)));
     }
 
 private:

include/maths/CTools.h

@@ -23,6 +23,7 @@
 
 #include <array>
 #include <cmath>
+#include <cstdint>
 #include <cstring>
 #include <iosfwd>
 #include <limits>
@@ -456,7 +457,7 @@ class MATHS_EXPORT CTools : private core::CNonInstantiatable {
             //      (interpreted as an integer) to the corresponding
             //      double value and fastLog uses the same approach
             //      to extract the mantissa.
-            uint64_t dx = 0x10000000000000ull / BINS;
+            std::uint64_t dx = 0x10000000000000ull / BINS;
             core::CIEEE754::SDoubleRep x;
             x.s_Sign = 0;
             x.s_Mantissa = (dx / 2) & core::CIEEE754::IEEE754_MANTISSA_MASK;
@@ -469,12 +470,12 @@ class MATHS_EXPORT CTools : private core::CNonInstantiatable {
                 // Use memcpy() rather than union to adhere to strict
                 // aliasing rules
                 std::memcpy(&value, &x, sizeof(double));
-                m_Table[i] = std::log2(value);
+                m_Table[i] = stable(std::log2(value));
             }
         }
 
         //! Lookup log2 for a given mantissa.
-        const double& operator[](uint64_t mantissa) const {
+        const double& operator[](std::uint64_t mantissa) const {
             return m_Table[mantissa >> FAST_LOG_SHIFT];
         }
 
@@ -494,7 +495,7 @@ class MATHS_EXPORT CTools : private core::CNonInstantiatable {
     //! \note This is taken from the approach given in
     //! http://www.icsi.berkeley.edu/pubs/techreports/TR-07-002.pdf
     static double fastLog(double x) {
-        uint64_t mantissa;
+        std::uint64_t mantissa;
         int log2;
         core::CIEEE754::decompose(x, mantissa, log2);
         return 0.693147180559945 * (FAST_LOG_TABLE[mantissa] + log2);
@@ -669,6 +670,15 @@ class MATHS_EXPORT CTools : private core::CNonInstantiatable {
     //! Compute \f$x^2\f$.
     static double pow2(double x) { return x * x; }
 
+    //! Compute a value from \p x which will be stable across platforms.
+    static double stable(double x) { return core::CIEEE754::dropbits(x, 1); }
+
+    //! A version of std::log which is stable across platforms.
+    static double stableLog(double x) { return stable(std::log(x)); }
+
+    //! A version of std::log which is stable across platforms.
+    static double stableExp(double x) { return stable(std::exp(x)); }
+
     //! Sigmoid function of \p p.
     static double sigmoid(double p) { return 1.0 / (1.0 + 1.0 / p); }
 
@@ -682,7 +692,7 @@ class MATHS_EXPORT CTools : private core::CNonInstantiatable {
     //! \param[in] sign Determines whether it's a step up or down.
     static double
     logisticFunction(double x, double width = 1.0, double x0 = 0.0, double sign = 1.0) {
-        return sigmoid(std::exp(std::copysign(1.0, sign) * (x - x0) / width));
+        return sigmoid(stableExp(std::copysign(1.0, sign) * (x - x0) / width));
     }
 
     //! Compute the softmax for the multinomial logit values \p logit.
@@ -696,7 +706,7 @@ class MATHS_EXPORT CTools : private core::CNonInstantiatable {
         double Z{0.0};
         double zmax{*std::max_element(z.begin(), z.end())};
         for (auto& zi : z) {
-            zi = std::exp(zi - zmax);
+            zi = stableExp(zi - zmax);
             Z += zi;
         }
         for (auto& zi : z) {

include/maths/CToolsDetail.h

@@ -316,7 +316,7 @@ void CTools::inplaceLogSoftmax(CDenseVector<SCALAR>& z) {
     double zmax{z.maxCoeff()};
     z.array() -= zmax;
     double Z{z.array().exp().sum()};
-    z.array() -= std::log(Z);
+    z.array() -= stableLog(Z);
 }
 }
 }

lib/api/unittest/CDataFrameAnalyzerFeatureImportanceTest.cc

@@ -485,8 +485,8 @@ BOOST_FIXTURE_TEST_CASE(testRegressionFeatureImportanceNoImportance, SFixture) {
             double c1{readShapValue(result, "c1")};
             double prediction{
                 result["row_results"]["results"]["ml"]["target_prediction"].GetDouble()};
-            // c1 explains 95% of the prediction value, i.e. the difference from the prediction is less than 2%.
-            BOOST_REQUIRE_CLOSE(c1, prediction, 5.0);
+            // c1 explains 94% of the prediction value, i.e. the difference from the prediction is less than 2%.
+            BOOST_REQUIRE_CLOSE(c1, prediction, 6.0);
             for (const auto& feature : {"c2", "c3", "c4"}) {
                 double c = readShapValue(result, feature);
                 BOOST_REQUIRE_SMALL(c, 2.0);

lib/core/CIEEE754.cc

@@ -7,14 +7,14 @@
 #include <core/CIEEE754.h>
 
 #include <cmath>
+#include <cstring>
 
 namespace ml {
 namespace core {
 
 double CIEEE754::round(double value, EPrecision precision) {
-    // This first decomposes the value into the mantissa
-    // and exponent to avoid the problem with overflow if
-    // the values are close to max double.
+    // First decomposes the value into the mantissa and exponent to avoid the
+    // problem with overflow if the values are close to max double.
 
     int exponent;
     double mantissa = std::frexp(value, &exponent);
@@ -39,5 +39,15 @@ double CIEEE754::round(double value, EPrecision precision) {
 
     return std::ldexp(mantissa, exponent);
 }
+
+double CIEEE754::dropbits(double value, int bits) {
+    SDoubleRep parsed;
+    static_assert(sizeof(double) == sizeof(SDoubleRep),
+                  "SDoubleRep definition unsuitable for memcpy to double");
+    std::memcpy(&parsed, &value, sizeof(double));
+    parsed.s_Mantissa &= ((IEEE754_MANTISSA_MASK << bits) & IEEE754_MANTISSA_MASK);
+    std::memcpy(&value, &parsed, sizeof(double));
+    return value;
+}
 }
 }