diff --git a/src/shogun/mathematics/linalg/linsolver/CGMShiftedFamilySolver.cpp b/src/shogun/mathematics/linalg/linsolver/CGMShiftedFamilySolver.cpp index b5ff482c38f..774602ac670 100644 --- a/src/shogun/mathematics/linalg/linsolver/CGMShiftedFamilySolver.cpp +++ b/src/shogun/mathematics/linalg/linsolver/CGMShiftedFamilySolver.cpp @@ -47,7 +47,7 @@ SGVector CCGMShiftedFamilySolver::solve( SGVector CCGMShiftedFamilySolver::solve_shifted_weighted( CLinearOperator* A, SGVector b, - SGVector shifts, SGVector weights) + SGVector shifts, SGVector weights, bool negate) { SG_DEBUG("Entering\n"); @@ -133,8 +133,9 @@ SGVector CCGMShiftedFamilySolver::solve_shifted_weighted( float64_t beta=-r_norm2/p_dot_Ap; // compute the zeta-shifted parameter of CG_M - compute_zeta_sh_new(zeta_sh_old, zeta_sh_cur, shifts, beta_old, beta, - alpha, zeta_sh_new); + compute_zeta_sh_new( + zeta_sh_old, zeta_sh_cur, shifts, beta_old, beta, alpha, + zeta_sh_new, negate); // compute beta-shifted parameter of CG_M compute_beta_sh(zeta_sh_new, zeta_sh_cur, beta, beta_sh); diff --git a/src/shogun/mathematics/linalg/linsolver/CGMShiftedFamilySolver.h b/src/shogun/mathematics/linalg/linsolver/CGMShiftedFamilySolver.h index b0e99236d20..d5e6243cfce 100644 --- a/src/shogun/mathematics/linalg/linsolver/CGMShiftedFamilySolver.h +++ b/src/shogun/mathematics/linalg/linsolver/CGMShiftedFamilySolver.h @@ -62,7 +62,8 @@ class CCGMShiftedFamilySolver */ virtual SGVector solve_shifted_weighted( CLinearOperator* A, SGVector b, - SGVector shifts, SGVector weights); + SGVector shifts, SGVector weights, + bool negate = false); /** @return object name */ virtual const char* get_name() const diff --git a/src/shogun/mathematics/linalg/linsolver/IterativeShiftedLinearFamilySolver.cpp b/src/shogun/mathematics/linalg/linsolver/IterativeShiftedLinearFamilySolver.cpp index 8138a0652c8..b14f404f64b 100644 --- a/src/shogun/mathematics/linalg/linsolver/IterativeShiftedLinearFamilySolver.cpp +++ b/src/shogun/mathematics/linalg/linsolver/IterativeShiftedLinearFamilySolver.cpp @@ -30,19 +30,24 @@ CIterativeShiftedLinearFamilySolver::~CIterativeShiftedLinearFamilySolver { } -template -void CIterativeShiftedLinearFamilySolver::compute_zeta_sh_new( - const SGVector& zeta_sh_old, const SGVector& zeta_sh_cur, const SGVector& shifts, - const T& beta_old, const T& beta_cur, const T& alpha, SGVector& zeta_sh_new) + template + void CIterativeShiftedLinearFamilySolver::compute_zeta_sh_new( + const SGVector& zeta_sh_old, const SGVector& zeta_sh_cur, + const SGVector& shifts, const T& beta_old, const T& beta_cur, + const T& alpha, SGVector& zeta_sh_new, bool negate) { // compute zeta_sh_new according to Jergerlehner, eq. 2.44 // [see IterativeShiftedLinearFamilySolver.h] for (index_t i=0; i(0.0)) diff --git a/src/shogun/mathematics/linalg/linsolver/IterativeShiftedLinearFamilySolver.h b/src/shogun/mathematics/linalg/linsolver/IterativeShiftedLinearFamilySolver.h index e6d4a7f4df6..1dc37f6287e 100644 --- a/src/shogun/mathematics/linalg/linsolver/IterativeShiftedLinearFamilySolver.h +++ b/src/shogun/mathematics/linalg/linsolver/IterativeShiftedLinearFamilySolver.h @@ -62,8 +62,9 @@ template class CIterativeShiftedLinearFamilySolver : public * @param weights the weights to be multiplied with each solution for each * shift */ - virtual SGVector solve_shifted_weighted(CLinearOperator* A, - SGVector b, SGVector shifts, SGVector weights) = 0; + virtual SGVector solve_shifted_weighted( + CLinearOperator* A, SGVector b, SGVector shifts, + SGVector weights, bool negate) = 0; /** @return object name */ virtual const char* get_name() const @@ -86,9 +87,10 @@ template class CIterativeShiftedLinearFamilySolver : public * @param alpha \f$\alpha\f$ non-shifted * @param zeta_sh_new \f$\zeta^{\sigma}_{n+1}\f$ to be computed */ - void compute_zeta_sh_new(const SGVector& zeta_sh_old, - const SGVector& zeta_sh_cur, const SGVector& shifts, - const T& beta_old, const T& beta_cur, const T& alpha, SGVector& zeta_sh_new); + void compute_zeta_sh_new( + const SGVector& zeta_sh_old, const SGVector& zeta_sh_cur, + const SGVector& shifts, const T& beta_old, const T& beta_cur, + const T& alpha, SGVector& zeta_sh_new, bool negate = false); /** * compute \f$\beta^{\sigma}_{n}\f$ as \f$\beta_{n}\frac{\zeta^{\sigma}_{n+1}} diff --git a/src/shogun/mathematics/linalg/ratapprox/logdet/LogDetEstimator.cpp b/src/shogun/mathematics/linalg/ratapprox/logdet/LogDetEstimator.cpp index 01f02aa3e61..57e42b574b8 100644 --- a/src/shogun/mathematics/linalg/ratapprox/logdet/LogDetEstimator.cpp +++ b/src/shogun/mathematics/linalg/ratapprox/logdet/LogDetEstimator.cpp @@ -125,9 +125,12 @@ SGVector CLogDetEstimator::sample(index_t num_estimates) index_t num_trace_samples=m_trace_sampler->get_num_samples(); SGVector samples(num_estimates); samples.zero(); + float64_t result = 0.0; +#pragma omp parallel for reduction(+ : result) for (index_t i = 0; i < num_estimates; ++i) { + result = 0.0; for (index_t j = 0; j < num_trace_samples; ++j) { SG_INFO( @@ -135,13 +138,11 @@ SGVector CLogDetEstimator::sample(index_t num_estimates) num_trace_samples); // get the trace sampler vector SGVector s = m_trace_sampler->sample(j); - // calculate the result for sample s - float64_t result = m_operator_log->compute(s); - { - samples[i] += result; - } - } + // calculate the result for sample s and add it to previous + result += m_operator_log->compute(s); } + samples[i] = result; + } SG_INFO("Finished computing %d log-det estimates\n", num_estimates); @@ -165,6 +166,7 @@ SGMatrix CLogDetEstimator::sample_without_averaging( index_t num_trace_samples = m_trace_sampler->get_num_samples(); SGMatrix samples(num_trace_samples, num_estimates); +#pragma omp parallel for for (index_t i = 0; i < num_estimates; ++i) { for (index_t j = 0; j < num_trace_samples; ++j) diff --git a/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/DenseMatrixExactLog.cpp b/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/DenseMatrixExactLog.cpp index 108495eb262..9a9ad6e55eb 100644 --- a/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/DenseMatrixExactLog.cpp +++ b/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/DenseMatrixExactLog.cpp @@ -75,7 +75,7 @@ void CDenseMatrixExactLog::precompute() } #endif // EIGEN_VERSION_AT_LEAST(3,1,0) -float64_t CDenseMatrixExactLog::compute(SGVector sample) +float64_t CDenseMatrixExactLog::compute(SGVector sample) const { SG_DEBUG("Entering...\n"); diff --git a/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/DenseMatrixExactLog.h b/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/DenseMatrixExactLog.h index a5a9ed1f5fc..038c9626fff 100644 --- a/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/DenseMatrixExactLog.h +++ b/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/DenseMatrixExactLog.h @@ -46,7 +46,7 @@ class CDenseMatrixExactLog : public COperatorFunction /** * method that solves the result for a sample */ - virtual float64_t compute(SGVector sample); + virtual float64_t compute(SGVector sample) const; /** @return object name */ virtual const char* get_name() const diff --git a/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationCGM.cpp b/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationCGM.cpp index f68968b82db..673047d060c 100644 --- a/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationCGM.cpp +++ b/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationCGM.cpp @@ -42,9 +42,6 @@ void CLogRationalApproximationCGM::init() SG_ADD((CSGObject**)&m_linear_solver, "linear_solver", "Linear solver for complex systems", MS_NOT_AVAILABLE); - - SG_ADD(&m_negated_shifts, "negated_shifts", - "Negated shifts", MS_NOT_AVAILABLE); } CLogRationalApproximationCGM::~CLogRationalApproximationCGM() @@ -52,24 +49,17 @@ CLogRationalApproximationCGM::~CLogRationalApproximationCGM() SG_UNREF(m_linear_solver); } -float64_t CLogRationalApproximationCGM::compute(SGVector sample) +float64_t +CLogRationalApproximationCGM::compute(SGVector sample) const { SG_DEBUG("Entering\n"); REQUIRE(sample.vector, "Sample is not initialized!\n"); REQUIRE(m_linear_operator, "Operator is not initialized!\n"); - // we need to take the negation of the shifts for this case - if (m_negated_shifts.vector == NULL) - { - m_negated_shifts = SGVector(m_shifts.vlen); - Map shifts(m_shifts.vector, m_shifts.vlen); - Map negated_shifts( - m_negated_shifts.vector, m_negated_shifts.vlen); - negated_shifts = -shifts; - } - + // we need to take the negation of the shifts for this case hence we set + // negate to true SGVector vec = m_linear_solver->solve_shifted_weighted( - m_linear_operator, sample, m_negated_shifts, m_weights); + m_linear_operator, sample, m_shifts, m_weights, true); // Take negative (see CRationalApproximation for the formula) Map v(vec.vector, vec.vlen); diff --git a/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationCGM.h b/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationCGM.h index 2ca5450755b..a55029c38c0 100644 --- a/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationCGM.h +++ b/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationCGM.h @@ -53,7 +53,7 @@ class CLogRationalApproximationCGM : public CRationalApproximation /** * method that solves the result for a sample */ - virtual float64_t compute(SGVector sample); + virtual float64_t compute(SGVector sample) const; /** @return object name */ virtual const char* get_name() const @@ -65,9 +65,6 @@ class CLogRationalApproximationCGM : public CRationalApproximation /** the linear solver for solving complex systems */ CCGMShiftedFamilySolver* m_linear_solver; - /** negated shifts to pass to CG-M linear solver */ - SGVector m_negated_shifts; - /** initialize with default values and register params */ void init(); }; diff --git a/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationIndividual.cpp b/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationIndividual.cpp index 3b982be7ff4..4786f3ac13c 100644 --- a/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationIndividual.cpp +++ b/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationIndividual.cpp @@ -53,7 +53,7 @@ CLogRationalApproximationIndividual::~CLogRationalApproximationIndividual() } float64_t -CLogRationalApproximationIndividual::compute(SGVector sample) +CLogRationalApproximationIndividual::compute(SGVector sample) const { SG_DEBUG("Entering..\n"); REQUIRE(sample.vector, "Sample is not initialized!\n"); diff --git a/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationIndividual.h b/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationIndividual.h index 5057e7cae8e..7f84f9047ca 100644 --- a/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationIndividual.h +++ b/src/shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationIndividual.h @@ -53,7 +53,7 @@ class CLogRationalApproximationIndividual : public CRationalApproximation /** *method that solves the result for a sample */ - virtual float64_t compute(SGVector sample); + virtual float64_t compute(SGVector sample) const; /** @return object name */ virtual const char* get_name() const diff --git a/src/shogun/mathematics/linalg/ratapprox/opfunc/OperatorFunction.h b/src/shogun/mathematics/linalg/ratapprox/opfunc/OperatorFunction.h index 353a21209ab..6786157f32c 100644 --- a/src/shogun/mathematics/linalg/ratapprox/opfunc/OperatorFunction.h +++ b/src/shogun/mathematics/linalg/ratapprox/opfunc/OperatorFunction.h @@ -79,7 +79,7 @@ template class COperatorFunction : public CSGObject /** * Method that computes for a sample and returns the final result */ - virtual float64_t compute(SGVector sample) = 0; + virtual float64_t compute(SGVector sample) const = 0; /** @return object name */ virtual const char* get_name() const diff --git a/src/shogun/mathematics/linalg/ratapprox/opfunc/RationalApproximation.h b/src/shogun/mathematics/linalg/ratapprox/opfunc/RationalApproximation.h index 4d5860ab99b..357cd3b25a8 100644 --- a/src/shogun/mathematics/linalg/ratapprox/opfunc/RationalApproximation.h +++ b/src/shogun/mathematics/linalg/ratapprox/opfunc/RationalApproximation.h @@ -106,7 +106,7 @@ class CRationalApproximation : public COperatorFunction /** * Method that computes for a particular sample */ - virtual float64_t compute(SGVector sample) = 0; + virtual float64_t compute(SGVector sample) const = 0; /** @return shifts */ SGVector get_shifts() const;