SHOGUN  6.1.3
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
13 #include <shogun/lib/SGVector.h>
14 #include <shogun/io/SGIO.h>
15 #include <shogun/lib/Time.h>
17
20
21 using namespace Eigen;
22
23 namespace shogun
24 {
25
28 {
29  SG_GCDEBUG("%s created (%p)\n", this->get_name(), this);
30 }
31
33  : CIterativeLinearSolver<float64_t>(store_residuals)
34 {
35  SG_GCDEBUG("%s created (%p)\n", this->get_name(), this);
36 }
37
39 {
40  SG_GCDEBUG("%s destroyed (%p)\n", this->get_name(), this);
41 }
42
45 {
47
48  // sanity check
49  REQUIRE(A, "Operator is NULL!\n");
50  REQUIRE(A->get_dimension()==b.vlen, "Dimension mismatch!\n");
51
52  // the final solution vector, initial guess is 0
53  SGVector<float64_t> result(b.vlen);
54  result.set_const(0.0);
55
56  // the rest of the part hinges on eigen3 for computing norms
57  Map<VectorXd> x(result.vector, result.vlen);
58  Map<VectorXd> b_map(b.vector, b.vlen);
59
60  // direction vector
61  SGVector<float64_t> p_(result.vlen);
62  Map<VectorXd> p(p_.vector, p_.vlen);
63
64  // residual r_i=b-Ax_i, here x_0=, so r_0=b
65  VectorXd r=b_map;
66
67  // initial direction is same as residual
68  p=r;
69
70  // the iterator for this iterative solver
73
74  // CG iteration begins
75  float64_t r_norm2=r.dot(r);
76
77  // start the timer
78  CTime time;
79  time.start();
80
81  // set the residuals to zero
84
85  for (it.begin(r); !it.end(r); ++it)
86  {
87  SG_DEBUG("CG iteration %d, residual norm %f\n",
88  it.get_iter_info().iteration_count,
89  it.get_iter_info().residual_norm);
90
92  {
93  m_residuals[it.get_iter_info().iteration_count]
94  =it.get_iter_info().residual_norm;
95  }
96
97  // apply linear operator to the direction vector
98  SGVector<float64_t> Ap_=A->apply(p_);
99  Map<VectorXd> Ap(Ap_.vector, Ap_.vlen);
100
101  // compute p^{T}Ap, if zero, failure
102  float64_t p_dot_Ap=p.dot(Ap);
103  if (p_dot_Ap==0.0)
104  break;
105
106  // compute the alpha parameter of CG
107  float64_t alpha=r_norm2/p_dot_Ap;
108
109  // update the solution vector and residual
110  // x_{i}=x_{i-1}+\alpha_{i}p
111  x+=alpha*p;
112
113  // r_{i}=r_{i-1}-\alpha_{i}p
114  r-=alpha*Ap;
115
116  // compute new ||r||_{2}, if zero, converged
117  float64_t r_norm2_i=r.dot(r);
118  if (r_norm2_i==0.0)
119  break;
120
121  // compute the beta parameter of CG
122  float64_t beta=r_norm2_i/r_norm2;
123
124  // update direction, and ||r||_{2}
125  r_norm2=r_norm2_i;
126  p=r+beta*p;
127  }
128
129  float64_t elapsed=time.cur_time_diff();
130
131  if (!it.succeeded(r))
132  SG_WARNING("Did not converge!\n");
133
134  SG_INFO("Iteration took %ld times, residual norm=%.20lf, time elapsed=%lf\n",
135  it.get_iter_info().iteration_count, it.get_iter_info().residual_norm, elapsed);
136
138  return result;
139 }
140
141 }
Class Time that implements a stopwatch based on either cpu time or wall clock time.
Definition: Time.h:42
#define SG_INFO(...)
Definition: SGIO.h:117
void begin(const VectorXt &residual)
virtual SGVector< float64_t > solve(CLinearOperator< float64_t > *A, SGVector< float64_t > b)
Definition: SGMatrix.h:25
#define REQUIRE(x,...)
Definition: SGIO.h:181
const bool end(const VectorXt &residual)
virtual const char * get_name() const
float64_t cur_time_diff(bool verbose=false)
Definition: Time.cpp:68
template class that is used as an iterator for an iterative linear solver. In the iteration of solvin...
#define SG_GCDEBUG(...)
Definition: SGIO.h:101
abstract template base for all iterative linear solvers such as conjugate gradient (CG) solvers...
double float64_t
Definition: common.h:60
float64_t start(bool verbose=false)
Definition: Time.cpp:59
virtual SGVector< T > apply(SGVector< T > b) const =0
void set_const(T const_elem)
Definition: SGVector.cpp:199
const index_t get_dimension() const
#define SG_DEBUG(...)
Definition: SGIO.h:106
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
const bool succeeded(const VectorXt &residual)
#define SG_WARNING(...)
Definition: SGIO.h:127
index_t vlen
Definition: SGVector.h:571

SHOGUN Machine Learning Toolbox - Documentation