SHOGUN  4.2.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 
14 #include <shogun/base/Parameter.h>
20 #include <vector>
21 #include <shogun/lib/SGVector.h>
22 
23 using namespace Eigen;
24 
25 namespace shogun
26 {
27 
28 CLanczosEigenSolver::CLanczosEigenSolver()
29  : CEigenSolver()
30 {
31  init();
32 }
33 
35  CLinearOperator<float64_t>* linear_operator)
36  : CEigenSolver(linear_operator)
37 {
38  init();
39 }
40 
41 void CLanczosEigenSolver::init()
42 {
43  m_max_iteration_limit=1000;
44  m_relative_tolerence=1E-6;
45  m_absolute_tolerence=1E-6;
46 
47  SG_ADD(&m_max_iteration_limit, "max_iteration_limit",
48  "Maximum number of iteration for the solver", MS_NOT_AVAILABLE);
49 
50  SG_ADD(&m_relative_tolerence, "relative_tolerence",
51  "Relative tolerence of solver", MS_NOT_AVAILABLE);
52 
53  SG_ADD(&m_absolute_tolerence, "absolute_tolerence",
54  "Absolute tolerence of solver", MS_NOT_AVAILABLE);
55 }
56 
58 {
59 }
60 
62 {
63  SG_DEBUG("Entering\n");
64 
66  {
67  SG_DEBUG("Minimum/maximum eigenvalues are already computed, exiting\n");
68  return;
69  }
70 
71  REQUIRE(m_linear_operator, "Operator is NULL!\n");
72 
73  // vector v_0
74  VectorXd v=VectorXd::Zero(m_linear_operator->get_dimension());
75 
76  // vector v_i, for i=1 this is random valued with norm 1
77  SGVector<float64_t> v_(v.rows());
78  Map<VectorXd> v_i(v_.vector, v_.vlen);
79  v_i=VectorXd::Random(v_.vlen);
80  v_i/=v_i.norm();
81 
82  // the iterator for this iterative solver
83  IterativeSolverIterator<float64_t> it(v_i, m_max_iteration_limit,
84  m_relative_tolerence, m_absolute_tolerence);
85 
86  // the diagonal for Lanczos T-matrix (tridiagonal)
87  std::vector<float64_t> alpha;
88 
89  // the subdiagonal for Lanczos T-matrix (tridiagonal)
90  std::vector<float64_t> beta;
91 
92  float64_t beta_i=0.0;
93  SGVector<float64_t> w_(v_.vlen);
94 
95  // CG iteration begins
96  for (it.begin(v_i); !it.end(v_i); ++it)
97  {
98  SG_DEBUG("CG iteration %d, residual norm %f\n",
99  it.get_iter_info().iteration_count,
100  it.get_iter_info().residual_norm);
101 
102  // apply linear operator to the direction vector
103  w_=m_linear_operator->apply(v_);
104  Map<VectorXd> w_i(w_.vector, w_.vlen);
105 
106  // compute v^{T}Av, if zero, failure
107  float64_t alpha_i=w_i.dot(v_i);
108  if (alpha_i==0.0)
109  break;
110 
111  // update w_i, v_(i-1) and find beta
112  w_i-=alpha_i*v_i+beta_i*v;
113  beta_i=w_i.norm();
114  v=v_i;
115  v_i=w_i/beta_i;
116 
117  // prepate Lanczos T-matrix from alpha and beta
118  alpha.push_back(alpha_i);
119  beta.push_back(beta_i);
120  }
121 
122  // solve Lanczos T-matrix to get the eigenvalues
123  int32_t M=0;
124  SGVector<float64_t> w(alpha.size());
125  int32_t info=0;
126 
127  // keeping copies of the diagonal and subdiagonal
128  // because subsequent call to dstemr destroys it
129  std::vector<float64_t> alpha_orig=alpha;
130  std::vector<float64_t> beta_orig=beta;
131 
132  if (!m_is_computed_min)
133  {
134  // computing min eigenvalue
135  wrap_dstemr('N', 'I', alpha.size(), &alpha[0], &beta[0],
136  0.0, 0.0, 1, 1, &M, w.vector, NULL, 1, 1, NULL, 0.0, &info);
137 
138  if (info==0)
139  {
140  SG_INFO("Iteration took %ld times, residual norm=%.20lf\n",
141  it.get_iter_info().iteration_count, it.get_iter_info().residual_norm);
142 
143  m_min_eigenvalue=w[0];
144  m_is_computed_min=true;
145  }
146  else
147  SG_WARNING("Some error occured while computing eigenvalues!\n");
148  }
149 
150  if (!m_is_computed_max)
151  {
152  // computing max eigenvalue
153  wrap_dstemr('N', 'I', alpha_orig.size(), &alpha_orig[0], &beta_orig[0],
154  0.0, 0.0, w.vlen, w.vlen, &M, w.vector, NULL, 1, 1, NULL, 0.0, &info);
155 
156  if (info==0)
157  {
158  SG_INFO("Iteration took %ld times, residual norm=%.20lf\n",
159  it.get_iter_info().iteration_count, it.get_iter_info().residual_norm);
160  m_max_eigenvalue=w[0];
161  m_is_computed_max=true;
162  }
163  else
164  SG_WARNING("Some error occured while computing eigenvalues!\n");
165  }
166 
167  SG_DEBUG("Leaving\n");
168 }
169 
170 }
171 #endif // HAVE_LAPACK
#define SG_INFO(...)
Definition: SGIO.h:118
void begin(const VectorXt &residual)
const index_t get_dimension() const
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
virtual SGVector< T > apply(SGVector< T > b) const =0
#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
#define SG_WARNING(...)
Definition: SGIO.h:128
#define SG_ADD(...)
Definition: SGObject.h:81
CLinearOperator< float64_t > * m_linear_operator
Definition: EigenSolver.h:99

SHOGUN Machine Learning Toolbox - Documentation