|
| 1 | +--- |
| 2 | +title: "How to Estimate Floating Point Errors Using Automatic Differentiation" |
| 3 | +layout: post |
| 4 | +excerpt: "This tutorial demostrates how automatic differentiation can be used |
| 5 | + to estimate floating point errors of algorithms." |
| 6 | +sitemap: false |
| 7 | +permalink: /tutorials/fp_error_estimation_clad_tutorial/ |
| 8 | +date: 24-08-2021 |
| 9 | +author: Garima Singh |
| 10 | +--- |
| 11 | + |
| 12 | +*Tutorial level: Intro* |
| 13 | + |
| 14 | +For the purpose of this tutorial, we will look at a very simple estimation |
| 15 | +model that is a slight modification of clad’s in-built Taylor approximation |
| 16 | +model. The first step of translating your model to code is to understand how |
| 17 | +you could represent it mathematically as a formula. For example, the |
| 18 | +mathematical formula behind the in-built Taylor approximation model is: |
| 19 | + |
| 20 | +$$ |
| 21 | +\begin{align*} |
| 22 | +A_f = \sum\limits_{i=1}^n \frac{\partial f}{\partial x_i} \ast |x_i| \ast |\varepsilon_M| |
| 23 | +\end{align*} |
| 24 | +$$ |
| 25 | + |
| 26 | +Where, |
| 27 | +- **Af** is the absolute error in the function that you want to analyse. |
| 28 | + We will refer to this as the “target function”. |
| 29 | +- **df/dx<sub>i</sub>** is the partial derivative of the target function |
| 30 | +with respect to the i<sup>th</sup> variable. This is where clad comes in: |
| 31 | +it provides us with the derivatives we need in our models. |
| 32 | +- **x<sub>i</sub>** is the i<sup>th</sup> variable of interest, here all |
| 33 | +floating point variables in the function are our variables of interest. |
| 34 | +All values, including temporary ones with limited scope, have to be |
| 35 | +analysed to get the most accurate estimation. |
| 36 | +- **E<sub>m</sub>** is the machine epsilon, which is a system dependent value |
| 37 | + giving the maximum error in the representation of a floating point number. |
| 38 | + This is the relative error in floating point values and has to be scaled |
| 39 | + up to be of any meaning. Hence, you see the x<sub>i</sub> in multiplication. |
| 40 | + |
| 41 | +Each i<sup>th</sup> term in the above summation becomes the error contribution |
| 42 | +to the total error by the i<sup>th</sup> assignment. The summation of all |
| 43 | +those values gives us the final absolute floating point error in the target |
| 44 | +function. |
| 45 | + |
| 46 | +Now that is what clad implements, and you might be able to tell that because |
| 47 | +of the use of the machine epsilon, the estimates we get from the default model |
| 48 | +are very strict (i.e., they report a tight upper bound on the error in the |
| 49 | +function). Which means while we can say that our target function will never |
| 50 | +exceed this specific amount of floating point error, it is sometimes necessary |
| 51 | +to get even tighter bounds. To achieve that, we can go one step further by |
| 52 | +incorporating the actual errors in the values themselves. This means that we |
| 53 | +can simply subtract the current value by the same value but with lower |
| 54 | +precision and use that as the error metric. So, our modified Taylor |
| 55 | +approximation model would look something like this |
| 56 | + |
| 57 | +$$ |
| 58 | +\begin{align*} |
| 59 | +A_f = \sum\limits_{i=1}^n \frac{\partial f}{\partial x_i} \ast |x_i - (float)x_i| |
| 60 | +\end{align*} |
| 61 | +$$ |
| 62 | + |
| 63 | +Here we have cast x<sub>i</sub> down to `float` from `double`. This model |
| 64 | +however imposes a constraint that all our floating point variables should |
| 65 | +be of type `double` or higher. |
| 66 | + |
| 67 | +Once you have a solid formula for each value of interest, the only thing |
| 68 | +that remains is translating that into a function. We have a formula that |
| 69 | +we just derived in the previous section and now we have to convert that to |
| 70 | +a function. So before we begin writing a function, it is important to |
| 71 | +understand how clad is internally able to generate error estimation code. |
| 72 | + |
| 73 | +Clad works on the very simple concept of a *Recursive AST Visitor*. |
| 74 | +For brevity's sake, we will not be discussing ASTs in depth in this tutorial. |
| 75 | +In a gist, these visitors recursively visit each statement in a program in a |
| 76 | +depth-first manner. To illustrate the visitation, consider the following simple |
| 77 | + example: |
| 78 | + |
| 79 | +Let’s suppose we have a semantically correct statement in our program as follows: |
| 80 | + |
| 81 | +$$ |
| 82 | +\begin{align*} |
| 83 | +x = a + b |
| 84 | +\end{align*} |
| 85 | +$$ |
| 86 | + |
| 87 | +The AST visitation of the above statement will look like this: |
| 88 | + |
| 89 | +<div align=center style="max-width:1095px; margin:0 auto;"> |
| 90 | + <img src="/images/tutorials/fp_error_estimation_clad_tutorial/ast.gif" |
| 91 | + style="max-width:90%;"><br/> |
| 92 | + <p align="center"> |
| 93 | + </p> |
| 94 | +</div> |
| 95 | + |
| 96 | +Now, clad will generate error estimation specific code on the leftmost leaf |
| 97 | +node of the above tree using the `AssignError()` function, which we will be |
| 98 | +overriding later in the tutorial. This also implies that clad will only |
| 99 | +generate error estimates for assignment operations, more specifically, it will |
| 100 | +generate errors for the LHS of all assignment operations. |
| 101 | + |
| 102 | +Now that we have an idea of how the error estimation framework works, all we |
| 103 | +have to do is: |
| 104 | + |
| 105 | +1. Inherit from the `clad::EstimationModel` class |
| 106 | +2. Override `AssignError()` and `SetError()` functions |
| 107 | +3. Register our derived class with the |
| 108 | + `clad::plugin::ErrorEstimationModelRegistry` |
| 109 | +4. Compile our new library into a shared object |
| 110 | + |
| 111 | +And we will be ready to use our custom estimation model! |
| 112 | + |
| 113 | +First, let us make a header and a `.cpp` file to get started. For this tutorial, |
| 114 | + let us refer to our custom model as CustomModel. Once you have created the |
| 115 | + files, populate them with the following code: |
| 116 | + |
| 117 | +**CustomModel.h** |
| 118 | +```cpp |
| 119 | +#include "clad/Differentiator/EstimationModel.h" |
| 120 | +#include "clad/tools/ClangPlugin.h" |
| 121 | + |
| 122 | +class CustomModel : public clad::EstimationModel { |
| 123 | +public: |
| 124 | + clang::Expr* AssignError(clad::StmtDiff refExpr); |
| 125 | + |
| 126 | + clang::Expr* SetError(clad::VarDeclDiff declStmt); |
| 127 | +}; |
| 128 | +``` |
| 129 | +
|
| 130 | +**CustomModel.cpp** |
| 131 | +```cpp |
| 132 | +#include "CustomModel.h" |
| 133 | +
|
| 134 | +clang::Expr* CustomModel::AssignError(clad::StmtDiff refExpr) { |
| 135 | + // Populate here with your code |
| 136 | +} |
| 137 | +clang::Expr* CustomModel::SetError(clad::VarDeclDiff declStmt) { |
| 138 | + // Populate here with your code |
| 139 | +} |
| 140 | +
|
| 141 | +static clad::plugin::ErrorEstimationModelRegistry::Add<CustomModel> |
| 142 | +CX("customModel", "Custom model for error estimation in clad"); |
| 143 | +``` |
| 144 | + |
| 145 | +Before we start filling in the definitions, let’s look at the last statement in `CustomModel.cpp` |
| 146 | + |
| 147 | +```cpp |
| 148 | +static clad::plugin::ErrorEstimationModelRegistry::Add<CustomModel> |
| 149 | +CX("customModel", "Custom model for error estimation in clad"); |
| 150 | +``` |
| 151 | +
|
| 152 | +This statement tells clad to register our custom model. Without this |
| 153 | +statement, clad will be unable to locate our derived class and its functions |
| 154 | +and will throw a symbol lookup error during its usage. |
| 155 | +
|
| 156 | +Now, the only thing left to do is to populate the two functions with our custom |
| 157 | +logic of error estimation. |
| 158 | +
|
| 159 | +First, let’s look at `SetError`, which is used to initialize the error associated |
| 160 | +with each variable of interest. This function might be suitable if we know of an |
| 161 | +intrinsic error for all the variables and want that to be a part of the final |
| 162 | +error. For now, let us return a `nullptr` which clad will interpret as a 0 |
| 163 | +literal. Which means that clad will initialize all error values as 0. |
| 164 | +
|
| 165 | +```cpp |
| 166 | +clang::Expr* CustomModel::SetError(clad::VarDeclDiff declStmt) { |
| 167 | + return nullptr; |
| 168 | +} |
| 169 | +``` |
| 170 | + |
| 171 | +Now, let’s look at the `AssignError` function. Recall that we want to translate |
| 172 | +the following mathematical formula into code: |
| 173 | + |
| 174 | +$$ |
| 175 | +\begin{align*} |
| 176 | +\triangle x_i = \frac{\partial f}{\partial x_i} \ast |x_i - (float)x_i| |
| 177 | +\end{align*} |
| 178 | +$$ |
| 179 | + |
| 180 | +First, we can obtain the partial derivative of x<sub>i</sub> using `getExpr_dx()`. |
| 181 | +Here, `refExpr` refers to x<sub>i</sub> |
| 182 | + |
| 183 | +```cpp |
| 184 | +clang::Expr* CustomModel::AssignError(clad::StmtDiff refExpr) { |
| 185 | + // Get the df/dx term |
| 186 | + auto derivative = refExpr.getExpr_dx(); |
| 187 | +} |
| 188 | +``` |
| 189 | +
|
| 190 | +Now, we need to create the float subtraction expression. So, we essentially |
| 191 | +need to get x<sub>i</sub> and cast it to float and then build a subtraction |
| 192 | +expression. First let's look at how we can build a cast expression. Here we |
| 193 | +use `BuildCStyleCastExpr` as follows: |
| 194 | +
|
| 195 | +```cpp |
| 196 | +clang::Expr* CustomModel::AssignError(clad::StmtDiff refExpr) { |
| 197 | + // Get the df/dx term |
| 198 | + auto derivative = refExpr.getExpr_dx(); |
| 199 | + auto expr = refExpr.getExpr(); |
| 200 | + // Get the cast expression |
| 201 | + auto castExpr = m_Sema |
| 202 | + .BuildCStyleCastExpr( |
| 203 | + expr->getBeginLoc(), |
| 204 | + m_Context.getTrivialTypeSourceInfo(m_Context.FloatTy), |
| 205 | + expr->getEndLoc(), expr) |
| 206 | + .get(); |
| 207 | +} |
| 208 | +``` |
| 209 | + |
| 210 | +Now, let’s use the `BuildOp` shorthand to build the subtraction expression. |
| 211 | +Note, the first argument to `BuildOp` is the type of operation and the next |
| 212 | +are the LHS and RHS of that operation. |
| 213 | + |
| 214 | +```cpp |
| 215 | +clang::Expr *CustomModel::AssignError(clad::StmtDiff refExpr) { |
| 216 | + // Get the df/dx term |
| 217 | + auto derivative = refExpr.getExpr_dx(); |
| 218 | + auto expr = refExpr.getExpr(); |
| 219 | + // Get the cast expression |
| 220 | + auto castExpr = m_Sema |
| 221 | + .BuildCStyleCastExpr( |
| 222 | + expr->getBeginLoc(), |
| 223 | + m_Context.getTrivialTypeSourceInfo(m_Context.FloatTy), |
| 224 | + expr->getEndLoc(), expr) |
| 225 | + .get(); |
| 226 | + // Build subtraction operator |
| 227 | + auto subExpr = BuildOp(BO_Sub, expr, castExpr); |
| 228 | +} |
| 229 | +``` |
| 230 | +
|
| 231 | +Since clad internally builds the assignment, we only need to return the |
| 232 | +multiplication expression. We can use the same `buildOp` to get that done! |
| 233 | +
|
| 234 | +```cpp |
| 235 | +clang::Expr* CustomModel::AssignError(clad::StmtDiff refExpr) { |
| 236 | + // Get the df/dx term |
| 237 | + auto derivative = refExpr.getExpr_dx(); |
| 238 | + auto expr = refExpr.getExpr(); |
| 239 | + // Get the cast expression |
| 240 | + auto castExpr = m_Sema |
| 241 | + .BuildCStyleCastExpr( |
| 242 | + expr->getBeginLoc(), |
| 243 | + m_Context.getTrivialTypeSourceInfo(m_Context.FloatTy), |
| 244 | + expr->getEndLoc(), expr) |
| 245 | + .get(); |
| 246 | + // Build subtraction operator |
| 247 | + auto subExpr = BuildOp(BO_Sub, expr, castExpr); |
| 248 | + // Return the final multiplication operation |
| 249 | + return BuildOp(BO_Mul, derivative, subExpr); |
| 250 | +
|
| 251 | +} |
| 252 | +``` |
| 253 | + |
| 254 | +And voila! We are all done with writing our custom estimation model! |
| 255 | + |
| 256 | +Lastly, we have to compile our source into a shared object. For this we shall |
| 257 | +first set up some environment variables to simplify following the rest of the |
| 258 | +tutorial. In a terminal, run the following: |
| 259 | + |
| 260 | +```bash |
| 261 | +export CLAD_INST=$PWD/../inst; |
| 262 | +export CLAD_BASE=$PWD/../clad; |
| 263 | +``` |
| 264 | + |
| 265 | +Then execute the following: |
| 266 | + |
| 267 | +```bash |
| 268 | +CLAD_INST/bin/clang -ICLAD_BASE/include -fPIC -shared -fno-rtti -Wl,-undefined -Wl,suppress ./CustomModel.cpp -o libCustomModel.so |
| 269 | +``` |
| 270 | + |
| 271 | +The above command creates a shared library called libCustomModel.so. Now to use |
| 272 | +this shared library with clad, simply pass it via CLI when invoking clad as |
| 273 | +follows: |
| 274 | + |
| 275 | +```bash |
| 276 | +-Xclang -plugin-arg-clad -Xclang -fcustom-estimation-model -Xclang -plugin-arg-clad -Xclang ./libCustomModel.so |
| 277 | +``` |
| 278 | + |
| 279 | +That is, a typical invocation to clad will now look like this |
| 280 | + |
| 281 | +```bash |
| 282 | +CLAD_INST/bin/clang -Xclang -add-plugin -Xclang clad -Xclang -load -Xclang CLAD_INST/clad.so -ICLAD_BASE/include -x c++ -lstdc++ -Xclang -plugin-arg-clad -Xclang -fcustom-estimation-model -Xclang -plugin-arg-clad -Xclang ./libCustomModel.so some.cpp |
| 283 | +``` |
| 284 | + |
| 285 | +Yay! You have successfully generated code instrumented with your custom |
| 286 | +estimation logic! |
| 287 | + |
| 288 | +To know about what all the error estimation framework is capable of, |
| 289 | +checkout the docs [here](https://clad.readthedocs.io/en/latest/index.html). |
0 commit comments