SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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

SHOGUN Machine Learning Toolbox - Documentation