TMB Documentation  v1.9.11
sparse_matrix_exponential.hpp
Go to the documentation of this file.
1 #ifdef TMBAD_FRAMEWORK
2 
13 
14 /* Shared data and methods for sparse matrix vector multiply: A^T * x */
15 template <class Type>
16 struct SparseATx {
17  // Eigens index type
18  typedef typename
19  Eigen::SparseMatrix<Type>::StorageIndex StorageIndex;
20  // Sparse matrix (A) dimension and number of non zeros
21  StorageIndex nrow, ncol, nnz;
22  // Sparsity pattern (A)
23  std::vector<StorageIndex> i; // Length nnz
24  std::vector<StorageIndex> p; // Length ncol+1
25  SparseATx(const Eigen::SparseMatrix<Type> &x) :
26  nrow(x.rows()),
27  ncol(x.cols()),
28  nnz(x.nonZeros()),
29  i(x.innerIndexPtr(), x.innerIndexPtr() + nnz),
30  p(x.outerIndexPtr(), x.outerIndexPtr() + ncol + 1) {}
31  StorageIndex rows() const { return nrow; }
32  StorageIndex cols() const { return ncol; }
33  /* Operation: y = A^T * x */
34  template<class T>
35  void f(const T *A, const T *x, T *y) {
36  for (StorageIndex j = 0; j < ncol; j++) {
37  y[j] = 0;
38  for (StorageIndex k = p[j]; k < p[j+1]; k++) {
39  y[j] += A[k] * x[i[k]]; // good access
40  }
41  }
42  }
43  /* Adjoint operation */
44  template<class T>
45  void df(const T *A, const T *x, const T *y,
46  T *dA, T *dx, const T *dy) {
47  for (StorageIndex j = 0; j < ncol; j++) {
48  for (StorageIndex k = p[j]; k < p[j+1]; k++) {
49  dA[k] += x[i[k]] * dy[j]; // good access
50  dx[i[k]] += A[k] * dy[j]; // bad access
51  }
52  }
53  }
54 };
55 
56 
57 template<class Type, bool transpose=false>
58 struct SpAxOp : TMBad::global::DynamicOperator< -1, -1> {
59  typedef typename
60  Eigen::SparseMatrix<Type>::StorageIndex StorageIndex;
61  std::shared_ptr<SparseATx<Type> > P;
62  SpAxOp() {}
63  SpAxOp(const Eigen::SparseMatrix<Type> &x) :
64  P(std::make_shared<SparseATx<Type> >( x )) { }
65  static const bool have_input_size_output_size = true;
66  TMBad::Index input_size() const {
67  return 2; // Two pointers
68  }
69  TMBad::Index output_size() const {
70  return (*P).cols();
71  }
72  // FIXME:
73  // static const bool forward_replay_copy = true;
74  template<class T>
75  void forward(TMBad::ForwardArgs<T> &args) {
76  const T* A = args.x_ptr(0);
77  const T* x = args.x_ptr(1);
78  T* y = args.y_ptr(0);
79  (*P).template f(A, x, y);
80  }
81  template<class T>
82  void reverse(TMBad::ReverseArgs<T> &args) {
83  const T* A = args.x_ptr(0);
84  const T* x = args.x_ptr(1);
85  const T* y = args.y_ptr(0);
86  T* dA = args.dx_ptr(0);
87  T* dx = args.dx_ptr(1);
88  const T* dy = args.dy_ptr(0);
89  (*P).template df(A, x, y, dA, dx, dy);
90  }
91  // ---- Dependencies ---- (copied from ad_blas)
92  // FIXME: Make general implementation for pointer based case
93  void dependencies(TMBad::Args<> &args, TMBad::Dependencies &dep) const {
94  dep.add_segment(args.input(0), (*P).nnz );
95  dep.add_segment(args.input(1), (*P).rows() );
96  }
97  static const bool have_dependencies = true;
99  static const bool implicit_dependencies = true;
101  static const bool allow_remap = false;
102  // Not yet implemented
103  void forward(TMBad::ForwardArgs<TMBad::Writer> &args) { ASSERT(false); }
104  void reverse(TMBad::ReverseArgs<TMBad::Writer> &args) { ASSERT(false); }
105  const char* op_name() { return "SpAxOp"; }
106 };
107 
108 // FIXME: Remember to transpose A on construction
109 
111 template<class T>
112 struct config {
114  bool normalize;
116  bool trace;
118  bool warn;
120  int Nmax;
121  config() : normalize(false), trace(true), warn(true), Nmax(100) {}
122 };
123 
134 template <class T>
135 struct expm_series {
136  typedef vectorize::vector<T> vec;
138  T N;
140  vec A_values;
148  void update(Eigen::SparseMatrix<T> &A) {
149  // FIXME: Assert same pattern
150  A_values = vec(A.valuePtr(), A.nonZeros());
151  }
152  expm_series() {}
153  expm_series(Eigen::SparseMatrix<T> &A, T N, config<T> cfg=config<T>()) :
154  N(N), A_values(A.valuePtr(), A.nonZeros()), multiply(A), cfg(cfg)
155  { }
157  vec operator()(vec x) {
158  N = min(N, T(cfg.Nmax));
159  std::vector<TMBad::ad_segment> args = {A_values, x, vec(N, 1)};
160  if (! F.initialized() ) {
161  struct Test {
162  config<T> cfg;
163  TMBad::Scalar Nold;
164  bool operator() (const std::vector<TMBad::Scalar*> &x) {
165  using TMBad::operator<<;
166  TMBad::Scalar N = x[2][0];
167  if ( (int) N == cfg.Nmax) {
168  if (cfg.warn)
169  Rf_warning("expm: N terms reduced to Nmax (%i)", (int) cfg.Nmax);
170  }
171  bool change = (Nold != N);
172  if (cfg.trace && change) {
173  Rcout << "Retaping:" << " Nold=" << Nold << " Nnew=" << N << "\n";
174  Nold = N;
175  }
176  return change;
177  }
178  };
179  Test N_changed = {cfg, N.Value()};
180  F = TMBad::ADFun_retaping(*this, args, N_changed);
181  }
182  return F(args);
183  }
184 private:
185  friend class TMBad::PackWrap<expm_series>;
186  // Packed: (A, x) -> exp(A) * x
187  TMBad::ad_segment operator() (const std::vector<TMBad::ad_segment> &args) {
188  // Unpack input
189  vec A = args[0];
190  vec x = args[1];
191  vec N_= args[2];
192  int N = (int) N_[0].Value();
193  // Evaluate series
194  vec term(x), y(x);
195  for (int n=1; n<N; n++) {
196  term = multiply(A, term) / n;
197  y += term;
198  }
199  return y;
200  }
201 };
202 
203 template<>
204 struct expm_series<double> {
206  int N;
207  Eigen::SparseMatrix<double> A;
208  config<double> cfg;
209  void update(Eigen::SparseMatrix<double> &A) {
210  this->A = A;
211  }
212  expm_series() {}
213  expm_series(const Eigen::SparseMatrix<double> &A, double N, config<double> cfg=config<double>()) : N(N), A(A), cfg(cfg)
214  {
215  if ((int) N > cfg.Nmax) {
216  if (cfg.warn)
217  Rf_warning("expm: N terms reduced to Nmax (%i)", (int) cfg.Nmax);
218  this->N = cfg.Nmax;
219  }
220  }
221  vec operator()(vec x) {
222  // Evaluate series
223  vec term(x), y(x);
224  for (int n=1; n<N; n++) {
225  term = A * term / n;
226  y += term;
227  }
228  return y;
229  }
230 };
231 
241 template<class T>
243  typedef vectorize::vector<T> vec;
244  expm_series<T> ExpA;
245  T rho;
249  T getN(T rho) {
250  return ceil(rho + 4*sqrt(rho) + 5);
251  }
252  // FIXME: Calculate rho=abs(max(diag(Q)))
253  expm_generator(Eigen::SparseMatrix<T> &Q, config<T> cfg=config<T>()) {
254  using TMBad::min;
255  Eigen::SparseMatrix<T> A(Q);
256  vector<T> d = Q.diagonal();
257  if (d.size() > 0) {
258  T M = d[0];
259  for (int i=1; i<d.size(); i++) {
260  M = min(M, d[i]);
261  }
262  rho = -M;
263  }
264  A.diagonal().array() += rho;
265  ExpA = expm_series<T>(A, getN(rho), cfg);
266  }
268  vec operator()(vec x) {
269  vec y = ExpA(x);
270  y = exp(-rho) * y;
271  if (ExpA.cfg.normalize) {
272  y = y / sum(y);
273  }
274  return y;
275  }
276 };
277 
278 } // End namespace
279 
280 #endif
Vector class used by TMB.
Definition: vector.hpp:17
Matrix exponential of generator matrix multiplied by vector.
Type * y_ptr(Index j)
pointer version - use with caution.
Definition: global.hpp:330
Type df(Type x, Type df1, Type df2, int give_log)
Probability density function of the Fisher distribution.
Access input/output values and derivatives during a reverse pass. Write access granted for the input ...
Definition: global.hpp:311
Shared configuration parameters for expm_series and expm_generator
Contiguous set of variables on the current tape.
Definition: global.hpp:2780
Type * y_ptr(Index j)
pointer version - use with caution.
Definition: global.hpp:291
int Nmax
Use no more than this number of terms in series.
Access input/output values during a forward pass. Write access granted for the output value only...
Definition: global.hpp:279
Transform a functor to have packed input/output.
Definition: checkpoint.hpp:253
Type * x_ptr(Index j)
pointer version - use with caution.
Definition: global.hpp:328
Vector of consequtive values on the AD tape.
Type min(const vector< Type > &x)
Type * x_ptr(Index j)
pointer version - use with caution.
Definition: global.hpp:289
bool normalize
Normalize result? (expm_generator only)
Type * dx_ptr(Index j)
pointer version - use with caution.
Definition: global.hpp:332
Sparse matrix exponentiation.
Operator auto-completion.
Definition: global.hpp:2129
Operator that requires dynamic allocation. Compile time known input/output size.
Definition: global.hpp:1590
std::vector< T > operator()(const std::vector< T > &xp)
Transformed functor assuming original maps std::vector<ad_segment> to ad_segment
Definition: checkpoint.hpp:259
Type sum(Vector< Type > x)
Definition: convenience.hpp:33
Matrix exponential multiplied by vector.
ADFun ADFun_retaping(Functor &F, const std::vector< ad_aug > &x, Test test=Test())
Construct ADFun that automatically retapes.
Definition: TMBad.hpp:1165
Index input(Index j) const
Get variable index of j&#39;th input to current operator.
Definition: global.hpp:265
TMBad::global::Complete< SpAxOp< T > > multiply
Type * dy_ptr(Index j)
pointer version - use with caution.
Definition: global.hpp:334
T getN(T rho)
Uniformization Grassmann (1977) eq. 10. => error < 1e-4 rho=abs(max(diag(Q)))
Container of ADFun object with packed input and output.
Definition: TMBad.hpp:1174
License: GPL v2