kernel_matrix_solver.h

Go to the documentation of this file.

/*
 * Copyright 2024 MusicScience37 (Kenta Kabashima)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once
 
#include <cmath>
#include <limits>
 
#include <Eigen/Core>
#include <Eigen/Eigenvalues>
#include <Eigen/IterativeLinearSolvers>
#include <Eigen/LU>
#include <Eigen/SparseCore>
#include <Eigen/src/Core/util/Constants.h>
 
#include "num_collect/base/concepts/dense_vector.h"
#include "num_collect/base/exception.h"
#include "num_collect/base/precondition.h"
#include "num_collect/logging/logging_macros.h"
#include "num_collect/rbf/kernel_matrix_type.h"
 
namespace num_collect::rbf::impl {

template <typename KernelValue, typename FunctionValue,
    kernel_matrix_type KernelMatrixType, bool UsesGlobalLengthParameter>
class kernel_matrix_solver;

template <typename KernelValue, typename FunctionValue>
class kernel_matrix_solver<KernelValue, FunctionValue,
    kernel_matrix_type::dense, true> {
public:
    using kernel_matrix_type = Eigen::MatrixX<KernelValue>;

    using vector_type = Eigen::VectorX<FunctionValue>;

    using scalar_type = KernelValue;

    kernel_matrix_solver() = default;

    void compute(
        const kernel_matrix_type& kernel_matrix, const vector_type& data) {
        decomposition_.compute(kernel_matrix, Eigen::ComputeEigenvectors);
        spectre_ = decomposition_.eigenvectors().adjoint() * data;
    }

    void solve(vector_type& coefficients, scalar_type reg_param,
        const vector_type& data) const {
        (void)data;  // not used in this version.
 
        reg_param = correct_reg_param_if_needed(reg_param);
 
        coefficients = decomposition_.eigenvectors() *
            (decomposition_.eigenvalues().array() + reg_param)
                .inverse()
                .matrix()
                .asDiagonal() *
            spectre_;
    }

    [[nodiscard]] auto calc_mle_objective(scalar_type reg_param) const
        -> scalar_type {
        reg_param = correct_reg_param_if_needed(reg_param);
 
        constexpr scalar_type limit = std::numeric_limits<scalar_type>::max() *
            static_cast<scalar_type>(1e-20);
        if (decomposition_.eigenvalues()(0) + reg_param <=
            static_cast<scalar_type>(0)) {
            return limit;
        }
 
        using std::log;
        const scalar_type value = static_cast<scalar_type>(spectre_.rows()) *
                log(calc_reg_term(reg_param)) +
            calc_log_determinant(reg_param);
        if (value < limit) {
            return value;
        }
        return limit;
    }

    [[nodiscard]] auto calc_common_coeff(scalar_type reg_param) const
        -> scalar_type {
        reg_param = correct_reg_param_if_needed(reg_param);
 
        return calc_reg_term(reg_param) /
            static_cast<scalar_type>(spectre_.rows());
    }

    template <base::concepts::dense_vector InputData>
    [[nodiscard]] auto calc_reg_term(
        const InputData& data, scalar_type reg_param) const -> scalar_type {
        reg_param = correct_reg_param_if_needed(reg_param);
 
        return ((decomposition_.eigenvectors().adjoint() * data)
                    .array()
                    .abs2()
                    .rowwise()
                    .sum() /
            (decomposition_.eigenvalues().array() + reg_param))
            .sum();
    }

    [[nodiscard]] auto eigenvalues() const -> decltype(auto) {
        return decomposition_.eigenvalues();
    }
 
private:
    [[nodiscard]] auto calc_reg_term(const scalar_type& reg_param) const
        -> scalar_type {
        return (spectre_.array().abs2().rowwise().sum() /
            (decomposition_.eigenvalues().array() + reg_param))
            .sum();
    }

    [[nodiscard]] auto calc_log_determinant(const scalar_type& reg_param) const
        -> scalar_type {
        return (decomposition_.eigenvalues().array() + reg_param).log().sum();
    }

    [[nodiscard]] auto correct_reg_param_if_needed(
        const scalar_type& reg_param) const noexcept -> scalar_type {
        const scalar_type smallest_eigenvalue = decomposition_.eigenvalues()(0);
        const scalar_type largest_eigenvalue = decomposition_.eigenvalues()(
            decomposition_.eigenvalues().size() - 1);
        const scalar_type eigenvalue_safe_limit =
            largest_eigenvalue * std::numeric_limits<scalar_type>::epsilon();
        const scalar_type reg_param_safe_limit =
            eigenvalue_safe_limit - smallest_eigenvalue;
 
        if (reg_param < reg_param_safe_limit) {
            return reg_param_safe_limit;
        }
        return reg_param;
    }

    Eigen::SelfAdjointEigenSolver<kernel_matrix_type> decomposition_;

    vector_type spectre_;
};

template <typename KernelValue, typename FunctionValue>
class kernel_matrix_solver<KernelValue, FunctionValue,
    kernel_matrix_type::dense, false> {
public:
    using kernel_matrix_type = Eigen::MatrixX<KernelValue>;

    using vector_type = Eigen::VectorX<FunctionValue>;

    using scalar_type = KernelValue;

    kernel_matrix_solver() = default;

    void compute(
        const kernel_matrix_type& kernel_matrix, const vector_type& data) {
        (void)data;  // not used in this version.
 
        solver_.compute(kernel_matrix);
    }

    void solve(vector_type& coefficients, scalar_type reg_param,
        const vector_type& data) const {
        NUM_COLLECT_PRECONDITION(reg_param == static_cast<scalar_type>(0),
            "Non-zero regularization parameter cannot be used in this "
            "configuration.");
 
        coefficients = solver_.solve(data);
    }
 
private:
    Eigen::PartialPivLU<kernel_matrix_type> solver_;
};

template <typename KernelValue, typename FunctionValue>
class kernel_matrix_solver<KernelValue, FunctionValue,
    kernel_matrix_type::sparse, false> {
public:
    using kernel_matrix_type =
        Eigen::SparseMatrix<KernelValue, Eigen::RowMajor>;

    using vector_type = Eigen::VectorX<FunctionValue>;

    using scalar_type = KernelValue;

    kernel_matrix_solver() = default;

    void compute(
        const kernel_matrix_type& kernel_matrix, const vector_type& data) {
        (void)data;  // not used in this version.
 
        solver_.compute(kernel_matrix);
    }

    void solve(vector_type& coefficients, scalar_type reg_param,
        const vector_type& data) const {
        NUM_COLLECT_PRECONDITION(reg_param == static_cast<scalar_type>(0),
            "Non-zero regularization parameter cannot be used in this "
            "configuration.");
 
        coefficients = solver_.solve(data);
    }
 
private:
    Eigen::BiCGSTAB<kernel_matrix_type> solver_;
};
 
}  // namespace num_collect::rbf::impl