Giter VIP home page Giter VIP logo

fastad's Introduction

FastAD

Build Status CircleCI Coverage Status Codacy Badge CII Best Practices GitHub tag (latest by date) GitHub release (latest by date) GitHub issue (latest by date) License

Table of contents

Overview

FastAD is a header-only C++ template library for automatic differentiation supporting both forward and reverse mode. It utilizes the latest features in C++17 and expression templates for efficient computation. FastAD is unique for the following:

Intuitive syntax

Syntax choice is very important for C++ developers. Our philosophy is that syntax should be as similar as possible to mathematical notation. This makes FastAD easy to use and allow users to write readable, intuitive, and simple code. See User Guide for more details.

Robustness

FastAD has been heavily unit-tested with high test coverage followed by a few integration tests. A variety of functions have been tested against analytical solutions; at machine-level precision, the derivatives coincide.

Memory Efficiency

FastAD is written to be incredibly efficient with memory usage and cache hits. The main overhead of most AD libraries is the tape, which stores adjoints. Using expression template techniques, and smarter memory management, we can significantly reduce this overhead.

Speed

Speed is the utmost critical aspect of any AD library. FastAD has been proven to be extremely fast, which inspired the name of this library. Benchmark shows over orders of magnitude improvement from existing libraries such as Adept, Stan Math Library, ADOL-C, CppAD, and Sacado (see our separate benchmark repositor ADBenchmark). Moreover, it also shows 10x improvement from the naive (and often inaccurate) finite-difference method.

Installation

First, clone the repo:

git clone https://github.com/JamesYang007/FastAD.git ~/FastAD

From here on, we will refer to the cloned directory as workspace_dir (in the example above, workspace_dir is ~/FastAD).

The library has the following dependencies:

  • Eigen3.3
  • GoogleTest (dev only)
  • Google Benchmark (dev only)

General Users

If the user already has Eigen3.3 installed in their system, they can omit the following step. For general users, if they wish to install Eigen locally, they can run

./setup.sh

from workspace_dir. This will install Eigen3.3 into workspace_dir/libs/eigen-3.3.7/build.

For those who want to install FastAD globally into the system, simply run:

./install.sh

This will build and install the header files into the system.

For users who want to install FastAD locally, run the following from workspace_dir:

mkdir -p build && cd build
cmake -DCMAKE_INSTALL_PREFIX=. ..
make install

One can set the CMAKE_INSTALL_PREFIX to anything. This example will install the library in workspace_dir/build.

For users that want to integrate FastAD in their own CMakeLists, use FetchContent (CMake >= 3.11)

include(FetchContent)
FetchContent_Declare(
        FastAD
        GIT_REPOSITORY https://github.com/JamesYang007/FastAD
        GIT_TAG v3.2.1
        GIT_SHALLOW TRUE
        GIT_PROGRESS TRUE)
FetchContent_MakeAvailable(FastAD)
# Further link target 'FastAD'

Developers

Run the following to install all of the dependencies locally:

./setup.sh dev

To build the library, run the following:

./clean-build.sh <debug/release> [other CMake flags...]

Here are the following options one can specify as a CMake flag -D...=ON (replace ... with any of the following):

  • FASTAD_ENABLE_TEST (builds tests)
  • FASTAD_ENABLE_BENCHMARK (builds benchmarks)
  • FASTAD_ENABLE_EXAMPLE (builds examples)

By default, the flags are OFF. Note that this only builds and does not install the library.

To run tests, execute the following:

cd build/<debug/release>
ctest -j6

To run benchmarks, change directory to build/<debug/release>/benchmark and run any one of the executables.

Integration

CMake

If your project is built using CMake, add the following to CMakeLists.txt in the root directory:

find_package(FastAD CONFIG REQUIRED)

If you installed the library locally, say path_to_install, then add the following:

find_package(FastAD CONFIG REQUIRED HINTS path_to_install/share)

For any program that requires FastAD, use target_link_libraries to link with FastAD::FastAD.

An example project that uses FastAD as a dependency may have a CMakeLists.txt that looks like this:

project("MyProject")
find_package(FastAD CONFIG REQUIRED)
add_executable(main src/main.cpp)
target_link_libraries(main FastAD::FastAD)

Others

Simply add the following flag when compiling your program:

-Ipath_to_install/include

An example build command would be:

g++ main.cpp -Ipath_to_install/include

If Eigen3.3 was installed locally, you must provide its path as well.

User Guide

The only header the user needs to include is fastad.

Forward Mode

Forward mode is extremely simple to use.

The only class a user will need to deal with is ForwardVar<T>, where T is the underlying data type (usually double). The API only exposes getters and setters:

ForwardVar<double> v;       // initialize value and adjoint to 0
v.set_value(1.);            // value is now 1.
double r = v.get_value();   // r is now 1.
v.set_adjoint(1.);          // adjoint is now 1.
double s = v.get_adjoint(); // s is now 1.

The rest of the work has already been done by the library with operator overloading.

Here is an example program that differentiates a complicated function:

#include <fastad>
#include <iostream>

int main()
{
    using namespace ad;

    ForwardVar<double> w1(0.), w2(1.);
    w1.set_adjoint(1.); // differentiate w.r.t. w1
    ForwardVar<double> w3 = w1 * sin(w2);
    ForwardVar<double> w4 = w3 + w1 * w2;
    ForwardVar<double> w5 = exp(w4 * w3);

    std::cout << "f(x, y) = exp((x * sin(y) + x * y) * x * sin(y))\n"
              << "df/dx = " << w5.get_adjoint() << std::endl;
    return 0;
}

We initialize w1 and w2 with values 0. and 1., respectively. We set the adjoint of w1 to 1. to indicate that we are differentiating w.r.t. w1. By default, all adjoints of ForwardVar are set to 0.. This indicates that we will be differentiating in the direction of (1.,0.), i.e. partial derivative w.r.t. w1. Note that user could also set the adjoint for w2, but this will compute the directional derivative multiplied by the norm of (w1, w2). After computing the desired expression, we get the directional derivative by calling get_adjoint() on the final ForwardVar object.

Reverse Mode

Basic Usage

The most basic usage simply requires users to create Var<T, ShapeType> objects. T denotes the underlying value type (usually double). ShapeType denotes the general shape of the variable. It must be one of ad::scl, ad::vec, ad::mat corresponding to scalar, (column) vector, and matrix, respectively.

Var<double, scl> x;
Var<double, vec> v(5);    // set size to 5
Var<double, mat> m(2, 3); // set shape to 2x3

From here, one can create complicated expressions by invoking a wide range of functions (see Quick Reference for a full list of expression builders).

As an example, here is an expression to differentiate sin(x) + cos(v):

auto expr = (sin(x) + cos(v));

Note that this represents a vector expression, since sin(x) is a scalar expression but cos(v) is a vectorized function on a vector, which is again a vector expression.

Before we differentiate, the expression is required to "bind" to a storage for the values and adjoints of intermediate expression nodes. The reason for this design is for speed purposes and cache hits. If the user wishes to manage this storage, they can do this:

auto size_pack = expr.bind_cache_size();
std::vector<double> val_buf(size_pack(0));
std::vector<double> adj_buf(size_pack(1));
expr.bind_cache({val_buf.data(), adj_buf.data()});

The bind_cache_size() will return exactly how many doubles are needed for values and adjoints, respectively, of type util::SizePack, which is an alias for Eigen::Array<size_t, 2, 1>. and bind(util::PtrPack<double>) will bind itself to that region of memory. It is encouraged to create the pointer pack object using initializer list as shown above. This pattern occurs so often that if the user does not care about managing this, they should use the following helper function:

auto expr_bound = ad::bind(sin(x) + cos(v));

ad::bind will return a wrapper class that wraps the expression and at construction binds it to a privately owned storage in the same way described above.

If the expression is not bound to any storage, it will lead to segfault!

To differentiate the expression, simply call the following:

auto f = ad::autodiff(expr_bound, seed);

// or if the raw expression is manually bound,
auto f = ad::autodiff(expr, seed);

where seed is the initial adjoint for the root of the expression. If the expression is scalar, seed is a literal (double) and the default value is 1, so the user does not have to input anything. If the expression is multi-dimensional, seed does not have a default value, must be of type Eigen::Array, and must have the same dimensions as the expression. autodiff will return the evaluated function value. This return value is T if it is a scalar expression, and otherwise, Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, ...>> where ... depends on the shape of the expression (1 if vector, Eigen::Dynamic if matrix).

You can retrieve the adjoints by calling get_adj(i,j) or get_adj() (with no arguments) from x, v like so:

x.get_adj(0,0); // (1) get adjoint for x 
v.get_adj(2,0); // (2) get adjoint for v at index 2
x.get_adj();    // (3) get full adjoint (same as (1))
v.get_adj();    // (4) get full adjoint

The full code for this example is the following:

#include <fastad>
#include <iostream>

int main()
{
    using namespace ad;

    Var<double, scl> x(2);
    Var<double, vec> v(5);
    
    // randomly generate values for v
    v.get().setRandom();

    // create AD expression bound to storage
    auto expr_bound = bind(sin(x) + cos(v));
    
    // seed to get gradient of function at index 2
    Eigen::Array<double, Eigen::Dynamic, 1> seed(v.size());
    seed.setZero();
    seed[2] = 1;

    // differentiate
    auto f = autodiff(expr_bound, seed);

    std::cout << x.get_adj() << std::endl;
    std::cout << v.get_adj(2,0) << std::endl;

    return 0;
}

Note: once you have differentiated an expression, you must reset the adjoints of all variables to 0 before differentiating again. This includes placeholder variables (see below). To that end, we provide a member function for Var called reset_adj().

Here is a more complicated example:

#include <fastad>

int main()
{
    Var<double, vec> v1(6);
    Var<double, vec> v2(5);
    Var<double, mat> M(5, 6);
    Var<double, vec> w(5);
    Var<double, scl> r;

    auto& v1_raw = v1.get();  // Eigen::Map
    auto& v2_raw = v2.get();  // Eigen::Map
    auto& M_raw = M.get();  // Eigen::Map

    // initialize...

    auto expr = bind((
        w = ad::dot(M, v1) + v2,
        r = sum(w) * sum(w * v2)
    ));

    autodiff(expr);

    std::cout << v1.get_adj(0,0) << std::endl;  // adjoint of v1 at index 0
    std::cout << v2.get_adj(1,0) << std::endl;  // adjoint of v2 at index 1
    std::cout << M.get_adj(1,2) << std::endl;   // adjoint of M at index (1,2)

    return 0;
}

Placeholder

In the previous example, we used a placeholder expression, which is of the form v = expr. You can use placeholders to greatly speed up the performance and also save a lot of memory.

Consider the following expression:

auto expr = (sin(x) + cos(v) + sum(cos(v)));

When there are common expressions (like cos(v)), they will be evaluated multiple times unnecessarily.

Placeholder expressions are created by using operator= with a Var and an expression:

Var<double, scl> x;
Var<double, vec> v(5);
Var<double, vec> w(v.size());
auto expr = (
    w = cos(v),
    sin(x) + w + sum(w)
);

This will only evaluate cos(v) once, and reuse the results for the subsequent expressions by using w.

While this is not specific to placeholder expressions, operator, is usually invoked to "glue" many placeholder expressions. However, one can certainly glue any kinds of expressions, if they wish.

Advanced Usage

For advanced users who need to get more low-level control over the memory for values and adjoints for all variables, they can use VarView<T, ShapeType>. All of the discussion above holds for VarView objects. In fact, when we build an expression out of Var of VarView, we convert all of them to VarViews so that the expression is solely a viewer.

VarView objects do not own the values and adjoints, but views them. Here is an example program that binds the viewers to a contiguous chunk of memory:

VarView<double, scl> x;
VarView<double, vec> v(3);
VarView<double, vec> w(3);

std::vector<double> vals(x.size() + v.size());
std::vector<double> adjs(x.size() + v.size());
std::vector<double> w_vals(w.size());
std::vector<double> w_adjs(w.size());

// x binds to the first element of storages
double* val_next = x.bind(vals.data());
double* adj_next = x.bind_adj(adjs.data());

// v binds starting from 2nd element of storages
v.bind(val_next);
v.bind_adj(adj_next);

// bind placeholders to a separate storage region
w.bind(w_vals.data());
w.bind_adj(w_adjs.data());

auto expr = (
    w = cos(v),
    sin(x) + w + sum(w)
);

Applications

Black-Scholes Put-Call Option Pricing

The following is an example of computing deltas in Black-Scholes model using FastAD. This example was taken from the autodiff library in boost.

#include <fastad>
#include <iostream>

enum class option_type {
    call, put
};

// Standard Normal CDF
template <class T>
inline auto Phi(const T& x)
{
    return 0.5 * (ad::erf(x / std::sqrt(2.)) + 1.);
}

// Generates expression that computes Black-Scholes option price
template <option_type cp, class Price, class Cache>
auto black_scholes_option_price(const Price& S,
                                double K,
                                double sigma,
                                double tau,
                                double r,
                                Cache& cache)
{
    cache.resize(3);
    double PV = K * std::exp(-r * tau);
    auto common_expr = (
            cache[0] = ad::log(S / K),
            cache[1] = (cache[0] + ((r + sigma * sigma / 2.) * tau)) / 
                            (sigma * std::sqrt(tau)),
            cache[2] = cache[1] - (sigma * std::sqrt(tau))
    );
    if constexpr (cp == option_type::call) {
        return (common_expr,
                Phi(cache[1]) * S - Phi(cache[2]) * PV);
    } else {
        return (common_expr,
                Phi(-cache[2]) * PV - Phi(-cache[1]) * S);
    }
}

int main()
{
    double K = 100.0;        
    double sigma = 5;        
    double tau = 30.0 / 365; 
    double r = 1.25 / 100;   
    ad::Var<double> S(105);  
    std::vector<ad::Var<double>> cache;

    auto call_expr = ad::bind(
            black_scholes_option_price<option_type::call>(
                S, K, sigma, tau, r, cache));

    double call_price = ad::autodiff(call_expr);

    std::cout << call_price << std::endl;
    std::cout << S.get_adj() << std::endl;

    // reset adjoints before differentiating again
    S.reset_adj();
    for (auto& c : cache) c.reset_adj();

    auto put_expr = ad::bind(
            black_scholes_option_price<option_type::put>(
                S, K, sigma, tau, r, cache));

    double put_price = ad::autodiff(put_expr);

    std::cout << put_price << std::endl;
    std::cout << S.get_adj() << std::endl;

    return 0;
}

We observed the same output as the one shown in boost for the prices and deltas (S's adjoints).

Quadratic Expression Differential

In ML applications, user usually provide a full parameter vector to an optimizer and write an objective function to calculate objective value and gradient. The gradient pointer is provided by the optimizer and user fill values. Then it will be more convinent to use VarView to bind the gradient pointer as buffer.

Here is an example of differentiating a quadratic expression x^T*Sigma*x using VarView.

#include <iostream>
#include "fastad"
#include <Eigen/src/Core/Matrix.h>

int main() {
    using namespace ad;

    // Generating buffer.
    Eigen::MatrixXd x_data(2, 1);
    x_data << 0.5, 0.6;
    Eigen::MatrixXd x_adj(2, 1);
    x_adj.setZero(); // Set adjoints to zeros.

    // Initialize variable.
    VarView<double, mat> x(x_data.data(), x_adj.data(), 2, 1);

    // Initialize matrix.
    Eigen::MatrixXd _Sigma(2, 2);
    _Sigma << 2, 3, 3, 6;
    std::cout << _Sigma << std::endl;
    auto Sigma = constant(_Sigma);

    // Quadratic expression: x^T*Sigma*x
    auto expr = bind(dot(dot(transpose(x), Sigma), x));
    // Seed
    Eigen::MatrixXd seed(1, 1);
    seed.setOnes(); // Usually seed is 1. DONT'T FORGET!
    // Auto differential.
    auto f = autodiff(expr, seed.array());

    // Print results.
    std::cout << "f: " << f << std::endl;
    std::cout << x.get() << std::endl;     //[0.5, 0.6]
    std::cout << x.get_adj() << std::endl; //[5.6, 10.2]

    return 0;
}

Simple Linear Regression Model

In a regression model, one has many rows of data. A loop is needed to calculate loss of each row.

#include "fastad"
#include <iostream>

int main() {
    using namespace ad;
    // Create data matrix.
    Eigen::MatrixXd X(5, 2);
    X << 1, 10, 2, 20, 3, 30, 4, 40, 5, 50;
    Eigen::VectorXd y(5);
    y << 32, 64, 96, 128, 160; // y=2*x1+3*x2

    // Generating buffer.
    Eigen::MatrixXd theta_data(2, 1);
    theta_data << 1, 2;
    Eigen::MatrixXd theta_adj(2, 1);
    theta_adj.setZero(); // Set adjoints to zeros.

    // Initialize variable.
    VarView<double, mat> theta(theta_data.data(), theta_adj.data(), 2, 1);

    // Create expr. Use a row buffer to store data. Then we only need to manipulate data when
    // looping.
    Eigen::MatrixXd x_row_buffer = X.row(0);
    auto xi = constant_view(x_row_buffer.data(), 1, X.cols());
    Eigen::MatrixXd y_row_buffer = y.row(0);
    auto yi = constant_view(y_row_buffer.data(), 1, y.cols());
    auto expr = bind(pow<2>(yi - dot(xi, theta)));

    // Seed
    Eigen::MatrixXd seed(1, 1);
    seed.setOnes(); // Usually seed is 1. DONT'T FORGET!

    // Loop over each row to calulate loss.
    double loss = 0;
    for (int i = 0; i < X.rows(); ++i) {
        x_row_buffer = X.row(i);
        y_row_buffer = y.row(i);

        auto f = autodiff(expr, seed.array());
        loss += f.coeff(0);
    }

    // Print results.
    std::cout << "loss: " << loss << std::endl; // 6655
    std::cout << theta.get() << std::endl;      //[1, 2]
    std::cout << theta.get_adj() << std::endl;  //[-1210, -12100]

    theta_adj.setZero(); // Reset differential to zero after one full pass.

    return 0;
}

3-layer Neural Network with Jump Connection

Here is an example of 3-layer neural network. The result has been confirmed by symbolic differential using Mathematica. See test/reverse/util/GenTestData.nb for Mathematica code used.

There is also a full independent example that using ceres to train this network in examples/ceres_3layer_neural_net. ceres has it's own automatic differential tool for non-linear least square problem which is unavailiable for the more versatile GradientProbelm. That's why we need FastAD. To run this example, you need to install and Eigen and ceres.

#include "fastad"
#include <iostream>

using namespace ad;
// A3*s(A2*s(A1*x+b1)+b2+(A1*x+b1))+b3
struct NN {
    Eigen::MatrixXd X, y;
    NN(const Eigen::MatrixXd &X, const Eigen::VectorXd &y)
        : X(X), y(y){

                };
    double loss(const double *parm, double *grad) {
        Eigen::Map<Eigen::VectorXd> g(grad, 25);
        g.setZero();

		//Create variables.
        VarView<double, mat> A1(const_cast<double *>(parm), grad, 3, 2);
        VarView<double, mat> b1(const_cast<double *>(parm + 6), grad + 6, 3, 1);
        VarView<double, mat> A2(const_cast<double *>(parm + 9), grad + 9, 3, 3);
        VarView<double, mat> b2(const_cast<double *>(parm + 18), grad + 18, 3, 1);
        VarView<double, mat> A3(const_cast<double *>(parm + 21), grad + 21, 1, 3);
        VarView<double, mat> b3(const_cast<double *>(parm + 24), grad + 24, 1, 1);

        // Data buffer
        Eigen::VectorXd x_row_buffer = X.row(0);
        auto xi = constant_view(x_row_buffer.data(), X.cols(), 1);
        Eigen::VectorXd y_row_buffer = y.row(0);
        auto yi = constant_view(y_row_buffer.data(), y.cols(), 1);
        
        // Expression
        auto x1 = dot(A1, xi) + b1;
        auto y1 = sigmoid(x1);
        auto y2 = sigmoid(dot(A2, y1) + b2) + x1;
        auto y3 = dot(A3, y2) + b3;
        auto residual_norm2 = pow<2>(yi - y3);
        auto expr = bind(residual_norm2);

        // Seed
        Eigen::MatrixXd seed(1, 1);
        seed.setOnes(); // Usually seed is 1. DONT'T FORGET!

        // Loop over each row to calulate loss.
        double loss = 0;
        for (int i = 0; i < X.rows(); ++i) {
            x_row_buffer = X.row(i);
            y_row_buffer = y.row(i);

            auto f = autodiff(expr, seed.array());
            loss += f.coeff(0);
        }
        return loss;
    };
};

int main() {

    // Create data matrix.
    Eigen::MatrixXd X(5, 2);
    X << 1, 10, 2, 20, 3, 30, 4, 40, 5, 50;
    Eigen::VectorXd y(5);
    y << 32, 64, 96, 128, 160; // y=2*x1+3*x2

    // Generating parameter buffer and NN.
    Eigen::VectorXd parm_data(25);
    parm_data << 0.043984, 0.960126, -0.520941, -0.800526, -0.0287914, 0.635809, 0.584603,
        -0.443382, 0.224304, 0.97505, -0.824084, 0.2363, 0.666392, -0.498828, -0.781428, -0.911053,
        -0.230156, -0.136367, 0.263425, 0.841535, 0.920342, 0.65629, 0.848248, -0.748697, 0.21522;

    Eigen::VectorXd grad_data(25);
    grad_data.setZero(); // Set adjoints to zeros.
    NN net(X, y);

    auto loss = net.loss(parm_data.data(), grad_data.data());

    // Print results.
    std::cout << "loss: " << loss << std::endl;      // 6655
    std::cout << "parm: " << parm_data << std::endl; //[1, 2]
    std::cout << "diff: " << grad_data << std::endl; //[-1210, -12100]

    return 0;
}

Quick Reference

Forward

ForwardVar:

  • class representing a variable for forward AD
  • set_value(T x): sets value to x
  • get_value(): gets underlying value
  • set_adjoint(T x): sets adjoint to x
  • get_adjoint(): gets underlying adjoint

Unary Functions:

  • unary minus: operator-
  • trig functions: sin, cos, tan, asin, acos, atan
  • others: exp, log, sqrt

Operators:

  • binary: +,-,*,/
  • increment: +=

Reverse

Shape Types:

  • ad::scl, ad::vec, ad::mat

VarView<T, ShapeType=scl>:

  • This is only useful for users who really want to optimize for performance
  • ShapeType must be one of the types listed above
  • T is the underlying value type
  • VarView(T* v, T* a, rows=1, cols=1):
    • constructs to view values starting from v, adjoints starting from a, and has the shape of rows x cols.
    • vector shapes must pass rows
    • matrix shapes must pass both rows and cols
  • VarView()
    • constructs with nullptrs
  • .bind(T* begin): views values starting from begin
  • .bind_adj(T* begin): views adjoints starting from begin

Var<T, ShapeType=scl>:

  • A Var is a VarView (views itself)
  • Main difference with VarView is that it owns the values and adjoints
  • Users will primarily use this class to represent AD variables.
  • API is same as VarView

Unary Functions (vectorized if multi-dimensional):

  • unary minus: operator-
  • trig functions: sin, cos, tan, asin, acos, atan
  • Hyperbolic: sinh, cosh, tanh
  • Neural network activations: sigmoid
  • others: exp, log, sqrt, erf

Operators:

  • binary: +,-,*,/
  • modification: +=, -=, *=, /=
  • comparison: <,<=,>,>=,==,!=,&&,||
    • Note: && and || are undefined behavior for multi-dimensional non-boolean expressions
  • placeholder: operator=
    • only overloaded for VarView expressions
  • glue: operator,
    • any expressions can be "glued" using this operator
  • unary minus: operator-
  • trig functions: sin, cos, tan, asin, acos, atan
  • others: exp, log, sqrt

Special Expressions:

  • ad::constant(T):
  • ad::constant(const Eigen::Vector<T, Eigen::Dynamic, 1>&):
  • ad::constant(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>&):
  • ad::constant_view(T*):
  • ad::constant_view(T*, rows):
  • ad::constant_view(T*, rows, cols):
  • ad::det<policy>(m):
    • determinant of matrix m
    • policy must be one of: DetFullPivLU, DetLDLT, DetLLT
      • see Eigen documentation for each of these (without the prefix Det) for when they apply.
  • ad::dot(m, v):
    • represents matrix product with a matrix and a (column) vector
  • ad::for_each(begin, end, f):
    • generalization of operator,
    • represents evaluating expressions generated by f when fed with elements from begin to end.
  • ad::if_else(cond, if, else):
    • represents an if-else statement
    • cond MUST be a scalar expression
    • if and else must have the exact same shape
  • ad::log_det<policy>(m)
    • same as det<policy>(m) but computes log-abs-determinant
    • policy must be one of: LogDetFullPivLU, LogDetLDLT, LogDetLLT
  • ad::norm(v):
    • represents the squared norm of a vector or Frobenius norm for matrix
  • ad::pow<n>(e):
    • compile-time known, integer-powered expression
  • ad::prod(begin, end, f):
    • represents the product of expressions generated by f when fed with elements from begin to end.
  • ad::prod(e):
    • represents the product of all elements of the expression e
    • e.g. if e is a vector expression, it represents the product of all its elements.
  • ad::sum(begin, end, f):
  • ad::sum(e):
    • same as prod but represents summation
  • ad::transpose(e):
    • matrix or vector transpose.

Stats Expressions: All log-pdfs are adjusted to omit constants. Parameters can have various combinations of shapes and follow the usual vectorized notion.

  • ad::bernoulli(x, p)
  • ad::cauchy_adj_log_pdf(x, loc, scale)
  • ad::normal_adj_log_pdf(x, mu, s)
  • ad::uniform_adj_log_pdf(x, min, max)
  • ad::wishart_adj_log_pdf(X, V, n)

Contact

If you have any questions about FastAD, please open an issue. When opening an issue, please describe in the fullest detail with a minimal example to recreate the problem.

For other general questions that cannot be resolved through opening issues, feel free to send me an email.

Contributors

James Yang Kent Hall Jean-Christophe Ruel ZhouYao
JamesYang007 Kent Jean-Christophe Ruel
github.com/JamesYang007 github.com/kentjhall github.com/jeanchristopheruel github.com/kilasuelika

Third Party Tools

Many third party tools were used for this project.

License

fastad's People

Contributors

jamesyang007 avatar kentjhall avatar kilasuelika avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar

fastad's Issues

using FastAD were Scalar type is a template parameter ?

Is it possible to use FastAD reverse mode to compute the derivative of an algorithm

template <class Scalar> Scalar f( const std::vector<Scalar>& x)

where f(x) only uses the operators + , -, *, /, and = operations with Scalar objects.
If so, is there an example that demonstrates this ?

LeafNode new member function: reset_value_ptr, reset_adj_ptr

Currently, when LeafNodes get copied, their internal value and adjoint pointers get copied to point to the original value and adjoint. This is usually the desired behavior, but sometimes they may be pointing to garbage memory. We should be able to reset these pointers to point to its own value and adjoint.

Supporting mat transpose.

mat transpose is very common in scientific computing. If we want to support that, can we write a function in reverse/core/unary.hpp to directly transpose underneath values and adjoints?

I would like to make a PR but I'm not sure if the above idea is correct.

Extending black-scholes example

It was pretty straightforward to extend the example to also return a derivative with respect to the volatitliy parameter ("vega") (see here for the Rcpp-wrapped example).

But when I (somewhat mechanically) try to roll it to the time (and rate) derivatives (after taking care of the actual expression, i.e. using ad::sqrt() and ad::exp()) I am hitting a wall in the brace initialiser. Even in the simpler case of just adding "tau":

black_scholes.cpp: In function ‘Rcpp::NumericMatrix black_scholes(double, double, double, double, double)’:
black_scholes.cpp:66:25: error: cannot convert ‘<brace-enclosed initializer list>’ to ‘ad::core::ValueAdjView<double, ad::scl>::ptr_pack_t’ {aka ‘ad::util::PtrPack<double>’}
   66 |     call_expr.bind_cache({del_buf_c.data(), veg_buf_c.data(), tau_buf_c.data()});
      |     ~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

biting with

   94 |     ptr_pack_t bind_cache(ptr_pack_t begin)
      |                           ~~~~~~~~~~~^~~~~

I must be doing something obvious wrong. I glanced at the documentation and saw the richer vector and matrix valued expression but here I think I just want a larger list of scalars. Is that possible?

Remove armadillo dependency

Remove armadillo dependency with self-implementation of matrix.

CMake, include/fastad_bits/ files, and test/ files, should not reference armadillo at all or USE_ARMA preprocessor flag.

Jacobian

Can we get any example when we have an array of functions, rather than a single one. This is quite a common case, when we need to compute Jacobian. I think I have some idea how to do that, but I am not sure this would be efficient.

Some general qs. about FastAD (new to C++)

I sent these questions in an email to Dr Yang but I thought I would post here too

I have read the paper on FastAD and I am very interested in using fastAD for my algorithm

I was wondering does FastAD work with Rcpp? If so how can I install it? I think it should be possible but I just wanted to check (i'm new to C++)

I have used "autodiff" library (https://autodiff.github.io/) however I have found it to not be much faster than numerical differentiation for my application - have you used this before? I noticed in the paper you didn't benchmark against it

Also I was wondering if it possible to compute a gradient w.r.t a std::vector filled with eigen matrices? (or any other 3D or higher dimensional structure)? or will all the parameters need to be input into a vector or matrix and then reformatted back into the container needed for the rest of the model afterwards?

Is it possible to do use fastAD just within a standard function (rather than using "int main()" etc)? Im new to C++ and have just been using functions for everything (also using it through R via Rcpp)

adding convenience operators

Hey,
I am trying out your lib and it seems nice. I think it would be nice to have some convenience operators like += etc. Also I was wondering why it is not possible to use both normal floating point values in arithmetic expressions involving ad::Var? I think it would be nice to be able to do something like
Var<double> a, b{2}; a = 5. * b * b;

Var<T> move constructor implementation

#18 this issue is related to current issue.

Should Var move its contents by copying base DualNum and setting the pointers current base members? Otherwise, default move ctor will make them point to de-allocated memory.

Reduce standard to C++14

The current library only supports from C++17. It would be nice to have it support from C++14.

Enable CMake integration with FetchContent

It would be very convenient to enable CMake integration with FetchContent like so

include(FetchContent)
FetchContent_Declare(
        FastAD
        GIT_REPOSITORY https://github.com/JamesYang007/FastAD
        GIT_TAG v3.2.1
        GIT_SHALLOW TRUE
        GIT_PROGRESS TRUE)
FetchContent_MakeAvailable(FastAD)

However, it would require to put FASTAD_ENABLE_TEST OFF by default in the CMakeLists.

Delta Function

In practice, it may be useful to have delta-functions.

Implement delta-function.

Problem using example in readme

I am getting the following messages when I try to compile one of the FastAD examples:

temp.cpp: In function ‘int main()’:
temp.cpp:12:22: error: no matching function for call to ‘autodiff(ad::core::ExprBind<ad::core::BinaryNode<ad::core::Add, ad::core::UnaryNode<ad::core::Sin, ad::VarView<double, ad::scl> >, ad::core::UnaryNode<ad::core::Cos, ad::VarView<double, ad::vec> > > >&)’
   12 |     auto f = autodiff(expr_bound);
      |              ~~~~~~~~^~~~~~~~~~~~
In file included from include/fastad_bits/reverse/core.hpp:7,
                 from include/fastad_bits/reverse.hpp:3,
                 from include/fastad:4,
                 from temp.cpp:1:

Below are the steps to reproduce this:

Step 1:

git clone https://github.com/JamesYang007/FastAD.git fastad.git

Step 2:

cd fastad.git ; ./setup.sh

Step 3:
Create the file temp.cpp with the following contents (from one of the FastAD readme examples):

#include <fastad>
#include <iostream>

int main()
{
    using namespace ad;

    Var<double, scl> x(2);
    Var<double, vec> v(5);

    auto expr_bound = bind(sin(x) + cos(v));
    auto f = autodiff(expr_bound);

    std::cout << x.get_adj() << std::endl;
    std::cout << v.get_adj(2,0) << std::endl;

    return 0;
}

Step 4:

g++ temp.cpp -o temp -I include -I libs/eigen-3.3.7/build/include/eigen3

New Features

  • Discrete distributions:
    • Binomial
    • Poisson
  • Continuous distributions:
    • Beta
    • Chi-squared
    • Exponential
    • F
    • Gamma
    • t
  • Simplex distributions:
    • Dirichlet
  • Covariance distributions:
    • Inverse wishart

Non-templated expression container?

Is there a standard way to store a FastAD expression in a "generic" expression object?

The BaseExpr class is templated on the derived type, which means it's not suitable for use in a class that can store any FastAD expression.

Maybe a variant of the std::any concept?

Remove CRTP and Replace with Concepts

There are many uses of CRTP such as ADExpression that can be simplified and made more robust using self-implemented Concepts.

Types to consider:

  • ADExpression
  • LeafNode
  • UnaryNode
  • BinaryNode
  • EqNode
  • GlueNode
  • SumNode
  • ForEach
  • ProdNode

There are meta-programming tools such as is_glue_eq, glue_size that should be changed/moved.

How to get nth derivative?

Vec<double> a, b {1,2,3};
auto expr = (a = b*b*b+b);

How do I get 3rd derivative of a wrt b? (d3a/db3)

Improvement to Forward mode

Hi @JamesYang007!

This repo seems pretty nice, especially your philosophy as laid out in your paper of arguing the benefit of having a pair of pointers of matrices, rather than a matrix of dual numbers.

I had been using https://github.com/autodiff/autodiff for a while which overloads Eigen's scalar type (i.e. the latter approach) to use a matrix of dual numbers, and I think there are quite a bit of overhead (and cache misses) when compared to the reference function (matrix of doubles). I wanted to test out your repo, but realised that this repo had been mainly focusing on reversed mode rather than forward mode (which is the focus of https://github.com/autodiff/autodiff). Do you have any plans to make some of the implementations in https://github.com/JamesYang007/FastAD/tree/master/include/fastad_bits/reverse/core applicable to both modes? More than that, it seems like forward mode right now only works with scalar types (float/double) rather than matrix/vector type?

Finally, one huge benefit of https://github.com/autodiff/autodiff is that it is immediately applicable to a lot of existing function/solvers (since it is using custom scalar type to do autodiff), while your approach requires manual work to implements different operators (e.g. all the ad::XXX custom operator rather than the Eigen::XXX operator) and https://github.com/autodiff/autodiff immediately works with any function that takes a templated Eigen argument. Do you have any thoughts on that? (I had thought of a possible approach could be that we can extend MatrixBase rather than using custom scalar type to keep track of the adjoints of the input variables.)

Conditional Statements

Some functions may not be continuous or differentiable at a point, but can be made so by analytic continuation. For example, sin(x)/x is undefined at x = 0, but can be made continuous by defining the function at x = 0 to be 1.

Implement conditional-statements such that an expression may be evaluated differently depending on input x.

Checkpoint

Create a new ad expression responsible for representing already-computed ad expression.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.