3
\$\begingroup\$

The most famous library for Support Vector Machine (SVM) algorithm is libsvm (https://github.com/cjlin1/libsvm/), but I felt that its code style is too old, I rewrote in newer C++ as a hobby project.

Link : https://github.com/frozenca/ML_practice/tree/main/ML/SVM

SVMData.h

#ifndef FROZENCA_SVMDATA_H
#define FROZENCA_SVMDATA_H

#include <array>
#include <charconv>
#include <ranges>
#include <iostream>
#include <string>
#include <string_view>
#include <variant>
#include <vector>

namespace frozenca {

using Samples = std::vector<std::vector<float>>;
using SubSamples = std::vector<std::reference_wrapper<const std::vector<float>>>;

class SVMTrainData {
private:
    std::variant<Samples, SubSamples> X_;
public:
    std::vector<float> y_;
    std::size_t n_;
    std::size_t C_;
    std::vector<std::pair<float, float>> feature_bounds_;
    friend class SVMModel;

    // full data construct
    explicit SVMTrainData(std::istream& is) {
        std::vector<std::vector<float>> X;
        is.seekg(0, std::ios::end);
        std::streamsize len = is.tellg();
        is.seekg(0, std::ios::beg);

        constexpr std::size_t buf_len = 1024;

        std::array<char, buf_len> buffer = {0};
        float target = 0.0f;
        std::vector<float> row;
        std::size_t col_index = 0;
        while (is.getline(&buffer[0], buf_len)) {
            auto curr = buffer.begin();
            auto end = buffer.end();
            while (curr != end) {
                float result = 0;
                auto delim = std::find_if(curr, end, [](int ch) {
                    return ch == ',' || ch == '\n' || ch == '\0';
                });
                std::string_view sv(curr, delim);
                std::from_chars(sv.begin(), sv.end(), result);
                if (!col_index) {
                    y_.push_back(result);
                } else {
                    row.push_back(result);
                }
                curr = delim + 1;
                ++col_index;
                if (*delim == '\n') {
                    C_ = row.size();
                    X.push_back(std::move(row));
                    row = {};
                    col_index = 0;
                    target = 0.0f;
                }
            }
        }
        assert(!X.empty());
        n_ = y_.size();
        if (n_ != X.size() ||
            std::ranges::any_of(X, [&](const auto& row) { return row.size() != C_; })) {
            throw std::invalid_argument("Number of features of sample does not match");
        }
        X_ = std::move(X);
        scale();
    }

    SVMTrainData(Samples X, std::vector<float> y) : X_(std::move(X)), y_(std::move(y)), n_(y_.size()), C_(std::get<Samples>(X_)[0].size()) {
        scale();
    }

    // constructor for subsamples
    SVMTrainData(SubSamples X, std::vector<float> y, std::vector<std::pair<float, float>> feature_bounds)
    : X_(std::move(X)), y_(std::move(y)), n_(y_.size()), C_(std::get<SubSamples>(X_)[0].get().size()),
    feature_bounds_(std::move(feature_bounds)) {}

    [[nodiscard]] Samples getX() const {
        if (std::holds_alternative<Samples>(X_)) {
            return std::get<Samples>(X_);
        } else {
            throw std::runtime_error("getX() called for subsample set");
        }
    }

    [[nodiscard]] SubSamples getXRef() const {
        if (std::holds_alternative<SubSamples>(X_)) {
            return std::get<SubSamples>(X_);
        } else {
            throw std::runtime_error("getXRef() called for full sample set");
        }
    }

    [[nodiscard]] std::vector<float> getX(std::size_t row_index) const {
        if (std::holds_alternative<Samples>(X_)) {
            return std::get<Samples>(X_)[row_index];
        } else {
            return std::get<SubSamples>(X_)[row_index].get();
        }
    }

    [[nodiscard]] std::reference_wrapper<const std::vector<float>> getXRef(std::size_t row_index) const {
        if (std::holds_alternative<Samples>(X_)) {
            return std::ref(std::get<Samples>(X_)[row_index]);
        } else {
            return std::get<SubSamples>(X_)[row_index];
        }
    }

private:
    void scale() {
        if (std::holds_alternative<SubSamples>(X_)) {
            return;
        }
        feature_bounds_ = std::vector<std::pair<float, float>>(C_, {std::numeric_limits<float>::max(),
              std::numeric_limits<float>::lowest()});
        auto& X = std::get<Samples>(X_);
        for (std::size_t i = 0; i < n_; ++i) {
            for (std::size_t c = 0; c < C_; ++c) {
                feature_bounds_[c].first = std::min(feature_bounds_[c].first, X[i][c]);
                feature_bounds_[c].second = std::max(feature_bounds_[c].second, X[i][c]);
            }
        }

        for (std::size_t c = 0; c < C_; ++c) {
            auto [fmin, fmax] = feature_bounds_[c];
            for (std::size_t i = 0; i < n_; ++i) {
                X[i][c] = 1.0f - 2.0f * (fmax - X[i][c])/(fmax - fmin);
            }
        }
    }

};

} // namespace frozenca

#endif //FROZENCA_SVMDATA_H

SVMKernel.h

#ifndef FROZENCA_SVMKERNEL_H
#define FROZENCA_SVMKERNEL_H

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <functional>
#include <ranges>
#include <vector>

namespace frozenca {

static float dot(const std::vector<float>& x, const std::vector<float>& y) {
    return std::inner_product(x.begin(), x.end(), y.begin(), 0.0f);
}

class Kernel {
public:
    virtual ~Kernel() = default;
    [[nodiscard]] virtual float operator()(const std::vector<float>& x, const std::vector<float>& y) const = 0;
    virtual void adjustParams(std::size_t f) = 0;
};

class KernelLinear final : public Kernel {
public:
    [[nodiscard]] float operator()(const std::vector<float>& x, const std::vector<float>& y) const {
        return dot(x, y);
    }
    void adjustParams(std::size_t f) final {
        // do nothing
    }
};

class KernelPoly final : public Kernel {
private:
    float gamma_ = 0.0f;
    float coef0_ = 0.0f;
    std::size_t degree_ = 3;
public:
    [[nodiscard]] float operator()(const std::vector<float>& x, const std::vector<float>& y) const{
        return std::pow(gamma_ * dot(x, y) + coef0_, degree_);
    }
    void adjustParams(std::size_t f) {
        gamma_ = 1.0f / f;
    }
};

class KernelRBF final : public Kernel {
private:
    float gamma_ = 0.0f;
public:
    [[nodiscard]] float operator()(const std::vector<float>& x, const std::vector<float>& y) const final {
        std::vector<float> diff = x;
        std::ranges::transform(diff, y, diff.begin(), std::minus{});
        return std::exp(-gamma_ * dot(diff, diff));
    }
    void adjustParams(std::size_t f) final {
        gamma_ = 1.0f / f;
    }
};

class KernelSigmoid final : public Kernel {
private:
    float gamma_ = 0.0f;
    float coef0_ = 0.0f;
public:
    [[nodiscard]] float operator()(const std::vector<float>& x, const std::vector<float>& y) const final {
        return std::tanh(gamma_ * dot(x, y) + coef0_);
    }
    void adjustParams(std::size_t f) final {
        gamma_ = 1.0f / f;
    }
};

} // namespace frozenca

#endif //FROZENCA_SVMKERNEL_H

SVMModel.h

#ifndef FROZENCA_SVMMODEL_H
#define FROZENCA_SVMMODEL_H

#include "SVMData.h"
#include "SVMKernel.h"
#include <algorithm>
#include <functional>
#include <memory>
#include <random>
#include <unordered_set>
#include <unordered_map>
#include <utility>
#include <vector>

namespace frozenca {

struct SVMParams {
    float C = 1.0f;
    float nu = 0.5f;
    float p = 0.1f;
    std::vector<float> weight_scale;
    bool shrinking = true; // use the shrinking heuristics
    float eps = 1e-3;
    bool compute_probability = false;
};

class SVMModel {
protected:
    // actual model parameters
    // support vectors for each class (size : k)
    std::unordered_map<std::size_t, std::vector<std::vector<float>>> SV_;
    // coefficients for SVs in decision functions. classifier between class (i, j) => SV_coeff_[i * k + j]
    std::unordered_map<std::size_t, std::vector<float>> SV_coeff_;
    // constants in decision functions (size: k * (k - 1) / 2)
    std::vector<float> rho_;

    std::vector<std::pair<float, float>> feature_bounds_;

    std::vector<float> prob_A_;
    std::vector<float> prob_B_;

public:
    std::unique_ptr<Kernel> kernel_ = nullptr;
    friend class SVMTwoClass;

    SVMModel(std::unique_ptr<Kernel> kernel) : kernel_(std::move(kernel)) {}
    virtual ~SVMModel() = default;

    virtual void fit(const SVMTrainData& svm_data, const SVMParams& params);
    virtual std::vector<float> crossValidate(const SVMTrainData& svm_data, const SVMParams& params, std::size_t count_fold) final;

    virtual std::vector<float> scaleInputs(const std::vector<float>& sample) const final;

    [[nodiscard]] virtual float predict(const std::vector<float>& sample) const = 0;
    [[nodiscard]] virtual std::vector<float> predict(const std::vector<std::vector<float>>& samples) const final {
        std::vector<float> results;
        for (const auto& sample: samples) {
            results.push_back(predict(sample));
        }
        return results;
    }
};

void SVMModel::fit(const SVMTrainData& svm_data, const SVMParams& params) {
    feature_bounds_ = svm_data.feature_bounds_;
    kernel_->adjustParams(svm_data.C_);
}

std::vector<float> SVMModel::scaleInputs(const std::vector<float>& sample) const {
    auto sample_scaled = sample;
    for (std::size_t c = 0; c < feature_bounds_.size(); ++c) {
        auto [fmin, fmax] = feature_bounds_[c];
        sample_scaled[c] = 1.0f - 2.0f * (fmax - sample[c])/(fmax - fmin);
    }
    return sample_scaled;
}

std::vector<float> SVMModel::crossValidate(const SVMTrainData& svm_data, const SVMParams& params, std::size_t count_fold) {
    const std::size_t l = svm_data.n_;
    if (count_fold > l) {
        count_fold = l;
        std::cerr << "WARNING: #folds > #data. Will use #folds = #data (LOOCV)\n";
    }
    std::vector<std::size_t> fold_indices(l);
    for (std::size_t f = 0; f < count_fold; ++f) {
        for (std::size_t i = f * l / count_fold; i < (f + 1) * l / count_fold; ++i) {
            fold_indices[i] = f;
        }
    }
    std::mt19937 gen(std::random_device{}());
    std::ranges::shuffle(fold_indices, gen);
    std::unordered_map<std::size_t, std::unordered_set<std::size_t>> fold_sets;
    for (std::size_t i = 0; i < l; ++i) {
        fold_sets[fold_indices[i]].insert(i);
    }

    std::vector<float> target(l);

    for (std::size_t f = 0; f < count_fold; ++f) {
        // leave f-th fold out
        SubSamples X_without_f;
        std::vector<float> y_without_f;
        for (std::size_t fd = 0; fd < count_fold; ++fd) {
            if (fd == f) {
                continue;
            }
            auto& fold_fd = fold_sets[fd];
            for (auto fd_sample : fold_fd) {
                X_without_f.push_back(svm_data.getXRef(fd_sample));
                y_without_f.push_back(svm_data.y_[fd_sample]);
            }
        }

        SVMTrainData CV_data_ffold (X_without_f, y_without_f, svm_data.feature_bounds_);
        fit(CV_data_ffold, params);
        auto& fold_f = fold_sets[f];
        for (auto f_sample : fold_f) {
            target[f_sample] = predict(svm_data.getX(f_sample));
        }
    }
    return target;
}

} // namespace frozenca

#endif //FROZENCA_SVMMODEL_H

SVMSolver.h

#ifndef FROZENCA_SVMSOLVER_H
#define FROZENCA_SVMSOLVER_H

#include "../../Matrix/Matrix.h"
#include <cstddef>
#include <tuple>
#include <unordered_set>
#include <utility>
#include <vector>

namespace frozenca {

static constexpr float TAU = 1e-12;

struct Solution {
    std::vector<float> alpha;
    float obj = 0.0f;
    float rho = 0.0f;
    float r = 0.0f; // for Solver_NU
    std::pair<float, float> upper_bound = {0.0f, 0.0f};
};

class Solver {
public:
    enum class Alpha {
        Lower,
        Upper,
        Free
    };

    Solver(std::size_t l, Mat<float>& Q, std::vector<float> QD, std::vector<float> p, std::vector<char> y,
           std::vector<float> alpha, std::pair<float, float> C, float eps, bool shrinking);

    virtual ~Solver() = default;

protected:
    const std::size_t l_;
    Mat<float>& Q_;
    std::vector<float> QD_;
    std::vector<float> p_;
    std::vector<char> y_;
    std::vector<float> alpha_;
    std::pair<float, float> C_;
    float eps_;
    std::vector<Alpha> alpha_status_;
    std::vector<float> G_; // gradient of objective function
    std::vector<float> G_bar_; // gradient, if we treat free variable as 0
    bool shrinking_;
    bool unshrink_ = false;
    std::unordered_set<std::size_t> active_set_;

public:
    [[nodiscard]] float getC(std::size_t i) const {
        return (y_[i] > 0) ? C_.first : C_.second;
    }

    [[nodiscard]] std::vector<float> getQ(std::size_t i) const {
        auto Q_i = Q_.col(i);
        std::vector<float> Qi;
        for (auto q : Q_i) {
            Qi.push_back(q);
        }
        return Qi;
    }

    void updateAlphaStatus(std::size_t i) {
        if (alpha_[i] >= getC(i)) {
            alpha_status_[i] = Alpha::Upper;
        } else if (alpha_[i] <= 0) {
            alpha_status_[i] = Alpha::Lower;
        } else {
            alpha_status_[i] = Alpha::Free;
        }
    }

public:
    virtual Solution Solve() final;
    virtual void updateAlphaGradientValues(std::size_t i, std::size_t j) final;
    virtual void reconstructGradient() final;

protected:
    virtual std::tuple<bool, std::size_t, std::size_t> selectWorkingSet() const = 0;
    virtual void calculateRho(Solution& sol) const = 0;
    virtual void doShrinking() = 0;

};

Solver::Solver(std::size_t l, Mat<float>& Q, std::vector<float> QD,
               std::vector<float> p, std::vector<char> y,
               std::vector<float> alpha, std::pair<float, float> C, float eps, bool shrinking) :
        l_{l}, Q_(Q), QD_(std::move(QD)), p_(std::move(p)), y_(std::move(y)),
        alpha_(std::move(alpha)), C_(std::move(C)), eps_{eps}, shrinking_{shrinking} {
    alpha_status_.resize(l_);
    for (std::size_t i = 0; i < l_; ++i) {
        active_set_.insert(i);
        updateAlphaStatus(i);
    }
    G_ = p_;
    G_bar_.resize(l_);
    for (std::size_t i = 0; i < l_; ++i) {
        if (alpha_status_[i] != Alpha::Lower) {
            auto Q_i = getQ(i);
            for (std::size_t j = 0; j < l_; ++j) {
                G_[j] += alpha_[i] * Q_i[j];
            }
            if (alpha_status_[i] == Alpha::Upper) {
                for (std::size_t j = 0; j < l_; ++j) {
                    G_bar_[j] += getC(i) * Q_i[j];
                }
            }
        }
    }
}

// Solver ordinary (non NU)

class SolverOrdinary final : public Solver {
public:
    SolverOrdinary(std::size_t l, Mat<float>& Q, std::vector<float> QD, std::vector<float> p, std::vector<char> y,
                   std::vector<float> alpha, std::pair<float, float> C, float eps, bool shrinking)
            : Solver(l, Q, QD, p, y, alpha, C, eps, shrinking) {}

    std::tuple<bool, std::size_t, std::size_t> selectWorkingSet() const;
    void calculateRho(Solution& sol) const;
    void doShrinking() final;
    bool beShrunk(std::size_t i, float g_max1, float g_max2) const;
};

Solution Solver::Solve() {
    std::size_t iter = 0;
    std::size_t max_iter = std::max(1'000'000lu, l_ > std::numeric_limits<std::size_t>::max() / 100 ?
                                                 std::numeric_limits<std::size_t>::max() : 100 * l_);
    std::size_t counter = std::min(l_, 1'000lu) + 1;

    while (iter < max_iter) {
        if (!--counter) {
            counter = std::min(l_, 1'000lu);
            if (shrinking_) {
                doShrinking();
            }
        }
        bool already_opt = false;
        std::size_t i = 0, j = 0;
        std::tie(already_opt, i, j) = selectWorkingSet();
        if (already_opt) {
            // reconstruct the whole gradient
            reconstructGradient();
            std::tie(already_opt, i, j) = selectWorkingSet();
            if (already_opt) {
                break;
            } else {
                counter = 1; // do shrinking in next iteration
            }
        }
        ++iter;

        updateAlphaGradientValues(i, j);
    }

    if (iter >= max_iter) {
        if (active_set_.size() < l_) {
            reconstructGradient();
        }
        std::cerr << "WARNING: reaching max number of iterations\n";
    }

    Solution s;
    calculateRho(s);
    for (std::size_t i = 0; i < l_; ++i) {
        s.obj += alpha_[i] * (G_[i] + p_[i]);
    }
    s.obj /= 2.0f;
    s.alpha = alpha_;
    s.upper_bound = C_;

    std::cout << "Optimization finished, #iter = " << iter << '\n';
    return s;
}

void Solver::updateAlphaGradientValues(std::size_t i, std::size_t j) {
    assert(active_set_.contains(i) && active_set_.contains(j));
    // update alpha
    auto Q_i = getQ(i);
    auto Q_j = getQ(j);
    auto C_i = getC(i);
    auto C_j = getC(j);
    auto old_alpha_i = alpha_[i];
    auto old_alpha_j = alpha_[j];

    if (y_[i] != y_[j]) {
        float quad_coeff = QD_[i] + QD_[j] + 2.0f * Q_i[j];
        if (quad_coeff <= 0) {
            quad_coeff = TAU;
        }
        float delta = (-G_[i] - G_[j]) / quad_coeff;
        float diff = alpha_[i] - alpha_[j];
        alpha_[i] += delta;
        alpha_[j] += delta;
        if (diff > 0) {
            if (alpha_[j] > 0) {
                alpha_[j] = 0;
                alpha_[i] = diff;
            }
        } else {
            if (alpha_[i] < 0) {
                alpha_[i] = 0;
                alpha_[j] = -diff;
            }
        }
        if (diff > C_i - C_j)  {
            if (alpha_[i] > C_i) {
                alpha_[i] = C_i;
                alpha_[j] = C_i - diff;
            }
        } else {
            if (alpha_[j] > C_j) {
                alpha_[j] = C_j;
                alpha_[i] = C_j + diff;
            }
        }
    } else {
        float quad_coeff = QD_[i] + QD_[j] - 2.0f * Q_i[j];
        if (quad_coeff <= 0) {
            quad_coeff = TAU;
        }
        float delta = (G_[i] - G_[j]) / quad_coeff;
        float sum = alpha_[i] + alpha_[j];
        alpha_[i] -= delta;
        alpha_[j] += delta;
        if (sum > C_i) {
            if (alpha_[i] > C_i) {
                alpha_[i] = C_i;
                alpha_[j] = sum - C_i;
            }
        } else {
            if (alpha_[j] < 0) {
                alpha_[j] = 0;
                alpha_[i] = sum;
            }
        }
        if (sum > C_j)  {
            if (alpha_[i] > C_j) {
                alpha_[i] = C_j;
                alpha_[j] = sum - C_j;
            }
        } else {
            if (alpha_[i] < 0) {
                alpha_[i] = 0;
                alpha_[j] = sum;
            }
        }
    }

    // update gradient
    float delta_alpha_i = alpha_[i] - old_alpha_i;
    float delta_alpha_j = alpha_[j] - old_alpha_j;

    for (auto k : active_set_) {
        G_[k] += Q_i[k] * delta_alpha_i + Q_j[k] * delta_alpha_j;
    }

    // update alpha status and G_bar
    bool ui = alpha_status_[i] == Alpha::Upper;
    bool uj = alpha_status_[j] == Alpha::Upper;
    updateAlphaStatus(i);
    updateAlphaStatus(j);
    if (ui != (alpha_status_[i] == Alpha::Upper)) {
        Q_i = getQ(i);
        if (ui) {
            for (std::size_t k = 0; k < l_; ++k) {
                G_bar_[k] -= C_i * Q_i[k];
            }
        } else {
            for (std::size_t k = 0; k < l_; ++k) {
                G_bar_[k] += C_i * Q_i[k];
            }
        }
    }
    if (uj != (alpha_status_[j] == Alpha::Upper)) {
        Q_j = getQ(j);
        if (uj) {
            for (std::size_t k = 0; k < l_; ++k) {
                G_bar_[k] -= C_j * Q_j[k];
            }
        } else {
            for (std::size_t k = 0; k < l_; ++k) {
                G_bar_[k] += C_j * Q_j[k];
            }
        }
    }
}

void Solver::reconstructGradient() {
    if (active_set_.size() == l_) {
        return;
    }
    std::unordered_set<std::size_t> inactive_set;
    for (std::size_t j = 0; j < l_; ++j) {
        if (!active_set_.contains(j)) {
            inactive_set.insert(j);
        }
    }

    std::size_t free_count = 0;

    for (std::size_t j = 0; j < l_; ++j) {
        if (active_set_.contains(j)) {
            if (alpha_status_[j] == Alpha::Free) {
                ++free_count;
            }
        } else {
            G_[j] = G_bar_[j] + p_[j];
        }
    }

    if (2 * free_count < active_set_.size()) {
        std::cerr << "WARNING: deactivating shrinking may be faster\n";
    }

    if (free_count * l_ > 2 * active_set_.size() * (l_ - active_set_.size())) {
        for (auto i : inactive_set) {
            auto Q_i = getQ(i);
            for (auto j : active_set_) {
                if (alpha_status_[j] == Alpha::Free) {
                    G_[i] += alpha_[j] * Q_i[j];
                }
            }
        }
    } else {
        for (auto i : active_set_) {
            if (alpha_status_[i] == Alpha::Free) {
                auto Q_i = getQ(i);
                for (auto j : inactive_set) {
                    G_[j] += alpha_[i] * Q_i[j];
                }
            }
        }
    }
    for (std::size_t i = 0; i < l_; ++i) {
        active_set_.insert(i);
    }
}

// [0] : return whether if already optimal
// [1], [2] : return i, j such that
// i : maximizes -y_i * grad(f)_i in I_high(alpha)
// j : minimizes the decreases of obj value
// (if quadratic coefficient <= 0, replace it with tau)
// -y_j * grad(f)_j < -y_i * grad(f)_i, j in I_low(alpha)
std::tuple<bool, std::size_t, std::size_t> SolverOrdinary::selectWorkingSet() const {
    float g_max1 = std::numeric_limits<float>::lowest();
    std::size_t g_max_idx = -1;
    for (auto i : active_set_) {
        if (y_[i] == +1) {
            if (alpha_status_[i] != Alpha::Upper) {
                if (g_max1 <= -G_[i]) {
                    g_max1 = -G_[i];
                    g_max_idx = i;
                }
            }
        } else {
            if (alpha_status_[i] != Alpha::Lower) {
                if (g_max1 <= G_[i]) {
                    g_max1 = G_[i];
                    g_max_idx = i;
                }
            }
        }
    }

    std::size_t i = g_max_idx;
    std::vector<float> Qi;
    if (i != -1) {
        Qi = getQ(i);
    }
    std::size_t g_min_idx = -1;
    float g_max2 = std::numeric_limits<float>::lowest();
    float obj_diff_min = std::numeric_limits<float>::max();
    for (auto j : active_set_) {
        if (y_[j] == +1) {
            if (alpha_status_[j] != Alpha::Lower) {
                float grad_diff = g_max1 + G_[j];
                if (G_[j] >= g_max2) {
                    g_max2 = G_[j];
                }
                if (grad_diff > 0) {
                    float obj_diff = 0.0f;
                    float quad_coeff = QD_[i] + QD_[j] - 2.0f * y_[i] * Qi[j];
                    if (quad_coeff > 0) {
                        obj_diff = -std::pow(grad_diff, 2.0f) / quad_coeff;
                    } else {
                        obj_diff = -std::pow(grad_diff, 2.0f) / TAU;
                    }
                    if (obj_diff <= obj_diff_min) {
                        g_min_idx = j;
                        obj_diff_min = obj_diff;
                    }
                }
            }
        } else {
            if (alpha_status_[j] != Alpha::Upper) {
                float grad_diff = g_max1 - G_[j];
                if (-G_[j] >= g_max2) {
                    g_max2 = -G_[j];
                }
                if (grad_diff > 0) {
                    float obj_diff = 0.0f;
                    float quad_coeff = QD_[i] + QD_[j] + 2.0f * y_[i] * Qi[j];
                    if (quad_coeff > 0) {
                        obj_diff = -std::pow(grad_diff, 2.0f) / quad_coeff;
                    } else {
                        obj_diff = -std::pow(grad_diff, 2.0f) / TAU;
                    }
                    if (obj_diff <= obj_diff_min) {
                        g_min_idx = j;
                        obj_diff_min = obj_diff;
                    }
                }
            }
        }
    }
    if (g_max1 + g_max2 < eps_ || g_min_idx == -1) {
        return {true, -1, -1};
    }
    return {false, g_max_idx, g_min_idx};
}

void SolverOrdinary::calculateRho(Solution& sol) const {
    float res = 0.0f;
    std::size_t count_free = 0;
    float ub = std::numeric_limits<float>::max();
    float lb = std::numeric_limits<float>::lowest();
    float sum_free = 0.0f;
    for (auto i : active_set_) {
        auto yG = y_[i] * G_[i];
        if (alpha_status_[i] == Alpha::Upper) {
            if (y_[i] == -1) {
                ub = std::min(ub, yG);
            } else {
                lb = std::max(lb, yG);
            }
        } else if (alpha_status_[i] == Alpha::Lower) {
            if (y_[i] == +1) {
                ub = std::min(ub, yG);
            } else {
                lb = std::max(lb, yG);
            }
        } else {
            ++count_free;
            sum_free += yG;
        }
    }
    if (count_free) {
        res = sum_free / count_free;
    } else {
        res = (ub + lb) / 2.0f;
    }
    sol.rho = res;
}

void SolverOrdinary::doShrinking() {
    float g_max1 = std::numeric_limits<float>::lowest();
    float g_max2 = std::numeric_limits<float>::lowest();

    for (auto i : active_set_) {
        if (y_[i] == +1) {
            if (alpha_status_[i] != Alpha::Upper) {
                g_max1 = std::max(g_max1, -G_[i]);
            }
            if (alpha_status_[i] != Alpha::Lower) {
                g_max2 = std::max(g_max2, G_[i]);
            }
        } else {
            if (alpha_status_[i] != Alpha::Upper) {
                g_max2 = std::max(g_max2, -G_[i]);
            }
            if (alpha_status_[i] != Alpha::Lower) {
                g_max1 = std::max(g_max1, G_[i]);
            }
        }
    }
    if (!unshrink_ && g_max1 + g_max2 <= eps_ * 10) {
        unshrink_ = true;
        reconstructGradient();
    }

    std::unordered_set<std::size_t> to_shrunk;
    for (auto i : active_set_) {
        if (beShrunk(i, g_max1, g_max2)) {
            to_shrunk.insert(i);
        }
    }
    for (auto i : to_shrunk) {
        active_set_.erase(i);
    }
}

bool SolverOrdinary::beShrunk(std::size_t i, float g_max1, float g_max2) const {
    assert(i < l_);
    if (alpha_status_[i] == Alpha::Upper) {
        if (y_[i] == +1) {
            return (-G_[i] > g_max1);
        } else {
            return (-G_[i] > g_max2);
        }
    } else if (alpha_status_[i] == Alpha::Lower) {
        if (y_[i] == +1) {
            return (G_[i] > g_max2);
        } else {
            return (G_[i] > g_max1);
        }
    }
    return false;
}


// Solver NU

class SolverNU final : public Solver {
public:
    SolverNU(std::size_t l, Mat<float>& Q, std::vector<float> QD, std::vector<float> p, std::vector<char> y,
             std::vector<float> alpha, std::pair<float, float> C, float eps, bool shrinking)
            : Solver(l, Q, std::move(QD), std::move(p), std::move(y), std::move(alpha), C, eps, shrinking) {}

    std::tuple<bool, std::size_t, std::size_t> selectWorkingSet() const;
    void calculateRho(Solution& sol) const;
    void doShrinking();
    bool beShrunk(std::size_t i, float g_max1, float g_max2, float g_max3, float g_max4) const;
};

// [0] : return whether if already optimal
// [1], [2] : return i, j such that
// i : maximizes -y_i * grad(f)_i in I_high(alpha)
// j : minimizes the decreases of obj value
// (if quadratic coefficient <= 0, replace it with tau)
// -y_j * grad(f)_j < -y_i * grad(f)_i, j in I_low(alpha)
std::tuple<bool, std::size_t, std::size_t> SolverNU::selectWorkingSet() const {
    float g_maxp1 = std::numeric_limits<float>::lowest();
    float g_maxn1 = std::numeric_limits<float>::lowest();
    std::size_t g_maxp_idx = -1;
    std::size_t g_maxn_idx = -1;
    for (auto i : active_set_) {
        if (y_[i] == +1) {
            if (alpha_status_[i] != Alpha::Upper) {
                if (g_maxp1 <= -G_[i]) {
                    g_maxp1 = -G_[i];
                    g_maxp_idx = i;
                }
            }
        } else {
            if (alpha_status_[i] != Alpha::Lower) {
                if (g_maxn1 <= G_[i]) {
                    g_maxn1 = G_[i];
                    g_maxn_idx = i;
                }
            }
        }
    }

    std::size_t ip = g_maxp_idx;
    std::vector<float> Qip;
    if (ip != -1) {
        Qip = getQ(ip);
    }
    std::size_t in = g_maxn_idx;
    std::vector<float> Qin;
    if (in != -1) {
        Qin = getQ(in);
    }
    std::size_t g_min_idx = -1;
    float g_maxp2 = std::numeric_limits<float>::lowest();
    float g_maxn2 = std::numeric_limits<float>::lowest();
    float obj_diff_min = std::numeric_limits<float>::max();
    for (auto j : active_set_) {
        if (y_[j] == +1) {
            if (alpha_status_[j] != Alpha::Lower) {
                float grad_diff = g_maxp1 + G_[j];
                if (G_[j] >= g_maxp2) {
                    g_maxp2 = G_[j];
                }
                if (grad_diff > 0) {
                    float obj_diff = 0.0f;
                    float quad_coeff = QD_[ip] + QD_[j] - 2.0f * Qip[j];
                    if (quad_coeff > 0) {
                        obj_diff = -std::pow(grad_diff, 2.0f) / quad_coeff;
                    } else {
                        obj_diff = -std::pow(grad_diff, 2.0f) / TAU;
                    }
                    if (obj_diff <= obj_diff_min) {
                        g_min_idx = j;
                        obj_diff_min = obj_diff;
                    }
                }
            }
        } else {
            if (alpha_status_[j] != Alpha::Upper) {
                float grad_diff = g_maxn1 - G_[j];
                if (-G_[j] >= g_maxn2) {
                    g_maxn2 = -G_[j];
                }
                if (grad_diff > 0) {
                    float obj_diff = 0.0f;
                    float quad_coeff = QD_[in] + QD_[j] - 2.0f * Qin[j];
                    if (quad_coeff > 0) {
                        obj_diff = -std::pow(grad_diff, 2.0f) / quad_coeff;
                    } else {
                        obj_diff = -std::pow(grad_diff, 2.0f) / TAU;
                    }
                    if (obj_diff <= obj_diff_min) {
                        g_min_idx = j;
                        obj_diff_min = obj_diff;
                    }
                }
            }
        }
    }
    if (std::max(g_maxp1 + g_maxp2, g_maxn1 + g_maxn2) < eps_ || g_min_idx == -1) {
        return {true, -1, -1};
    }
    return {false, (y_[g_min_idx] == +1) ? g_maxp_idx : g_maxn_idx, g_min_idx};
}

void SolverNU::calculateRho(Solution& sol) const {
    float res = 0.0f;
    std::size_t count_free1 = 0;
    std::size_t count_free2 = 0;
    float ub1 = std::numeric_limits<float>::max();
    float ub2 = std::numeric_limits<float>::max();
    float lb1 = std::numeric_limits<float>::lowest();
    float lb2 = std::numeric_limits<float>::lowest();
    float sum_free1 = 0.0f;
    float sum_free2 = 0.0f;
    for (auto i : active_set_) {
        if (y_[i] == +1) {
            if (alpha_status_[i] == Alpha::Upper) {
                lb1 = std::max(lb1, G_[i]);
            } else if (alpha_status_[i] == Alpha::Lower) {
                ub1 = std::min(ub1, G_[i]);
            } else {
                ++count_free1;
                sum_free1 += G_[i];
            }
        } else {
            if (alpha_status_[i] == Alpha::Upper) {
                lb2 = std::max(lb2, G_[i]);
            } else if (alpha_status_[i] == Alpha::Lower) {
                ub2 = std::min(ub2, G_[i]);
            } else {
                ++count_free2;
                sum_free2 += G_[i];
            }
        }
    }
    float r1 = 0.0f;
    float r2 = 0.0f;
    if (count_free1) {
        r1 = sum_free1 / count_free1;
    } else {
        r1 = (ub1 + lb1) / 2.0f;
    }
    if (count_free2) {
        r2 = sum_free2 / count_free2;
    } else {
        r2 = (ub2 + lb2) / 2.0f;
    }
    sol.r = (r1 + r2) / 2.0f;
    sol.rho = (r1 - r2) / 2.0f;
}

void SolverNU::doShrinking() {
    float g_max1 = std::numeric_limits<float>::lowest();
    float g_max2 = std::numeric_limits<float>::lowest();
    float g_max3 = std::numeric_limits<float>::lowest();
    float g_max4 = std::numeric_limits<float>::lowest();

    for (auto i : active_set_) {
        if (alpha_status_[i] != Alpha::Upper) {
            if (y_[i] == +1) {
                g_max1 = std::max(g_max1, -G_[i]);
            } else {
                g_max4 = std::max(g_max4, -G_[i]);
            }
        }
        if (alpha_status_[i] != Alpha::Lower) {
            if (y_[i] == +1) {
                g_max2 = std::max(g_max2, G_[i]);
            } else {
                g_max3 = std::max(g_max3, G_[i]);
            }
        }
    }
    if (!unshrink_ && std::max(g_max1 + g_max2, g_max3 + g_max4) <= eps_ * 10) {
        unshrink_ = true;
        reconstructGradient();
    }

    std::unordered_set<std::size_t> to_shrunk;
    for (auto i : active_set_) {
        if (beShrunk(i, g_max1, g_max2, g_max3, g_max4)) {
            to_shrunk.insert(i);
        }
    }
    for (auto i : to_shrunk) {
        active_set_.erase(i);
    }
}

bool SolverNU::beShrunk(std::size_t i, float g_max1, float g_max2,
                        float g_max3, float g_max4) const {
    assert(i < l_);
    if (alpha_status_[i] == Alpha::Upper) {
        if (y_[i] == +1) {
            return (-G_[i] > g_max1);
        } else {
            return (-G_[i] > g_max4);
        }
    } else if (alpha_status_[i] == Alpha::Lower) {
        if (y_[i] == +1) {
            return (G_[i] > g_max2);
        } else {
            return (G_[i] > g_max3);
        }
    }
    return false;
}

} // namespace frozenca

#endif //FROZENCA_SVMSOLVER_H

SVMClassification.h

#ifndef FROZENCA_SVMCLASSIFICATION_H
#define FROZENCA_SVMCLASSIFICATION_H

#include "SVMModel.h"
#include "SVMSolver.h"
#include <algorithm>
#include <functional>
#include <ranges>
#include <unordered_set>
#include <unordered_map>

namespace frozenca {

class SVMClassification : public SVMModel {
private:
    std::unordered_map<int, std::unordered_set<std::size_t>> groups_;
    std::vector<char> nonzero_;
    std::vector<int> data_label_;
    std::vector<int> labels_;
public:
    SVMClassification(std::unique_ptr<Kernel> kernel) : SVMModel(std::move(std::move(kernel))) {}
    virtual ~SVMClassification() = default;
    virtual void fit(const SVMTrainData& svm_data, const SVMParams& params) final;
    virtual void constructGroups(const SVMTrainData& svm_data) final;
    virtual void constructModel(const SVMTrainData& svm_data,
                        const std::vector<std::pair<std::vector<float>, float>>& decision_functions) final;
    virtual std::pair<Mat<float>, std::vector<float>> computeQs(std::size_t l, const SubSamples& X, const std::vector<char>& y) const final;
    [[nodiscard]] virtual float predict(const std::vector<float>& sample) const final;
protected:
    virtual std::pair<std::vector<float>, float> fitOne(const SVMTrainData& svm_sub_data, const SVMParams& params,
                                                        float Cp, float Cn) = 0;
};

void SVMClassification::fit(const SVMTrainData& svm_data, const SVMParams& params) {
    SVMModel::fit(svm_data, params);
    constructGroups(svm_data);
    const std::size_t l = svm_data.n_;
    const std::size_t k = groups_.size();

    // train k * (k - 1) / 2 binary classifiers
    nonzero_.clear();
    nonzero_.resize(l);

    std::vector<SubSamples> class_X (k);

    for (std::size_t i = 0; i < k; ++i) {
        auto label_i = labels_[i];
        auto& group_i = groups_[label_i];
        for (auto index_i : group_i) {
            class_X[i].push_back(svm_data.getXRef(index_i));
        }
    }

    std::vector<std::pair<std::vector<float>, float>> decision_functions;
    for (std::size_t i = 0; i < k; ++i) {
        auto& group_i = groups_[labels_[i]];
        for (std::size_t j = i + 1; j < k; ++j) {
            auto& group_j = groups_[labels_[j]];
            SubSamples sample_ij;
            std::ranges::copy(class_X[i], std::back_inserter(sample_ij));
            std::ranges::copy(class_X[j], std::back_inserter(sample_ij));
            std::vector<float> label_ij;
            for (std::size_t ci = 0; ci < class_X[i].size(); ++ci) {
                label_ij.push_back(+1);
            }
            for (std::size_t cj = 0; cj < class_X[j].size(); ++cj) {
                label_ij.push_back(-1);
            }
            SVMTrainData sub_data(sample_ij, label_ij, svm_data.feature_bounds_);
            auto wi = params.weight_scale.empty() ? 1.0f : params.weight_scale[i];
            auto wj = params.weight_scale.empty() ? 1.0f : params.weight_scale[j];
            auto f = fitOne(sub_data, params, params.C * wi, params.C * wj);
            std::size_t alpha_index = 0;
            for (auto index_i : group_i) {
                if (!nonzero_[index_i] && std::fabs(f.first[alpha_index]) > 0) {
                    nonzero_[index_i] = true;
                }
                alpha_index++;
            }
            for (auto index_j : group_j) {
                if (!nonzero_[index_j] && std::fabs(f.first[alpha_index]) > 0) {
                    nonzero_[index_j] = true;
                }
                alpha_index++;
            }
            decision_functions.push_back(std::move(f));
        }
    }
    constructModel(svm_data, decision_functions);
}

void SVMClassification::constructGroups(const SVMTrainData& svm_data) {
    if (!groups_.empty()) {
        return;
    }
    const std::size_t l = svm_data.n_;
    data_label_.resize(l);
    for (std::size_t i = 0; i < l; ++i) {
        int this_label = static_cast<int>(svm_data.y_[i]);
        groups_[this_label].insert(i);
        data_label_[i] = this_label;
    }
    for (const auto& [label, _] : groups_) {
        labels_.push_back(label);
    }
    if (groups_.size() == 1) {
        std::cerr << "WARNING: training data in only one class\n";
    }
}

void SVMClassification::constructModel(const SVMTrainData& svm_data,
                    const std::vector<std::pair<std::vector<float>, float>>& decision_functions) {
    const std::size_t k = groups_.size();

    SV_.clear();
    SV_coeff_.clear();
    rho_.clear();

    for (const auto& [alpha, rho] : decision_functions) {
        rho_.push_back(rho);
    }
    for (std::size_t i = 0; i < k; ++i) {
        for (auto idx : groups_[labels_[i]]) {
            if (nonzero_[idx]) {
                SV_[i].push_back(svm_data.getX(idx));
            }
        }
    }
    std::size_t decision_function_index = 0;
    for (std::size_t i = 0; i < k; ++i) {
        auto& group_i = groups_[labels_[i]];
        for (std::size_t j = i + 1; j < k; ++j) {
            auto& group_j = groups_[labels_[j]];
            auto& alpha = decision_functions[decision_function_index].first;
            std::size_t alpha_index = 0;
            for (auto index_i : group_i) {
                if (nonzero_[index_i]) {
                    SV_coeff_[i * k + j].push_back(alpha[alpha_index]);
                }
                ++alpha_index;
            }
            for (auto index_j : group_j) {
                if (nonzero_[index_j]) {
                    SV_coeff_[i * k + j].push_back(alpha[alpha_index]);
                }
                ++alpha_index;
            }
            ++decision_function_index;
        }
    }
}

std::pair<Mat<float>, std::vector<float>> SVMClassification::computeQs(std::size_t l,
                                                                       const SubSamples& X,
                                                                       const std::vector<char>& y) const {
    Mat<float> Q (l, l);
    std::vector<float> QD (l);
    for (std::size_t i = 0; i < l; ++i) {
        for (std::size_t j = 0; j < l; ++j) {
            Q[{i, j}] = y[i] * y[j] * (*kernel_)(X[i], X[j]);
        }
        QD[i] = Q[{i, i}];
    }
    return {Q, QD};
}

float SVMClassification::predict(const std::vector<float>& sample) const {
    auto sample_scaled = scaleInputs(sample);
    const std::size_t k = groups_.size();
    std::vector<std::size_t> vote(k);

    std::unordered_map<std::size_t, std::vector<float>> k_values;
    for (std::size_t i = 0; i < k; ++i) {
        std::size_t SV_index = 0;
        for (auto idx : groups_.at(labels_[i])) {
            if (nonzero_[idx]) {
                k_values[i].push_back((*kernel_)(sample_scaled, SV_.at(i)[SV_index++]));
            }
        }
    }
    std::size_t rho_index = 0;
    for (std::size_t i = 0; i < k; ++i) {
        auto& group_i = groups_.at(labels_[i]);
        for (std::size_t j = i + 1; j < k; ++j) {
            float sum = 0.0f;
            auto& group_j = groups_.at(labels_[j]);
            auto& curr_coeffs = SV_coeff_.at(i * k + j);
            std::size_t idx_i = 0;
            for (auto idx : group_i) {
                if (nonzero_[idx]) {
                    sum += curr_coeffs[idx_i] * k_values[i][idx_i];
                    ++idx_i;
                }
            }
            std::size_t idx_j = 0;
            for (auto idx : group_j) {
                if (nonzero_[idx]) {
                    sum += curr_coeffs[idx_i + idx_j] * k_values[j][idx_j];
                    ++idx_j;
                }
            }
            sum -= rho_[rho_index++];
            if (sum > 0.0f) {
                ++vote[i];
            } else {
                ++vote[j];
            }
        }
    }
    auto vote_max_index = std::distance(vote.begin(), std::ranges::max_element(vote));
    return static_cast<float>(labels_[vote_max_index]);
}

class SVMCSVC final : public SVMClassification {
public:
    SVMCSVC(std::unique_ptr<Kernel> kernel) : SVMClassification(std::move(kernel)) {}

    std::pair<std::vector<float>, float> fitOne(const SVMTrainData& svm_data, const SVMParams& params,
                                                float Cp, float Cn);
};

std::pair<std::vector<float>, float> SVMCSVC::fitOne(const SVMTrainData& svm_data,
                                                     const SVMParams& params,
                                                     float Cp, float Cn) {
    const std::size_t l = svm_data.n_;
    std::vector<float> minus_ones (l, -1.0f);
    std::vector<char> y (l);
    std::vector<float> alpha (l);
    for (std::size_t i = 0; i < l; ++i) {
        if (svm_data.y_[i] > 0) {
            y[i] = +1;
        } else {
            y[i] = -1;
        }
    }

    auto [Q, QD] = computeQs(l, svm_data.getXRef(), y);
    SolverOrdinary s(l, Q, QD, minus_ones, y, alpha, {Cp, Cn}, params.eps, params.shrinking);
    auto solution = s.Solve();

    if (Cp == Cn) {
        float sum_alpha = std::accumulate(solution.alpha.begin(), solution.alpha.end(), 0.0f);
        std::cout << "nu = " << sum_alpha / (Cp * l) << '\n';
    }
    std::ranges::transform(solution.alpha, y, solution.alpha.begin(), std::multiplies{});
    return {solution.alpha, solution.rho};
}

class SVMNuSVC final : public SVMClassification {
public:
    SVMNuSVC(std::unique_ptr<Kernel> kernel, const SVMParams& params) : SVMClassification(std::move(kernel)) {}
    std::pair<std::vector<float>, float> fitOne(const SVMTrainData& svm_data, const SVMParams& params,
                                                float Cp, float Cn);
};


std::pair<std::vector<float>, float> SVMNuSVC::fitOne(const SVMTrainData& svm_data, const SVMParams& params,
                                                     float Cp, float Cn) {
    const std::size_t l = svm_data.n_;

    std::vector<float> alpha (l);
    std::vector<char> y (l);
    auto sum_pos = params.nu * l / 2.0f;
    auto sum_neg = params.nu * l / 2.0f;
    for (std::size_t i = 0; i < l; ++i) {
        if (svm_data.y_[i] > 0) {
            y[i] = +1;
            alpha[i] = std::min(1.0f, sum_pos);
            sum_pos -= alpha[i];
        } else {
            y[i] = -1;
            alpha[i] = std::min(1.0f, sum_neg);
            sum_neg -= alpha[i];
        }
    }

    std::vector<float> zeros(l);

    auto [Q, QD] = computeQs(l, svm_data.getXRef(), y);
    SolverOrdinary s(l, Q, QD, zeros, y, alpha, {1.0f, 1.0f}, params.eps, params.shrinking);
    auto solution = s.Solve();

    auto r = solution.r;
    std::cout << "C = " << 1.0f / r << '\n';

    std::ranges::transform(solution.alpha, y, solution.alpha.begin(), [&r](auto a, auto y) {
        return a * y / r;
    });
    solution.rho /= r;
    solution.obj /= std::pow(r, 2.0f);
    solution.upper_bound = {1.0f / r, 1.0f / r};
    return {solution.alpha, solution.rho};
}

} // namespace frozenca

Test code: Generated male/female height, weights and trained. Predicted newly generated male/female heights and weights, and it predicts correctly.

#include "Matrix/Matrix.h"
#include "ML/SVM/SupportVectorMachine.h"
#include <iostream>
#include <memory>
#include <random>

namespace fc = frozenca;

int main() {
    std::unique_ptr<fc::SVMModel> sm = std::make_unique<fc::SVMCSVC>(std::make_unique<fc::KernelRBF>());

    std::mt19937 gen(std::random_device{}());

    constexpr std::size_t num_samples = 200;

    std::normal_distribution<float> male_height (177.0f, 8.0f);
    std::normal_distribution<float> male_weight (70.0f, 6.0f);

    std::normal_distribution<float> female_height (162.0f, 7.0f);
    std::normal_distribution<float> female_weight (53.0f, 5.5f);

    std::vector<std::vector<float>> X;
    std::vector<float> y;

    for (std::size_t i = 0; i < num_samples; ++i) {
        auto height = male_height(gen);
        auto weight = male_weight(gen);
        X.push_back(std::vector<float>{height, weight});
        y.push_back(+1);
    }
    for (std::size_t i = 0; i < num_samples; ++i) {
        auto height = female_height(gen);
        auto weight = female_weight(gen);
        X.push_back(std::vector<float>{height, weight});
        y.push_back(-1);
    }
    fc::SVMTrainData data(X, y);
    fc::SVMParams params;
    params.shrinking = false;
    sm->fit(data, params);

    std::vector<std::vector<float>> X1;

    for (std::size_t i = 0; i < 10; ++i) {
        auto height = male_height(gen);
        auto weight = male_weight(gen);
        X1.push_back({height, weight});
    }

    auto y1 = sm->predict(X1);

    for (auto y_pred : y1) {
        std::cout << y_pred << ' ';
    }
    std::cout << '\n';

    std::vector<std::vector<float>> X2;

    for (std::size_t i = 0; i < 10; ++i) {
        auto height = female_height(gen);
        auto weight = female_weight(gen);
        X2.push_back({height, weight});
    }

    auto y2 = sm->predict(X2);

    for (auto y_pred : y2) {
        std::cout << y_pred << ' ';
    }
    std::cout << '\n';

}

Test Output

Optimization finished, #iter = 39
nu = 0.18
1 1 1 1 1 1 1 1 1 1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1

Pretty good!

Feel free to comment anything!

\$\endgroup\$
1
  • \$\begingroup\$ I would suggest the you use templates so that the same code can be used for float and double. \$\endgroup\$
    – jdt
    Nov 7 at 12:28
2
\$\begingroup\$

Split SVMTrainData instead of using a variant

I think you can avoid the variant by splitting SVMTrainData into a class that holds the full training data, and one that holds a subset of the data using references, perhaps named SVMTrainSet. Have an easy way to generate the latter from the first (like how a std::string_view or a std::span works), and use the latter exclusively in the rest of the code.

Missing error checking

Errors seeking or reading from the input stream are ignored by your code. There's a only a check for X being empty and rows not all having the same length. The seeking is unnecessary since you don't seem to be using the variable len at all (ensure you enable compiler warnings, which should have caught this). After the while-loop, check if is.eof() is true, if not there was an error while reading.

Be consistent reporting errors

I see both assert() and throw being using in the constructor that reads from a std::istream. Be consistent, and use throw for any recoverable error.

I also would use std::runtime_error instead of std::invalid_argument, the latter is normally meant when you pass an invalid argument to a function, not for I/O errors.

Avoid unnecessary moving

You don't need to create X first and then move it into X_. In the constructor that scales existing samples, just take X by const reference, and have the constructor make the copy:

SVMTrainData(const Samples &X, const std::vector<float> &y) :
        X_(X), y_(y), n_(y_.size()), C_(X[0].size()) {
    scale();
}

In the constructor that reads from the input stream, you can initialize X_ to the right type and then use it directly inside the constructor body:

explicit SVMTrainData(std::istream& is): X_(Samples{}) {
    ...
    X_.push_back(std::move(row));
    ...
}

Similarly, you can avoid moving the row by adding an entry to X_ and getting a reference to it:

auto &row = X_.emplace_back();

However, that needs some other change to your code to make it possible:

Read lines using std::getline()

The member function getline() of std::istream is a bit unfortunate; it requires you to provide a buffer long enough to hold the line, and then there are all sorts of issues to deal with, like what to do if the line is actually longer than the buffer. A simpler way to read lines is to use std::getline() to read a whole line into a std::string:

std::string buffer;

while (std::getline(is, buffer)) {
    // We now know we have read a full line
    auto &row = X_.emplace_back();
    std::size_t curr{};
    while (...) {
         auto delim = buffer.find(',', curr);
         ...
         row.push_back(result);
    }
}

Don't unnecessarily store duplicate information

There is no need for n_ and C_ in SVMTrainData, just use X_.size() and X_[0].size() directly, create member functions that return X_.size() and X_[0].size(), or use range-for loops to avoid needing to use the sizes explicitly.

Duplicating information always has the issue that you have to keep it in sync, otherwise bad things can happen.

Use range-for loops where appropriate

You can avoid a lot of manual iteration by using range-for loops. This simplifies code and reduces the chance of mistakes, like accidentily swapping n_ and C_. For example, you can rewrite scale() like so:

void scale() {
    feature_bounds_ = ...;

    for (auto &row: X_) {
         for (std::size_t c = 0; c < row.size(); ++c) {
              feature_bounds[c] = {
                  std::min(feature_bounds_[c].first, row[c]),
                  std::max(feature_bounds_[c].second, row[c])
              };
         }
    }

    ...
}

If you wait for C++23 or can use the range-v3 library, then you can use zip() for the inner loop:

using ranges::views::zip;

for (auto &row: X_) {
     for (auto [value, bounds]: zip(row, feature_bounds_) {
          bounds = {
              std::min(bounds.first, value),
              std::max(bounds.second, value)
          };
     }
}

Use of virtual classes

I think most of the use of virtual in your code does not make sense. For the Kernel, I think it would be better to store it as a std::function inside SVMModel:

class SVMModel {
    ...
    std::function<float(const std::vector<float>&, const std::vector<float>&)> kernel_;
    ...
};

The only issue is how to handle adjustParams. Perhaps you can just pass f as a third argument to the kernel function.

For the SVMModel, there are lots of virtual final functions in the base class. That doesn't make sense, why not make those regular functions? As for fit(), this also doesn't need to be virtual, as the only users of the base class's fit() are the derived classes themselves.

Another option is to use template parameters instead of virtual classes. Consider being able to write this in main():

SVMCSVC<KernelRBF> sm;
...
sm.fit(data, params);
...
for (auto y_pred: sm.predict(X1)) {
    std::cout << y_pred << ' ';
}
std::cout << '\n';

Use your own type aliases consistently

You define the type alias Samples, but in most of the code you are writing std::vector<std::vector<float>> explicitly. If you do define a type alias, use it consistently.

Consider going for a more functional style

Your classes are written such that you have an object, and then you call fit() on that object, and now the object is modified to hold the result of the fit. A more functional approach would be to have fit() not modify its own object, but rather return the result of the fit. In fact, fit() could just be a free function, so you can write:

fc::SVMTrainData data(X, y);
fc::SVMParams params;
params.shrinking = false;
params.kernel = KernelRBF;

auto sm = fc::SVMFit(data, params);

The same goes for the Solver class. Instead of first constructing the class with all the parameters, and then calling solver.Solve(), it would be much more natural to write:

auto solution = SolveOrdinary(l, Q, QD, ...);
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.