smosolver.h

Go to the documentation of this file.

/*
 * SPDX-FileCopyrightText: Copyright (c) 2019-2026, NVIDIA CORPORATION.
 * SPDX-License-Identifier: Apache-2.0
 */
 
#pragma once
 
#include <cuml/common/logger.hpp>
#include <cuml/svm/svm_model.h>
 
#include <raft/core/handle.hpp>
 
#include <rmm/device_scalar.hpp>
#include <rmm/device_uvector.hpp>
 
#include <thrust/device_ptr.h>
 
#include <cuvs/distance/distance.hpp>
#include <cuvs/distance/grammian.hpp>
 
#include <cassert>
#include <chrono>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <sstream>
#include <string>
#include <type_traits>
 
namespace ML {
namespace SVM {
 
template <typename math_t>
class SmoSolver {
 public:
  SmoSolver(const raft::handle_t& handle,
            SvmParameter param,
            cuvs::distance::kernels::KernelType kernel_type,
            cuvs::distance::kernels::GramMatrixBase<math_t>* kernel,
            bool is_precomputed = false)
    : handle(handle),
      C(param.C),
      tol(param.tol),
      kernel(kernel),
      kernel_type(kernel_type),
      is_precomputed(is_precomputed),
      cache_size(param.cache_size),
      nochange_steps(param.nochange_steps),
      epsilon(param.epsilon),
      svmType(param.svmType),
      stream(handle.get_stream()),
      return_buff(2, stream),
      alpha(0, stream),
      C_vec(0, stream),
      delta_alpha(0, stream),
      f(0, stream),
      y_label(0, stream)
  {
    ML::default_logger().set_level(param.verbosity);
  }
 
  void GetNonzeroDeltaAlpha(const math_t* vec,
                            int n_ws,
                            const int* idx,
                            math_t* nz_vec,
                            int* n_nz,
                            int* nz_idx,
                            cudaStream_t stream);
  template <typename MatrixViewType>
  void Solve(MatrixViewType matrix,
             int n_rows,
             int n_cols,
             math_t* y,
             const math_t* sample_weight,
             math_t** dual_coefs,
             int* n_support,
             SupportStorage<math_t>* support_matrix,
             int** idx,
             math_t* b,
             int max_iter       = -1,
             int max_outer_iter = -1,
             int max_inner_iter = 10000);
 
  void UpdateF(math_t* f, int n_rows, const math_t* delta_alpha, int n_ws, const math_t* cacheTile);
 
  void Initialize(math_t** y, const math_t* sample_weight, int n_rows, int n_cols);
 
  void InitPenalty(math_t* C_vec, const math_t* sample_weight, int n_rows);
 
  void SvcInit(const math_t* y);
 
  void SvrInit(const math_t* yr, int n_rows, math_t* yc, math_t* f);
 
  int GetNIter() { return n_iter; };
 
 private:
  const raft::handle_t& handle;
  cudaStream_t stream;
 
  int n_rows  = 0;  
  int n_cols  = 0;  
  int n_ws    = 0;  
  int n_train = 0;  
 
  // Buffers for the domain [n_train]
  rmm::device_uvector<math_t> alpha;    
  rmm::device_uvector<math_t> f;        
  rmm::device_uvector<math_t> y_label;  
 
  rmm::device_uvector<math_t> C_vec;  
 
  // Buffers for the working set [n_ws]
  rmm::device_uvector<math_t> delta_alpha;
 
  // Buffers to return some parameters from the kernel (iteration number, and
  // convergence information)
  rmm::device_uvector<math_t> return_buff;
  math_t host_return_buff[2];
 
  math_t C;
  math_t tol;      
  math_t epsilon;  
 
  cuvs::distance::kernels::GramMatrixBase<math_t>* kernel;
  cuvs::distance::kernels::KernelType kernel_type;
  bool is_precomputed;  
  float cache_size;     
 
  SvmType svmType;  
 
  // Variables to track convergence of training
  math_t diff_prev;
  int n_small_diff;
  int nochange_steps;
  int n_increased_diff;
  int n_outer_iter;
  int n_iter;
  bool report_increased_diff;
 
  bool CheckStoppingCondition(math_t diff)
  {
    if (diff > diff_prev * 1.5 && n_outer_iter > 0) {
      // Ideally, diff should decrease monotonically. In practice we can have
      // small fluctuations (10% increase is not uncommon). Here we consider a
      // 50% increase in the diff value large enough to indicate a problem.
      // The 50% value is an educated guess that triggers the convergence debug
      // message for problematic use cases while avoids false alarms in many
      // other cases.
      n_increased_diff++;
    }
    if (report_increased_diff && n_outer_iter > 100 && n_increased_diff > n_outer_iter * 0.1) {
      CUML_LOG_DEBUG(
        "Solver is not converging monotonically. This might be caused by "
        "insufficient normalization of the feature columns. In that case "
        "MinMaxScaler((0,1)) could help. Alternatively, for nonlinear kernels, "
        "you can try to increase the gamma parameter. To limit execution time, "
        "you can also adjust the number of iterations using the max_iter "
        "parameter.");
      report_increased_diff = false;
    }
    bool keep_going = true;
    if (abs(diff - diff_prev) < 0.001 * tol) {
      n_small_diff++;
    } else {
      diff_prev    = diff;
      n_small_diff = 0;
    }
    if (n_small_diff > nochange_steps) {
      CUML_LOG_ERROR(
        "SMO error: Stopping due to unchanged diff over %d"
        " consecutive steps",
        nochange_steps);
      keep_going = false;
    }
    if (diff < tol) keep_going = false;
    if (isnan(diff)) {
      std::string txt;
      if (std::is_same<float, math_t>::value) {
        txt +=
          " This might be caused by floating point overflow. In such case using"
          " fp64 could help. Alternatively, try gamma='scale' kernel"
          " parameter.";
      }
      THROW("SMO error: NaN found during fitting.%s", txt.c_str());
    }
    return keep_going;
  }
 
  int GetDefaultMaxIter(int n_train, int max_outer_iter)
  {
    if (max_outer_iter == -1) {
      max_outer_iter = n_train < std::numeric_limits<int>::max() / 100
                         ? n_train * 100
                         : std::numeric_limits<int>::max();
      max_outer_iter = max(100000, max_outer_iter);
    }
    // else we have user defined iteration count which we do not change
    return max_outer_iter;
  }
 
  void ResizeBuffers(int n_train, int n_cols)
  {
    // This needs to know n_train, therefore it can be only called during solve
    alpha.resize(n_train, stream);
    C_vec.resize(n_train, stream);
    f.resize(n_train, stream);
    delta_alpha.resize(n_ws, stream);
    if (svmType == EPSILON_SVR) y_label.resize(n_train, stream);
  }
 
  void ReleaseBuffers()
  {
    alpha.resize(0, stream);
    delta_alpha.resize(0, stream);
    f.resize(0, stream);
    y_label.resize(0, stream);
  }
};
 
};  // end namespace SVM
};  // end namespace ML