parrallel computation of sample in LogDetEstimator (#4235)

* computing sample_without_averaging in parrallel * compute CLogDetEstimator::sample in parallel * thread safety of compute * remove extra vector from LogRationalApproximationCGM::compute
shogun-toolbox · May 21, 2018 · 41755fb · 41755fb
1 parent ccf4f84
commit 41755fb
Show file tree

Hide file tree

Showing 13 changed files with 44 additions and 46 deletions.
diff --git a/src/shogun/mathematics/linalg/linsolver/CGMShiftedFamilySolver.cpp b/src/shogun/mathematics/linalg/linsolver/CGMShiftedFamilySolver.cpp
@@ -47,7 +47,7 @@ SGVector<float64_t> CCGMShiftedFamilySolver::solve(
 
 SGVector<complex128_t> CCGMShiftedFamilySolver::solve_shifted_weighted(
 	CLinearOperator<float64_t>* A, SGVector<float64_t> b,
-	SGVector<complex128_t> shifts, SGVector<complex128_t> weights)
+	SGVector<complex128_t> shifts, SGVector<complex128_t> weights, bool negate)
 {
 	SG_DEBUG("Entering\n");
 
@@ -133,8 +133,9 @@ SGVector<complex128_t> 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
@@ -62,7 +62,8 @@ class CCGMShiftedFamilySolver
 	 */
 	virtual SGVector<complex128_t> solve_shifted_weighted(
 		CLinearOperator<float64_t>* A, SGVector<float64_t> b,
-		SGVector<complex128_t> shifts, SGVector<complex128_t> weights);
+		SGVector<complex128_t> shifts, SGVector<complex128_t> 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
@@ -30,19 +30,24 @@ CIterativeShiftedLinearFamilySolver<T, ST>::~CIterativeShiftedLinearFamilySolver
 	{
 	}
 
-template <class T, class ST>
-void CIterativeShiftedLinearFamilySolver<T, ST>::compute_zeta_sh_new(
-		const SGVector<ST>& zeta_sh_old, const SGVector<ST>& zeta_sh_cur, const SGVector<ST>& shifts,
-		const T& beta_old, const T& beta_cur, const T& alpha, SGVector<ST>& zeta_sh_new)
+	template <class T, class ST>
+	void CIterativeShiftedLinearFamilySolver<T, ST>::compute_zeta_sh_new(
+	    const SGVector<ST>& zeta_sh_old, const SGVector<ST>& zeta_sh_cur,
+	    const SGVector<ST>& shifts, const T& beta_old, const T& beta_cur,
+	    const T& alpha, SGVector<ST>& zeta_sh_new, bool negate)
 	{
 		// compute zeta_sh_new according to Jergerlehner, eq. 2.44
 		// [see IterativeShiftedLinearFamilySolver.h]
 		for (index_t i=0; i<zeta_sh_new.vlen; ++i)
 		{
+			ST shift_value = shifts[i];
+			if (negate)
+				shift_value = -shifts[i];
 			ST numer=zeta_sh_old[i]*zeta_sh_cur[i]*beta_old;
 
-			ST denom=beta_cur*alpha*(zeta_sh_old[i]-zeta_sh_cur[i])
-				+beta_old*zeta_sh_old[i]*(1.0-beta_cur*shifts[i]);
+			ST denom =
+			    beta_cur * alpha * (zeta_sh_old[i] - zeta_sh_cur[i]) +
+			    beta_old * zeta_sh_old[i] * (1.0 - beta_cur * shift_value);
 
 			// handle division by zero
 			if (denom==static_cast<ST>(0.0))

diff --git a/src/shogun/mathematics/linalg/linsolver/IterativeShiftedLinearFamilySolver.h b/src/shogun/mathematics/linalg/linsolver/IterativeShiftedLinearFamilySolver.h
@@ -62,8 +62,9 @@ template<class T, class ST=T> class CIterativeShiftedLinearFamilySolver : public
 	 * @param weights the weights to be multiplied with each solution for each
 	 * shift
 	 */
-	virtual SGVector<ST> solve_shifted_weighted(CLinearOperator<T>* A,
-		SGVector<T> b, SGVector<ST> shifts, SGVector<ST> weights) = 0;
+	virtual SGVector<ST> solve_shifted_weighted(
+		CLinearOperator<T>* A, SGVector<T> b, SGVector<ST> shifts,
+		SGVector<ST> weights, bool negate) = 0;
 
 	/** @return object name */
 	virtual const char* get_name() const
@@ -86,9 +87,10 @@ template<class T, class ST=T> 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<ST>& zeta_sh_old,
-		const SGVector<ST>& zeta_sh_cur, const SGVector<ST>& shifts,
-		const T& beta_old, const T& beta_cur, const T& alpha, SGVector<ST>& zeta_sh_new);
+	void compute_zeta_sh_new(
+		const SGVector<ST>& zeta_sh_old, const SGVector<ST>& zeta_sh_cur,
+		const SGVector<ST>& shifts, const T& beta_old, const T& beta_cur,
+		const T& alpha, SGVector<ST>& 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
@@ -125,23 +125,24 @@ SGVector<float64_t> CLogDetEstimator::sample(index_t num_estimates)
 	index_t num_trace_samples=m_trace_sampler->get_num_samples();
 	SGVector<float64_t> 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(
 				"Computing log-determinant trace sample %d/%d\n", j,
 				num_trace_samples);
 			// get the trace sampler vector
 			SGVector<float64_t> 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<float64_t> CLogDetEstimator::sample_without_averaging(
 	index_t num_trace_samples = m_trace_sampler->get_num_samples();
 	SGMatrix<float64_t> 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
@@ -75,7 +75,7 @@ void CDenseMatrixExactLog::precompute()
 }
 #endif // EIGEN_VERSION_AT_LEAST(3,1,0)
 
-float64_t CDenseMatrixExactLog::compute(SGVector<float64_t> sample)
+float64_t CDenseMatrixExactLog::compute(SGVector<float64_t> 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
@@ -46,7 +46,7 @@ class CDenseMatrixExactLog : public COperatorFunction<float64_t>
 	/**
 	 * method that solves the result for a sample
 	 */
-	virtual float64_t compute(SGVector<float64_t> sample);
+	virtual float64_t compute(SGVector<float64_t> 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
@@ -42,34 +42,24 @@ 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()
 {
 	SG_UNREF(m_linear_solver);
 }
 
-float64_t CLogRationalApproximationCGM::compute(SGVector<float64_t> sample)
+float64_t
+CLogRationalApproximationCGM::compute(SGVector<float64_t> 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<complex128_t>(m_shifts.vlen);
-		Map<VectorXcd> shifts(m_shifts.vector, m_shifts.vlen);
-		Map<VectorXcd> 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<complex128_t> 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<VectorXcd> 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
@@ -53,7 +53,7 @@ class CLogRationalApproximationCGM : public CRationalApproximation
 	/**
 	 * method that solves the result for a sample
 	 */
-	virtual float64_t compute(SGVector<float64_t> sample);
+	virtual float64_t compute(SGVector<float64_t> 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<complex128_t> 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
@@ -53,7 +53,7 @@ CLogRationalApproximationIndividual::~CLogRationalApproximationIndividual()
 }
 
 float64_t
-CLogRationalApproximationIndividual::compute(SGVector<float64_t> sample)
+CLogRationalApproximationIndividual::compute(SGVector<float64_t> 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
@@ -53,7 +53,7 @@ class CLogRationalApproximationIndividual : public CRationalApproximation
 	/**
 	 *method that solves the result for a sample
 	 */
-	virtual float64_t compute(SGVector<float64_t> sample);
+	virtual float64_t compute(SGVector<float64_t> 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
@@ -79,7 +79,7 @@ template<class T> class COperatorFunction : public CSGObject
 	/**
 	 * Method that computes for a sample and returns the final result
 	 */
-	virtual float64_t compute(SGVector<T> sample) = 0;
+	virtual float64_t compute(SGVector<T> 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
@@ -106,7 +106,7 @@ class CRationalApproximation : public COperatorFunction<float64_t>
 	/**
 	 * Method that computes for a particular sample
 	 */
-	virtual float64_t compute(SGVector<float64_t> sample) = 0;
+	virtual float64_t compute(SGVector<float64_t> sample) const = 0;
 
 	/** @return shifts */
 	SGVector<complex128_t> get_shifts() const;