SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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 
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 }

SHOGUN Machine Learning Toolbox - Documentation