SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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>
21 
22 namespace shogun
23 {
24 
27 {
28  init();
29 
30  SG_GCDEBUG("%s created (%p)\n", this->get_name(), this)
31 }
32 
34  CLinearOperator<float64_t>* linear_operator,
35  CIndependentComputationEngine* computation_engine,
36  CEigenSolver* eigen_solver,
37  float64_t desired_accuracy,
38  EOperatorFunction function_type)
39  : COperatorFunction<float64_t>(linear_operator, computation_engine,
40  function_type)
41 {
42  init();
43 
44  m_eigen_solver=eigen_solver;
46 
47  m_desired_accuracy=desired_accuracy;
48 
49  SG_GCDEBUG("%s created (%p)\n", this->get_name(), this)
50 }
51 
53 {
55 
56  SG_GCDEBUG("%s destroyed (%p)\n", this->get_name(), this)
57 }
58 
59 void CRationalApproximation::init()
60 {
61  m_eigen_solver=NULL;
63  m_num_shifts=0;
65 
66  SG_ADD((CSGObject**)&m_eigen_solver, "eigen_solver",
67  "Eigen solver for computing extremal eigenvalues", MS_NOT_AVAILABLE);
68 
69  SG_ADD(&m_shifts, "complex_shifts", "Complex shifts in the linear system",
71 
72  SG_ADD(&m_weights, "complex_weights", "Complex weights of the linear system",
74 
75  SG_ADD(&m_constant_multiplier, "constant_multiplier",
76  "Constant multiplier in the rational approximation",
78 
79  SG_ADD(&m_num_shifts, "num_shifts",
80  "Number of shifts in the quadrature rule", MS_NOT_AVAILABLE);
81 
82  SG_ADD(&m_desired_accuracy, "desired_accuracy",
83  "Desired accuracy of the rational approximation", MS_NOT_AVAILABLE);
84 }
85 
87 {
88  return m_shifts;
89 }
90 
92 {
93  return m_weights;
94 }
95 
97 {
98  return m_constant_multiplier;
99 }
100 
102 {
103  return m_num_shifts;
104 }
105 
107 {
108  m_num_shifts=num_shifts;
109 }
110 
112 {
113  // compute extremal eigenvalues
115  SG_INFO("max_eig=%.15lf\n", m_eigen_solver->get_max_eigenvalue());
116  SG_INFO("min_eig=%.15lf\n", m_eigen_solver->get_min_eigenvalue());
117 
119  "Minimum eigenvalue is negative, please provide a Hermitian matrix\n");
120 
121  // compute number of shifts from accuracy if shifts are not set yet
122  if (m_num_shifts==0)
124 
125  SG_INFO("Computing %d shifts\n", m_num_shifts);
126  compute_shifts_weights_const();
127 }
128 
130 {
131  REQUIRE(m_desired_accuracy>0, "Desired accuracy must be positive but is %f\n",
133 
136 
137  float64_t log_cond_number=CMath::log(max_eig)-CMath::log(min_eig);
138  float64_t two_pi_sq=2.0*M_PI*M_PI;
139 
140  int32_t num_shifts=static_cast<index_t>(-1.5*(log_cond_number+6.0)
141  *CMath::log(m_desired_accuracy)/two_pi_sq);
142 
143  return num_shifts;
144 }
145 
146 void CRationalApproximation::compute_shifts_weights_const()
147 {
148  float64_t PI=M_PI;
149 
150  // eigenvalues are always real in this case
153 
154  // we need $(\frac{\lambda_{M}}{\lambda_{m}})^{frac{1}{4}}$ and
155  // $(\lambda_{M}*\lambda_{m})^{frac{1}{4}}$ for the rest of the
156  // calculation where $lambda_{M}$ and $\lambda_{m} are the maximum
157  // and minimum eigenvalue respectively
158  float64_t m_div=CMath::pow(max_eig/min_eig, 0.25);
159  float64_t m_mult=CMath::pow(max_eig*min_eig, 0.25);
160 
161  // k=$\frac{(\frac{\lambda_{M}}{\lambda_{m}})^\frac{1}{4}-1}
162  // {(\frac{\lambda_{M}}{\lambda_{m}})^\frac{1}{4}+1}$
163  float64_t k=(m_div-1)/(m_div+1);
164  float64_t L=-CMath::log(k)/PI;
165 
166  // compute K and K'
167  float64_t K=0.0, Kp=0.0;
169 
170  // compute constant multiplier
171  m_constant_multiplier=-8*K*m_mult/(k*PI*m_num_shifts);
172 
173  // compute Jacobi elliptic functions sn, cn, dn and compute shifts, weights
174  // using conformal mapping of the quadrature rule for discretization of the
175  // contour integral
176  float64_t m=CMath::sq(k);
177 
178  // allocate memory for shifts
181 
182  for (index_t i=0; i<m_num_shifts; ++i)
183  {
184  complex128_t t=complex128_t(0.0, 0.5*Kp)-K+(0.5+i)*2*K/m_num_shifts;
185 
186  complex128_t sn, cn, dn;
187  CJacobiEllipticFunctions::ellipJC(t, m, sn, cn, dn);
188 
189  complex128_t w=m_mult*(1.0+k*sn)/(1.0-k*sn);
190  complex128_t dzdt=cn*dn/CMath::sq(1.0/k-sn);
191 
192  // compute shifts and weights as per Hale et. al. (2008) [2]
193  m_shifts[i]=CMath::sq(w);
194 
195  switch (m_function_type)
196  {
197  case OF_SQRT:
198  m_weights[i]=dzdt;
199  break;
200  case OF_LOG:
201  m_weights[i]=2.0*CMath::log(w)*dzdt/w;
202  break;
203  case OF_POLY:
205  m_weights[i]=complex128_t(0.0);
206  break;
207  case OF_UNDEFINED:
208  SG_WARNING("Operator function is undefined!\n")
209  m_weights[i]=complex128_t(0.0);
210  break;
211  }
212  }
213 }
214 
215 }
virtual void compute()=0
#define SG_INFO(...)
Definition: SGIO.h:118
std::complex< float64_t > complex128_t
Definition: common.h:67
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.
SGVector< complex128_t > get_weights() const
int32_t index_t
Definition: common.h:62
static T sq(T x)
Definition: Math.h:450
#define REQUIRE(x,...)
Definition: SGIO.h:206
#define SG_NOTIMPLEMENTED
Definition: SGIO.h:139
const float64_t get_min_eigenvalue() const
Definition: EigenSolver.h:68
static void ellipKKp(Real L, Real &K, Real &Kp)
virtual const char * get_name() const
SGVector< complex128_t > m_shifts
#define SG_REF(x)
Definition: SGObject.h:51
static void ellipJC(Complex u, Real m, Complex &sn, Complex &cn, Complex &dn)
#define SG_GCDEBUG(...)
Definition: SGIO.h:102
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:112
double float64_t
Definition: common.h:50
#define M_PI
workaround for log2 being a define on cygwin
Definition: Math.h:59
SGVector< complex128_t > m_weights
void set_num_shifts(index_t num_shifts)
#define SG_UNREF(x)
Definition: SGObject.h:52
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:922
SGVector< complex128_t > get_shifts() const
const float64_t get_max_eigenvalue() const
Definition: EigenSolver.h:81
#define SG_WARNING(...)
Definition: SGIO.h:128
#define SG_ADD(...)
Definition: SGObject.h:81
static int32_t pow(bool x, int32_t n)
Definition: Math.h:535

SHOGUN Machine Learning Toolbox - Documentation