SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
LanczosEigenSolver.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  */
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation