1 /* 2 * autodiff.h 3 * 4 * Created on: 16 Apr 2013 5 * Author: s0965328 6 */ 7 8 #ifndef AUTODIFF_H_ 9 #define AUTODIFF_H_ 10 #include <boost/unordered_set.hpp> 11 #include <boost/numeric/ublas/matrix_proxy.hpp> 12 #include <boost/numeric/ublas/matrix.hpp> 13 #include <boost/numeric/ublas/matrix_sparse.hpp> 14 #include <boost/numeric/ublas/io.hpp> 15 #include "auto_diff_types.h" 16 #include "Node.h" 17 #include "VNode.h" 18 #include "OPNode.h" 19 #include "PNode.h" 20 #include "ActNode.h" 21 #include "EdgeSet.h" 22 23 24 /* 25 * + Function and Gradient Evaluation 26 * The tapeless implementation for function and derivative evaluation 27 * Advantage for tapeless: 28 * Low memory usage 29 * Function evaluation use one stack 30 * Gradient evaluation use two stack. 31 * Disadvantage for tapeless: 32 * Inefficient if the expression tree have repeated nodes. 33 * for example: 34 * root 35 * / \ 36 * * * 37 * / \ / \ 38 * x1 x1 x1 x1 39 * Tapeless implementation will go through all the edges. 40 * ie. adjoint of x will be updated 4 times for the correct 41 * gradient of x1.While the tape implemenation can discovery this 42 * dependence and update adjoint of x1 just twice. The computational 43 * graph (DAG) for a taped implemenation will be look like bellow. 44 * root 45 * /\ 46 * * 47 * /\ 48 * x1 49 * 50 * + Forward Hessian Evaluation: 51 * This is an inefficient implementation of the forward Hessian method. It will evaluate the diagonal 52 * and upper triangular of the Hessian. The gradient is also evaluation in the same routine. The result 53 * will be returned in an array. 54 * To use this method, one have to provide a len parameter. len = (nvar+3)*nvar/2 where nvar is the number 55 * of independent variables. ie. x_1 x_2 ... x_nvar. And the varaible id need to be a consequent integer start 56 * with 0. 57 * ret_vec will contains len number of doubles. Where the first nvar elements is the gradient vector, 58 * and the rest of (nvar+1)*nvar/2 elements are the upper/lower plus the diagonal part of the Hessian matrix 59 * in row format. 60 * This algorithm is inefficient, because at each nodes, it didn't check the dependency of the independent 61 * variables up to the current node. (or it is hard to do so for this setup). Therefore, it computes a full 62 * for loops over each independent variable (ie. assume they are all dependent), for those independent 63 * variables that are not dependent at the current node, zero will be produced by computation. 64 * By default the forward mode hessian routing is disabled. To enable the forward hessian interface, the 65 * compiler marco FORWARD_ENABLED need to be set equal to 1 in auto_diff_types.h 66 * 67 * + Reverse Hessian*Vector Evaluation: 68 * Simple, building a tape in the forward pass, and a reverse pass will evaluate the Hessian*vector. The implemenation 69 * also discovery the repeated subexpression and use one piece of memory on the tape for the same subexpression. This 70 * allow efficient evaluation, because the repeated subexpression only evaluate once in the forward and reverse pass. 71 * This algorithm can be called n times to compute a full Hessian, where n equals the number of independent 72 * variables. 73 * */ 74 75 typedef boost::numeric::ublas::compressed_matrix<double,boost::numeric::ublas::column_major,0,std::vector<std::size_t>,std::vector<double> > col_compress_matrix; 76 typedef boost::numeric::ublas::matrix_row<col_compress_matrix > col_compress_matrix_row; 77 typedef boost::numeric::ublas::matrix_column<col_compress_matrix > col_compress_matrix_col; 78 79 namespace AutoDiff{ 80 81 //node creation methods 82 extern PNode* create_param_node(double value); 83 extern VNode* create_var_node(double v=NaN_Double); 84 extern OPNode* create_uary_op_node(OPCODE code, Node* left); 85 extern OPNode* create_binary_op_node(OPCODE code, Node* left,Node* right); 86 87 //single constraint version 88 extern double eval_function(Node* root); 89 extern unsigned int nzGrad(Node* root); 90 extern double grad_reverse(Node* root,vector<Node*>& nodes, vector<double>& grad); 91 extern unsigned int nzHess(EdgeSet&); 92 extern double hess_reverse(Node* root, vector<Node*>& nodes, vector<double>& dhess); 93 94 //multiple constraints version 95 extern unsigned int nzGrad(Node* root, boost::unordered_set<Node*>& vnodes); 96 extern double grad_reverse(Node* root, vector<Node*>& nodes, col_compress_matrix_row& rgrad); 97 extern unsigned int nzHess(EdgeSet&,boost::unordered_set<Node*>& set1, boost::unordered_set<Node*>& set2); 98 extern double hess_reverse(Node* root, vector<Node*>& nodes, col_compress_matrix_col& chess); 99 100 #if FORWARD_ENDABLED 101 //forward methods 102 extern void hess_forward(Node* root, unsigned int nvar, double** hess_mat); 103 #endif 104 105 //utiliy methods 106 extern void nonlinearEdges(Node* root, EdgeSet& edges); 107 extern unsigned int numTotalNodes(Node*); 108 extern string tree_expr(Node* root); 109 extern void print_tree(Node* root); 110 extern void autodiff_setup(); 111 extern void autodiff_cleanup(); 112 }; 113 114 #endif /* AUTODIFF_H_ */ 115