SHOGUN  v3.0.0
LanczosEigenSolver.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
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>
18
21 #include <vector>
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
48  "Maximum number of iteration for the solver", MS_NOT_AVAILABLE);
49
51  "Relative tolerence of solver", MS_NOT_AVAILABLE);
52
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_EIGEN3
172 #endif // HAVE_LAPACK

SHOGUN Machine Learning Toolbox - Documentation