SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
LanczosEigenSolver.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  */
9 
10 #include <shogun/lib/common.h>
11 
12 #ifdef HAVE_LAPACK
13 #ifdef HAVE_EIGEN3
14 
15 #include <shogun/base/Parameter.h>
21 #include <vector>
22 #include <shogun/lib/SGVector.h>
23 
24 using namespace Eigen;
25 
26 namespace shogun
27 {
28 
29 CLanczosEigenSolver::CLanczosEigenSolver()
30  : CEigenSolver()
31 {
32  init();
33 }
34 
37  : CEigenSolver(linear_operator)
38 {
39  init();
40 }
41 
42 void CLanczosEigenSolver::init()
43 {
44  m_max_iteration_limit=1000;
45  m_relative_tolerence=1E-6;
46  m_absolute_tolerence=1E-6;
47 
48  SG_ADD(&m_max_iteration_limit, "max_iteration_limit",
49  "Maximum number of iteration for the solver", MS_NOT_AVAILABLE);
50 
51  SG_ADD(&m_relative_tolerence, "relative_tolerence",
52  "Relative tolerence of solver", MS_NOT_AVAILABLE);
53 
54  SG_ADD(&m_absolute_tolerence, "absolute_tolerence",
55  "Absolute tolerence of solver", MS_NOT_AVAILABLE);
56 }
57 
59 {
60 }
61 
63 {
64  SG_DEBUG("Entering\n");
65 
67  {
68  SG_DEBUG("Minimum/maximum eigenvalues are already computed, exiting\n");
69  return;
70  }
71 
72  REQUIRE(m_linear_operator, "Operator is NULL!\n");
73 
74  // vector v_0
75  VectorXd v=VectorXd::Zero(m_linear_operator->get_dimension());
76 
77  // vector v_i, for i=1 this is random valued with norm 1
78  SGVector<float64_t> v_(v.rows());
79  Map<VectorXd> v_i(v_.vector, v_.vlen);
80  v_i=VectorXd::Random(v_.vlen);
81  v_i/=v_i.norm();
82 
83  // the iterator for this iterative solver
84  IterativeSolverIterator<float64_t> it(v_i, m_max_iteration_limit,
85  m_relative_tolerence, m_absolute_tolerence);
86 
87  // the diagonal for Lanczos T-matrix (tridiagonal)
88  std::vector<float64_t> alpha;
89 
90  // the subdiagonal for Lanczos T-matrix (tridiagonal)
91  std::vector<float64_t> beta;
92 
93  float64_t beta_i=0.0;
94  SGVector<float64_t> w_(v_.vlen);
95 
96  // CG iteration begins
97  for (it.begin(v_i); !it.end(v_i); ++it)
98  {
99  SG_DEBUG("CG iteration %d, residual norm %f\n",
100  it.get_iter_info().iteration_count,
101  it.get_iter_info().residual_norm);
102 
103  // apply linear operator to the direction vector
104  w_=m_linear_operator->apply(v_);
105  Map<VectorXd> w_i(w_.vector, w_.vlen);
106 
107  // compute v^{T}Av, if zero, failure
108  float64_t alpha_i=w_i.dot(v_i);
109  if (alpha_i==0.0)
110  break;
111 
112  // update w_i, v_(i-1) and find beta
113  w_i-=alpha_i*v_i+beta_i*v;
114  beta_i=w_i.norm();
115  v=v_i;
116  v_i=w_i/beta_i;
117 
118  // prepate Lanczos T-matrix from alpha and beta
119  alpha.push_back(alpha_i);
120  beta.push_back(beta_i);
121  }
122 
123  // solve Lanczos T-matrix to get the eigenvalues
124  int32_t M=0;
125  SGVector<float64_t> w(alpha.size());
126  int32_t info=0;
127 
128  // keeping copies of the diagonal and subdiagonal
129  // because subsequent call to dstemr destroys it
130  std::vector<float64_t> alpha_orig=alpha;
131  std::vector<float64_t> beta_orig=beta;
132 
133  if (!m_is_computed_min)
134  {
135  // computing min eigenvalue
136  wrap_dstemr('N', 'I', alpha.size(), &alpha[0], &beta[0],
137  0.0, 0.0, 1, 1, &M, w.vector, NULL, 1, 1, NULL, 0.0, &info);
138 
139  if (info==0)
140  {
141  SG_INFO("Iteration took %ld times, residual norm=%.20lf\n",
142  it.get_iter_info().iteration_count, it.get_iter_info().residual_norm);
143 
144  m_min_eigenvalue=w[0];
145  m_is_computed_min=true;
146  }
147  else
148  SG_WARNING("Some error occured while computing eigenvalues!\n");
149  }
150 
151  if (!m_is_computed_max)
152  {
153  // computing max eigenvalue
154  wrap_dstemr('N', 'I', alpha_orig.size(), &alpha_orig[0], &beta_orig[0],
155  0.0, 0.0, w.vlen, w.vlen, &M, w.vector, NULL, 1, 1, NULL, 0.0, &info);
156 
157  if (info==0)
158  {
159  SG_INFO("Iteration took %ld times, residual norm=%.20lf\n",
160  it.get_iter_info().iteration_count, it.get_iter_info().residual_norm);
161  m_max_eigenvalue=w[0];
162  m_is_computed_max=true;
163  }
164  else
165  SG_WARNING("Some error occured while computing eigenvalues!\n");
166  }
167 
168  SG_DEBUG("Leaving\n");
169 }
170 
171 }
172 #endif // HAVE_EIGEN3
173 #endif // HAVE_LAPACK
#define SG_INFO(...)
Definition: SGIO.h:118
void begin(const VectorXt &residual)
float64_t m_min_eigenvalue
Definition: EigenSolver.h:93
Definition: SGMatrix.h:20
float64_t m_max_eigenvalue
Definition: EigenSolver.h:96
#define REQUIRE(x,...)
Definition: SGIO.h:206
const bool end(const VectorXt &residual)
void wrap_dstemr(char jobz, char range, int n, double *diag, double *subdiag, double vl, double vu, int il, int iu, int *m, double *w, double *z__, int ldz, int nzc, int *isuppz, int tryrac, int *info)
Definition: lapack.cpp:373
template class that is used as an iterator for an iterative linear solver. In the iteration of solvin...
double float64_t
Definition: common.h:50
#define SG_DEBUG(...)
Definition: SGIO.h:107
Abstract base class that provides an abstract compute method for computing eigenvalues of a real valu...
Definition: EigenSolver.h:24
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
Abstract template base class that represents a linear operator, e.g. a matrix.
CLinearOperator< SGVector< float64_t >, SGVector< float64_t > > * m_linear_operator
Definition: EigenSolver.h:99
#define SG_WARNING(...)
Definition: SGIO.h:128
#define SG_ADD(...)
Definition: SGObject.h:81

SHOGUN Machine Learning Toolbox - Documentation