SHOGUN
v3.2.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2013 Soumyajit De 00008 * Written (W) 2013 Heiko Strathmann 00009 */ 00010 00011 #include <shogun/lib/config.h> 00012 #include <shogun/lib/SGVector.h> 00013 #include <shogun/base/Parameter.h> 00014 #include <shogun/mathematics/Math.h> 00015 #include <shogun/mathematics/JacobiEllipticFunctions.h> 00016 #include <shogun/mathematics/linalg/linop/LinearOperator.h> 00017 #include <shogun/mathematics/linalg/linsolver/LinearSolver.h> 00018 #include <shogun/mathematics/linalg/eigsolver/EigenSolver.h> 00019 #include <shogun/mathematics/linalg/ratapprox/opfunc/RationalApproximation.h> 00020 #include <shogun/lib/computation/engine/IndependentComputationEngine.h> 00021 00022 namespace shogun 00023 { 00024 00025 CRationalApproximation::CRationalApproximation() 00026 : COperatorFunction<float64_t>() 00027 { 00028 init(); 00029 00030 SG_GCDEBUG("%s created (%p)\n", this->get_name(), this) 00031 } 00032 00033 CRationalApproximation::CRationalApproximation( 00034 CLinearOperator<float64_t>* linear_operator, 00035 CIndependentComputationEngine* computation_engine, 00036 CEigenSolver* eigen_solver, 00037 float64_t desired_accuracy, 00038 EOperatorFunction function_type) 00039 : COperatorFunction<float64_t>(linear_operator, computation_engine, 00040 function_type) 00041 { 00042 init(); 00043 00044 m_eigen_solver=eigen_solver; 00045 SG_REF(m_eigen_solver); 00046 00047 m_desired_accuracy=desired_accuracy; 00048 00049 SG_GCDEBUG("%s created (%p)\n", this->get_name(), this) 00050 } 00051 00052 CRationalApproximation::~CRationalApproximation() 00053 { 00054 SG_UNREF(m_eigen_solver); 00055 00056 SG_GCDEBUG("%s destroyed (%p)\n", this->get_name(), this) 00057 } 00058 00059 void CRationalApproximation::init() 00060 { 00061 m_eigen_solver=NULL; 00062 m_constant_multiplier=0.0; 00063 m_num_shifts=0; 00064 m_desired_accuracy=0.0; 00065 00066 SG_ADD((CSGObject**)&m_eigen_solver, "eigen_solver", 00067 "Eigen solver for computing extremal eigenvalues", MS_NOT_AVAILABLE); 00068 00069 SG_ADD(&m_shifts, "complex_shifts", "Complex shifts in the linear system", 00070 MS_NOT_AVAILABLE); 00071 00072 SG_ADD(&m_weights, "complex_weights", "Complex weights of the linear system", 00073 MS_NOT_AVAILABLE); 00074 00075 SG_ADD(&m_constant_multiplier, "constant_multiplier", 00076 "Constant multiplier in the rational approximation", 00077 MS_NOT_AVAILABLE); 00078 00079 SG_ADD(&m_num_shifts, "num_shifts", 00080 "Number of shifts in the quadrature rule", MS_NOT_AVAILABLE); 00081 00082 SG_ADD(&m_desired_accuracy, "desired_accuracy", 00083 "Desired accuracy of the rational approximation", MS_NOT_AVAILABLE); 00084 } 00085 00086 SGVector<complex128_t> CRationalApproximation::get_shifts() const 00087 { 00088 return m_shifts; 00089 } 00090 00091 SGVector<complex128_t> CRationalApproximation::get_weights() const 00092 { 00093 return m_weights; 00094 } 00095 00096 float64_t CRationalApproximation::get_constant_multiplier() const 00097 { 00098 return m_constant_multiplier; 00099 } 00100 00101 index_t CRationalApproximation::get_num_shifts() const 00102 { 00103 return m_num_shifts; 00104 } 00105 00106 void CRationalApproximation::set_num_shifts(index_t num_shifts) 00107 { 00108 m_num_shifts=num_shifts; 00109 } 00110 00111 void CRationalApproximation::precompute() 00112 { 00113 // compute extremal eigenvalues 00114 m_eigen_solver->compute(); 00115 SG_INFO("max_eig=%.15lf\n", m_eigen_solver->get_max_eigenvalue()); 00116 SG_INFO("min_eig=%.15lf\n", m_eigen_solver->get_min_eigenvalue()); 00117 00118 REQUIRE(m_eigen_solver->get_min_eigenvalue()>0, 00119 "Minimum eigenvalue is negative, please provide a Hermitian matrix\n"); 00120 00121 // compute number of shifts from accuracy if shifts are not set yet 00122 if (m_num_shifts==0) 00123 m_num_shifts=compute_num_shifts_from_accuracy(); 00124 00125 SG_INFO("Computing %d shifts\n", m_num_shifts); 00126 compute_shifts_weights_const(); 00127 } 00128 00129 int32_t CRationalApproximation::compute_num_shifts_from_accuracy() 00130 { 00131 REQUIRE(m_desired_accuracy>0, "Desired accuracy must be positive but is %f\n", 00132 m_desired_accuracy); 00133 00134 float64_t max_eig=m_eigen_solver->get_max_eigenvalue(); 00135 float64_t min_eig=m_eigen_solver->get_min_eigenvalue(); 00136 00137 float64_t log_cond_number=CMath::log(max_eig)-CMath::log(min_eig); 00138 float64_t two_pi_sq=2.0*M_PI*M_PI; 00139 00140 int32_t num_shifts=static_cast<index_t>(-1.5*(log_cond_number+6.0) 00141 *CMath::log(m_desired_accuracy)/two_pi_sq); 00142 00143 return num_shifts; 00144 } 00145 00146 void CRationalApproximation::compute_shifts_weights_const() 00147 { 00148 float64_t PI=M_PI; 00149 00150 // eigenvalues are always real in this case 00151 float64_t max_eig=m_eigen_solver->get_max_eigenvalue(); 00152 float64_t min_eig=m_eigen_solver->get_min_eigenvalue(); 00153 00154 // we need $(\frac{\lambda_{M}}{\lambda_{m}})^{frac{1}{4}}$ and 00155 // $(\lambda_{M}*\lambda_{m})^{frac{1}{4}}$ for the rest of the 00156 // calculation where $lambda_{M}$ and $\lambda_{m} are the maximum 00157 // and minimum eigenvalue respectively 00158 float64_t m_div=CMath::pow(max_eig/min_eig, 0.25); 00159 float64_t m_mult=CMath::pow(max_eig*min_eig, 0.25); 00160 00161 // k=$\frac{(\frac{\lambda_{M}}{\lambda_{m}})^\frac{1}{4}-1} 00162 // {(\frac{\lambda_{M}}{\lambda_{m}})^\frac{1}{4}+1}$ 00163 float64_t k=(m_div-1)/(m_div+1); 00164 float64_t L=-CMath::log(k)/PI; 00165 00166 // compute K and K' 00167 float64_t K=0.0, Kp=0.0; 00168 CJacobiEllipticFunctions::ellipKKp(L, K, Kp); 00169 00170 // compute constant multiplier 00171 m_constant_multiplier=-8*K*m_mult/(k*PI*m_num_shifts); 00172 00173 // compute Jacobi elliptic functions sn, cn, dn and compute shifts, weights 00174 // using conformal mapping of the quadrature rule for discretization of the 00175 // contour integral 00176 float64_t m=CMath::sq(k); 00177 00178 // allocate memory for shifts 00179 m_shifts=SGVector<complex128_t>(m_num_shifts); 00180 m_weights=SGVector<complex128_t>(m_num_shifts); 00181 00182 for (index_t i=0; i<m_num_shifts; ++i) 00183 { 00184 complex128_t t=complex128_t(0.0, 0.5*Kp)-K+(0.5+i)*2*K/m_num_shifts; 00185 00186 complex128_t sn, cn, dn; 00187 CJacobiEllipticFunctions::ellipJC(t, m, sn, cn, dn); 00188 00189 complex128_t w=m_mult*(1.0+k*sn)/(1.0-k*sn); 00190 complex128_t dzdt=cn*dn/CMath::sq(1.0/k-sn); 00191 00192 // compute shifts and weights as per Hale et. al. (2008) [2] 00193 m_shifts[i]=CMath::sq(w); 00194 00195 switch (m_function_type) 00196 { 00197 case OF_SQRT: 00198 m_weights[i]=dzdt; 00199 break; 00200 case OF_LOG: 00201 m_weights[i]=2.0*CMath::log(w)*dzdt/w; 00202 break; 00203 case OF_POLY: 00204 SG_NOTIMPLEMENTED 00205 m_weights[i]=complex128_t(0.0); 00206 break; 00207 case OF_UNDEFINED: 00208 SG_WARNING("Operator function is undefined!\n") 00209 m_weights[i]=complex128_t(0.0); 00210 break; 00211 } 00212 } 00213 } 00214 00215 }