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 */ 00009 00010 #include <shogun/lib/common.h> 00011 00012 #ifdef HAVE_LAPACK 00013 #ifdef HAVE_EIGEN3 00014 00015 #include <shogun/base/Parameter.h> 00016 #include <shogun/mathematics/lapack.h> 00017 #include <shogun/mathematics/eigen3.h> 00018 #include <shogun/mathematics/linalg/linop/LinearOperator.h> 00019 #include <shogun/mathematics/linalg/linsolver/IterativeSolverIterator.h> 00020 #include <shogun/mathematics/linalg/eigsolver/LanczosEigenSolver.h> 00021 #include <vector> 00022 00023 using namespace Eigen; 00024 00025 namespace shogun 00026 { 00027 00028 CLanczosEigenSolver::CLanczosEigenSolver() 00029 : CEigenSolver() 00030 { 00031 init(); 00032 } 00033 00034 CLanczosEigenSolver::CLanczosEigenSolver( 00035 CLinearOperator<float64_t>* linear_operator) 00036 : CEigenSolver(linear_operator) 00037 { 00038 init(); 00039 } 00040 00041 void CLanczosEigenSolver::init() 00042 { 00043 m_max_iteration_limit=1000; 00044 m_relative_tolerence=1E-6; 00045 m_absolute_tolerence=1E-6; 00046 00047 SG_ADD(&m_max_iteration_limit, "max_iteration_limit", 00048 "Maximum number of iteration for the solver", MS_NOT_AVAILABLE); 00049 00050 SG_ADD(&m_relative_tolerence, "relative_tolerence", 00051 "Relative tolerence of solver", MS_NOT_AVAILABLE); 00052 00053 SG_ADD(&m_absolute_tolerence, "absolute_tolerence", 00054 "Absolute tolerence of solver", MS_NOT_AVAILABLE); 00055 } 00056 00057 CLanczosEigenSolver::~CLanczosEigenSolver() 00058 { 00059 } 00060 00061 void CLanczosEigenSolver::compute() 00062 { 00063 SG_DEBUG("Entering\n"); 00064 00065 if (m_is_computed_min && m_is_computed_max) 00066 { 00067 SG_DEBUG("Minimum/maximum eigenvalues are already computed, exiting\n"); 00068 return; 00069 } 00070 00071 REQUIRE(m_linear_operator, "Operator is NULL!\n"); 00072 00073 // vector v_0 00074 VectorXd v=VectorXd::Zero(m_linear_operator->get_dimension()); 00075 00076 // vector v_i, for i=1 this is random valued with norm 1 00077 SGVector<float64_t> v_(v.rows()); 00078 Map<VectorXd> v_i(v_.vector, v_.vlen); 00079 v_i=VectorXd::Random(v_.vlen); 00080 v_i/=v_i.norm(); 00081 00082 // the iterator for this iterative solver 00083 IterativeSolverIterator<float64_t> it(v_i, m_max_iteration_limit, 00084 m_relative_tolerence, m_absolute_tolerence); 00085 00086 // the diagonal for Lanczos T-matrix (tridiagonal) 00087 std::vector<float64_t> alpha; 00088 00089 // the subdiagonal for Lanczos T-matrix (tridiagonal) 00090 std::vector<float64_t> beta; 00091 00092 float64_t beta_i=0.0; 00093 SGVector<float64_t> w_(v_.vlen); 00094 00095 // CG iteration begins 00096 for (it.begin(v_i); !it.end(v_i); ++it) 00097 { 00098 SG_DEBUG("CG iteration %d, residual norm %f\n", 00099 it.get_iter_info().iteration_count, 00100 it.get_iter_info().residual_norm); 00101 00102 // apply linear operator to the direction vector 00103 w_=m_linear_operator->apply(v_); 00104 Map<VectorXd> w_i(w_.vector, w_.vlen); 00105 00106 // compute v^{T}Av, if zero, failure 00107 float64_t alpha_i=w_i.dot(v_i); 00108 if (alpha_i==0.0) 00109 break; 00110 00111 // update w_i, v_(i-1) and find beta 00112 w_i-=alpha_i*v_i+beta_i*v; 00113 beta_i=w_i.norm(); 00114 v=v_i; 00115 v_i=w_i/beta_i; 00116 00117 // prepate Lanczos T-matrix from alpha and beta 00118 alpha.push_back(alpha_i); 00119 beta.push_back(beta_i); 00120 } 00121 00122 // solve Lanczos T-matrix to get the eigenvalues 00123 int32_t M=0; 00124 SGVector<float64_t> w(alpha.size()); 00125 int32_t info=0; 00126 00127 // keeping copies of the diagonal and subdiagonal 00128 // because subsequent call to dstemr destroys it 00129 std::vector<float64_t> alpha_orig=alpha; 00130 std::vector<float64_t> beta_orig=beta; 00131 00132 if (!m_is_computed_min) 00133 { 00134 // computing min eigenvalue 00135 wrap_dstemr('N', 'I', alpha.size(), &alpha[0], &beta[0], 00136 0.0, 0.0, 1, 1, &M, w.vector, NULL, 1, 1, NULL, 0.0, &info); 00137 00138 if (info==0) 00139 { 00140 SG_INFO("Iteration took %ld times, residual norm=%.20lf\n", 00141 it.get_iter_info().iteration_count, it.get_iter_info().residual_norm); 00142 00143 m_min_eigenvalue=w[0]; 00144 m_is_computed_min=true; 00145 } 00146 else 00147 SG_WARNING("Some error occured while computing eigenvalues!\n"); 00148 } 00149 00150 if (!m_is_computed_max) 00151 { 00152 // computing max eigenvalue 00153 wrap_dstemr('N', 'I', alpha_orig.size(), &alpha_orig[0], &beta_orig[0], 00154 0.0, 0.0, w.vlen, w.vlen, &M, w.vector, NULL, 1, 1, NULL, 0.0, &info); 00155 00156 if (info==0) 00157 { 00158 SG_INFO("Iteration took %ld times, residual norm=%.20lf\n", 00159 it.get_iter_info().iteration_count, it.get_iter_info().residual_norm); 00160 m_max_eigenvalue=w[0]; 00161 m_is_computed_max=true; 00162 } 00163 else 00164 SG_WARNING("Some error occured while computing eigenvalues!\n"); 00165 } 00166 00167 SG_DEBUG("Leaving\n"); 00168 } 00169 00170 } 00171 #endif // HAVE_EIGEN3 00172 #endif // HAVE_LAPACK