SHOGUN  6.1.3
RationalApproximation.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Soumyajit De
8  * Written (W) 2013 Heiko Strathmann
9  */
10 
11 #include <shogun/lib/config.h>
12 #include <shogun/lib/SGVector.h>
13 #include <shogun/base/Parameter.h>
15 #ifdef USE_GPL_SHOGUN
16 #include <shogun/mathematics/JacobiEllipticFunctions.h>
17 #endif //USE_GPL_SHOGUN
23 
24 namespace shogun
25 {
26 
29 {
30  init();
31 
32  SG_GCDEBUG("%s created (%p)\n", this->get_name(), this)
33 }
34 
36  CLinearOperator<float64_t>* linear_operator,
37  CIndependentComputationEngine* computation_engine,
39  float64_t desired_accuracy,
40  EOperatorFunction function_type)
41  : COperatorFunction<float64_t>(linear_operator, computation_engine,
42  function_type)
43 {
44  init();
45 
48 
49  m_desired_accuracy=desired_accuracy;
50 
51  SG_GCDEBUG("%s created (%p)\n", this->get_name(), this)
52 }
53 
55 {
57 
58  SG_GCDEBUG("%s destroyed (%p)\n", this->get_name(), this)
59 }
60 
61 void CRationalApproximation::init()
62 {
63  m_eigen_solver=NULL;
65  m_num_shifts=0;
67 
68  SG_ADD((CSGObject**)&m_eigen_solver, "eigen_solver",
69  "Eigen solver for computing extremal eigenvalues", MS_NOT_AVAILABLE);
70 
71  SG_ADD(&m_shifts, "complex_shifts", "Complex shifts in the linear system",
73 
74  SG_ADD(&m_weights, "complex_weights", "Complex weights of the linear system",
76 
77  SG_ADD(&m_constant_multiplier, "constant_multiplier",
78  "Constant multiplier in the rational approximation",
80 
81  SG_ADD(&m_num_shifts, "num_shifts",
82  "Number of shifts in the quadrature rule", MS_NOT_AVAILABLE);
83 
84  SG_ADD(&m_desired_accuracy, "desired_accuracy",
85  "Desired accuracy of the rational approximation", MS_NOT_AVAILABLE);
86 }
87 
89 {
90  return m_shifts;
91 }
92 
94 {
95  return m_weights;
96 }
97 
99 {
100  return m_constant_multiplier;
101 }
102 
104 {
105  return m_num_shifts;
106 }
107 
109 {
110  m_num_shifts=num_shifts;
111 }
112 
114 {
115  // compute extremal eigenvalues
117  SG_INFO("max_eig=%.15lf\n", m_eigen_solver->get_max_eigenvalue());
118  SG_INFO("min_eig=%.15lf\n", m_eigen_solver->get_min_eigenvalue());
119 
121  "Minimum eigenvalue is negative, please provide a Hermitian matrix\n");
122 
123  // compute number of shifts from accuracy if shifts are not set yet
124  if (m_num_shifts==0)
126 
127  SG_INFO("Computing %d shifts\n", m_num_shifts);
128  compute_shifts_weights_const();
129 }
130 
132 {
133  REQUIRE(m_desired_accuracy>0, "Desired accuracy must be positive but is %f\n",
135 
138 
139  float64_t log_cond_number=CMath::log(max_eig)-CMath::log(min_eig);
140  float64_t two_pi_sq=2.0*M_PI*M_PI;
141 
142  int32_t num_shifts=static_cast<index_t>(-1.5*(log_cond_number+6.0)
143  *CMath::log(m_desired_accuracy)/two_pi_sq);
144 
145  return num_shifts;
146 }
147 
148 void CRationalApproximation::compute_shifts_weights_const()
149 {
150  float64_t PI=M_PI;
151 
152  // eigenvalues are always real in this case
155 
156  // we need $(\frac{\lambda_{M}}{\lambda_{m}})^{frac{1}{4}}$ and
157  // $(\lambda_{M}*\lambda_{m})^{frac{1}{4}}$ for the rest of the
158  // calculation where $lambda_{M}$ and $\lambda_{m} are the maximum
159  // and minimum eigenvalue respectively
160  float64_t m_div=CMath::pow(max_eig/min_eig, 0.25);
161  float64_t m_mult=CMath::pow(max_eig*min_eig, 0.25);
162 
163  // k=$\frac{(\frac{\lambda_{M}}{\lambda_{m}})^\frac{1}{4}-1}
164  // {(\frac{\lambda_{M}}{\lambda_{m}})^\frac{1}{4}+1}$
165  float64_t k=(m_div-1)/(m_div+1);
166  float64_t L=-CMath::log(k)/PI;
167 
168  // compute K and K'
169  float64_t K=0.0, Kp=0.0;
170 #ifdef USE_GPL_SHOGUN
171  CJacobiEllipticFunctions::ellipKKp(L, K, Kp);
172 #else
174 #endif //USE_GPL_SHOGUN
175 
176  // compute constant multiplier
177  m_constant_multiplier=-8*K*m_mult/(k*PI*m_num_shifts);
178 
179  // compute Jacobi elliptic functions sn, cn, dn and compute shifts, weights
180  // using conformal mapping of the quadrature rule for discretization of the
181  // contour integral
182  float64_t m=CMath::sq(k);
183 
184  // allocate memory for shifts
187 
188  for (index_t i=0; i<m_num_shifts; ++i)
189  {
190  complex128_t t=complex128_t(0.0, 0.5*Kp)-K+(0.5+i)*2*K/m_num_shifts;
191 
192  complex128_t sn, cn, dn;
193 #ifdef USE_GPL_SHOGUN
194  CJacobiEllipticFunctions::ellipJC(t, m, sn, cn, dn);
195 #else
197 #endif //USE_GPL_SHOGUN
198 
199  complex128_t w=m_mult*(1.0+k*sn)/(1.0-k*sn);
200  complex128_t dzdt=cn*dn/CMath::sq(1.0/k-sn);
201 
202  // compute shifts and weights as per Hale et. al. (2008) [2]
203  m_shifts[i]=CMath::sq(w);
204 
205  switch (m_function_type)
206  {
207  case OF_SQRT:
208  m_weights[i]=dzdt;
209  break;
210  case OF_LOG:
211  m_weights[i]=2.0*CMath::log(w)*dzdt/w;
212  break;
213  case OF_POLY:
215  m_weights[i]=complex128_t(0.0);
216  break;
217  case OF_UNDEFINED:
218  SG_WARNING("Operator function is undefined!\n")
219  m_weights[i]=complex128_t(0.0);
220  break;
221  }
222  }
223 }
224 
225 }
virtual void compute()=0
#define SG_INFO(...)
Definition: SGIO.h:117
std::complex< float64_t > complex128_t
Definition: common.h:77
SGVector< complex128_t > get_weights() const
void eigen_solver(const SGMatrix< T > &A, SGVector< T > &eigenvalues, SGMatrix< T > &eigenvectors)
int32_t index_t
Definition: common.h:72
static T sq(T x)
Definition: Math.h:418
#define REQUIRE(x,...)
Definition: SGIO.h:181
#define SG_NOTIMPLEMENTED
Definition: SGIO.h:138
const float64_t get_min_eigenvalue() const
Definition: EigenSolver.h:68
virtual const char * get_name() const
SGVector< complex128_t > m_shifts
#define SG_REF(x)
Definition: SGObject.h:52
#define SG_GCDEBUG(...)
Definition: SGIO.h:101
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:124
#define SG_GPL_ONLY
Definition: SGIO.h:139
Abstract template base class for computing for a linear operator C and a vector s. submit_jobs method creates a bunch of jobs needed to solve for this particular and attaches one unique job aggregator to each of them, then submits them all to the computation engine.
double float64_t
Definition: common.h:60
#define M_PI
Definition: Math.h:50
SGVector< complex128_t > m_weights
void set_num_shifts(index_t num_shifts)
#define SG_UNREF(x)
Definition: SGObject.h:53
Abstract base class that provides an abstract compute method for computing eigenvalues of a real valu...
Definition: EigenSolver.h:24
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
Abstract base class for solving multiple independent instances of CIndependentJob. It has one method, submit_job, which may add the job to an internal queue and might block if there is yet not space in the queue. After jobs are submitted, it might not yet be ready. wait_for_all waits until all jobs are completed, which *must be* called to guarantee that all jobs are finished.
static float64_t log(float64_t v)
Definition: Math.h:714
SGVector< complex128_t > get_shifts() const
const float64_t get_max_eigenvalue() const
Definition: EigenSolver.h:81
#define SG_WARNING(...)
Definition: SGIO.h:127
#define SG_ADD(...)
Definition: SGObject.h:93
static int32_t pow(bool x, int32_t n)
Definition: Math.h:474

SHOGUN Machine Learning Toolbox - Documentation