problems/sparse-logistic-regression.cpp

This is an example that builds a problem in C++ that can be loaded dynamically by the alpaqa solvers.

This is an example that builds a problem in C++ that can be loaded dynamically by the alpaqa solvers.The problem is a simple sparse logistic regression task.

All problem functions are implemented as member functions of the Problem class. This class also contains a member funcs of type alpaqa_problem_functions_t, which is used to expose those functions to the solvers. The provided problem functions are added to this funcs struct by the Problem constructor. To adapt the member functions to ordinary function pointers, the alpaqa::member_caller function is used.
The constructor also loads the actual classification data from the given CSV file.

The main entry point is the alpaqa_problem_register function at the bottom. It is the only function that is exported from the shared module, and it is called by the alpaqa::dl::DLProblem class when loading the problem for the solvers.
Its purpose is to create a problem instance, and return a pointer to it, as well as additional information such as how to clean up the created problem, and a pointer to the provided functions (the funcs member from earlier).

The entry point also accepts a pointer to type-erased user data. You can use this to pass any additional data to the problem constructor, such as problem parameters. When using the alpaqa-driver program, the type of this user data is always std::any *. Inside of the std::any, there is a span of string views containing the problem-specific command line options.
In this example, the datafile and λ_factor options are supported. The former selects the CSV file to load the data from, the latter controls the regularization parameter.

To run this example, first build the alpaqa-driver program, and then build the sparse-logistic-regression module. The regression data is generated by executing the generate-data.py script (requires the scikit-learn package), it outputs a CSV file that is loaded in the load_data function defined below.
Finally, execute the driver with the appropriate arguments as shown below.

export CFLAGS="-march=native"; export CXXFLAGS="-march=native"; export FFLAGS="-march=native"
cmake -Bbuild -S. -G "Ninja Multi-Config" -D CMAKE_POSITION_INDEPENDENT_CODE=On
cmake --build build -j --config RelWithDebInfo -t driver
cmake --build build -j --config RelWithDebInfo -t sparse-logistic-regression
python3 -m pip install -U scikit-learn
python3 examples/problems/generate-data.py
./build/src/RelWithDebInfo/alpaqa-driver \
    build/examples/problems/RelWithDebInfo/sparse-logistic-regression.so \
    problem.datafile=breast_cancer.csv \
    problem.λ_factor=1e-5 \
    sol=output/ \
    method=panoc.lbfgs

The solution is written to a CSV file in the output directory.

Run ./build/src/RelWithDebInfo/alpaqa-driver without arguments for usage and argument information.

#include <sparse-logistic-regression-export.h>
 
#include <alpaqa/config/config.hpp>
#include <alpaqa/dl/dl-problem.h>
#include <alpaqa/params/params.hpp>
#include <alpaqa/util/io/csv.hpp>
USING_ALPAQA_CONFIG(alpaqa::DefaultConfig);
 
#include <algorithm>
#include <any>
#include <cassert>
#include <cstdint>
#include <filesystem>
#include <fstream>
#include <iostream>
#include <random>
#include <span>
#include <stdexcept>
#include <string_view>
#include <utility>
namespace fs = std::filesystem;
 
struct Problem {
    alpaqa_problem_functions_t funcs{};
    length_t n;         ///< Number of features
    length_t m;         ///< Number of data points
    real_t λ;           ///< Regularization factor
    real_t μ;           ///< Scaling factor
    mat A;              ///< Data matrix (m×n)
    vec b;              ///< Binary labels (m)
    vec Aᵀb;            ///< Work vector (n)
    mutable vec Ax;     ///< Work vector (m)
    fs::path data_file; ///< File we loaded the data from
 
    /// φ(x) = ∑ ln(1 + exp(-b x))
    real_t logistic_loss(crvec x) const {
        auto &&xa = x.array();
        auto &&ba = b.array();
        return ((-ba * xa).exp() + 1).log().sum();
    }
 
    /// -dφ(x)/dx = b / (exp(b x) + 1)
    void neg_deriv_logistic_loss(crvec x, rvec g) const {
        auto &&xa = x.array();
        auto &&ba = b.array();
        g         = ba / ((ba * xa).exp() + 1);
    }
 
    /// σ₂(x) = b exp(b x) / (exp(b x) + 1)²
    /// @note Assumes that b is a binary vector (i.e. b² = b).
    void sigmoid2(crvec x, rvec g) const {
        auto &&xa = x.array();
        auto &&ba = b.array();
        g         = (ba * xa).exp();
        g         = ba * g.array() / ((g.array() + 1) * (g.array() + 1));
    }
 
    /// Objective function.
    real_t eval_f(const real_t *x_) const {
        cmvec x{x_, n};
        Ax.noalias() = A * x;
        return μ * logistic_loss(Ax);
    }
 
    /// Gradient of objective.
    void eval_grad_f(const real_t *x_, real_t *g_) const {
        cmvec x{x_, n};
        mvec g{g_, n};
        Ax.noalias() = A * x;
        // ∇(f∘A)(x) = Aᵀ∇f(Ax)
        neg_deriv_logistic_loss(Ax, Ax);         // Ax ← -∇f(Ax)
        g.noalias() = -μ * (A.transpose() * Ax); // g ← μAᵀ∇f(Ax)
    }
 
    /// Hessian-vector product of objective.
    void eval_hess_f_prod(const real_t *x_, const real_t *v_,
                          real_t *Hv_) const {
        cmvec x{x_, n};
        cmvec v{v_, n};
        mvec Hv{Hv_, n};
        Ax.noalias() = A * x;
        sigmoid2(Ax, Ax);
        Ax.noalias() = Ax.asDiagonal() * (A * v);
        Hv.noalias() = μ * (A.transpose() * Ax);
    }
 
    /// Hessian of objective.
    void eval_hess_f(const real_t *x_, real_t *H_) const {
        cmvec x{x_, n};
        mmat H{H_, n, n};
        Ax.noalias() = A * x;
        sigmoid2(Ax, Ax);
        H.noalias() = A.transpose() * (μ * Ax.asDiagonal()) * A;
    }
 
    /// Hessian-vector product of Lagrangian.
    void eval_hess_L_prod(const real_t *x, [[maybe_unused]] const real_t *y,
                          real_t scale, const real_t *v, real_t *Hv) const {
        eval_hess_f_prod(x, v, Hv);
        if (scale != 1)
            mvec{Hv, n} *= scale;
    }
 
    /// Hessian-vector product of augmented Lagrangian.
    void eval_hess_ψ_prod(const real_t *x, const real_t *y,
                          [[maybe_unused]] const real_t *Σ, real_t scale,
                          [[maybe_unused]] const real_t *zl,
                          [[maybe_unused]] const real_t *zu, const real_t *v,
                          real_t *Hv) const {
        eval_hess_L_prod(x, y, scale, v, Hv);
    }
 
    /// Hessian of Lagrangian.
    void eval_hess_L(const real_t *x, [[maybe_unused]] const real_t *y,
                     real_t scale, real_t *H) const {
        eval_hess_f(x, H);
        if (scale != 1)
            mmat{H, n, n} *= scale;
    }
 
    /// Hessian of augmented Lagrangian.
    void eval_hess_ψ(const real_t *x, const real_t *y,
                     [[maybe_unused]] const real_t *Σ, real_t scale,
                     [[maybe_unused]] const real_t *zl,
                     [[maybe_unused]] const real_t *zu, real_t *H) const {
        eval_hess_L(x, y, scale, H);
    }
 
    /// Both the objective and its gradient.
    real_t eval_f_grad_f(const real_t *x_, real_t *g_) const {
        cmvec x{x_, n};
        mvec g{g_, n};
        Ax.noalias() = A * x;
        real_t f     = μ * logistic_loss(Ax);
        // ∇(f∘A)(x) = Aᵀ∇f(Ax)
        neg_deriv_logistic_loss(Ax, Ax);         // Ax ← -∇f(Ax)
        g.noalias() = -μ * (A.transpose() * Ax); // g ← μAᵀ∇f(Ax)
        return f;
    }
 
    /// Constraints function (unconstrained).
    void eval_g(const real_t *, real_t *) const {}
 
    /// Gradient-vector product of constraints.
    void eval_grad_g_prod(const real_t *, const real_t *, real_t *gr_) const {
        mvec{gr_, n}.setZero();
    }
 
    /// Jacobian of constraints.
    void eval_jac_g(const real_t *, real_t *) const {}
 
    /// ℓ₁-regularization term.
    void initialize_l1_reg(real_t *lambda, length_t *size) const {
        if (!lambda) {
            *size = 1;
        } else {
            assert(*size == 1);
            *lambda = λ;
        }
    }
 
    /// Loads classification data from a CSV file.
    /// The first row contains the number of data points and the number of
    /// features, separated by a space.
    /// The second row contains the binary labels for all data points.
    /// Every following row contains the values of one feature for all data
    /// points.
    void load_data() {
        std::ifstream csv_file{data_file};
        if (!csv_file)
            throw std::runtime_error("Unable to open file '" +
                                     data_file.string() + "'");
        // Load dimensions (#data points, #features)
        csv_file >> m >> n;
        csv_file.ignore(1, '\n');
        if (!csv_file)
            throw std::runtime_error(
                "Unable to read dimensions from data file");
        b.resize(m);
        A.resize(m, n);
        Aᵀb.resize(n);
        Ax.resize(m);
        // Read the target labels
        alpaqa::csv::read_row(csv_file, b);
        // Read the data
        for (length_t i = 0; i < n; ++i)
            alpaqa::csv::read_row(csv_file, A.col(i));
    }
 
    std::string get_name() const {
        return "sparse logistic regression (\"" + data_file.string() + "\")";
    }
 
    /// Constructor loads CSV data file and exposes the problem functions by
    /// initializing the @c funcs member.
    Problem(fs::path csv_filename, real_t λ_factor)
        : data_file(std::move(csv_filename)) {
        load_data();
        Aᵀb.noalias() = A.transpose() * b;
        real_t λ_max = Aᵀb.lpNorm<Eigen::Infinity>() / static_cast<real_t>(m);
        λ            = λ_factor * λ_max;
        μ            = 1. / static_cast<real_t>(m);
 
        using P = Problem;
        using alpaqa::member_caller;
        funcs.n                = n;
        funcs.m                = 0;
        funcs.eval_f           = member_caller<&P::eval_f>();
        funcs.eval_grad_f      = member_caller<&P::eval_grad_f>();
        funcs.eval_f_grad_f    = member_caller<&P::eval_f_grad_f>();
        funcs.eval_g           = member_caller<&P::eval_g>();
        funcs.eval_grad_g_prod = member_caller<&P::eval_grad_g_prod>();
        funcs.eval_jac_g       = member_caller<&P::eval_jac_g>();
        funcs.eval_hess_L_prod = member_caller<&P::eval_hess_L_prod>();
        funcs.eval_hess_ψ_prod = member_caller<&P::eval_hess_ψ_prod>();
        funcs.eval_hess_L      = member_caller<&P::eval_hess_L>();
        funcs.eval_hess_ψ      = member_caller<&P::eval_hess_ψ>();
        if (λ > 0)
            funcs.initialize_l1_reg = member_caller<&P::initialize_l1_reg>();
    }
};
 
/// Main entry point of this file, it is called by the
/// @ref alpaqa::dl::DLProblem class.
extern "C" SPARSE_LOGISTIC_REGRESSION_EXPORT alpaqa_problem_register_t
alpaqa_problem_register(void *user_data_v) noexcept try {
    // Check and convert user arguments
    if (!user_data_v)
        throw std::invalid_argument("Missing user data");
    const auto &user_data = *reinterpret_cast<std::any *>(user_data_v);
    using param_t         = std::span<std::string_view>;
    auto opts             = std::any_cast<param_t>(user_data);
    std::vector<unsigned> used(opts.size());
    // CSV file to load dataset from
    std::string_view datafilename;
    alpaqa::params::set_params(datafilename, "datafile", opts, used);
    if (datafilename.empty())
        throw std::invalid_argument("Missing option problem.datafile");
    // Regularization factor
    real_t λ_factor = 0.1;
    alpaqa::params::set_params(λ_factor, "λ_factor", opts, used);
    // Check any unused options
    auto unused_opt = std::find(used.begin(), used.end(), 0);
    auto unused_idx = static_cast<size_t>(unused_opt - used.begin());
    if (unused_opt != used.end())
        throw std::invalid_argument("Unused problem option: " +
                                    std::string(opts[unused_idx]));
    // Build and expose problem
    auto problem = std::make_unique<Problem>(datafilename, λ_factor);
    alpaqa_problem_register_t result;
    alpaqa::register_member_function(result, "get_name", &Problem::get_name);
    result.functions = &problem->funcs;
    result.instance  = problem.release();
    result.cleanup   = [](void *instance) {
        delete reinterpret_cast<Problem *>(instance);
    };
    return result;
} catch (...) {
    return {.exception = new alpaqa_exception_ptr_t{std::current_exception()}};
}