SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
RationalApproximation.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation