From 2b788e78fb2ea5492a1f7a9f6e66349655486669 Mon Sep 17 00:00:00 2001 From: CarloLucibello Date: Tue, 20 Feb 2018 23:19:48 -0500 Subject: [PATCH 1/3] add jacobian, hessian, jvp, vjp --- src/AutoGrad.jl | 2 ++ src/highorder.jl | 61 +++++++++++++++++++++++++++++++++++++++++++++++ test/highorder.jl | 10 ++++++++ 3 files changed, 73 insertions(+) create mode 100644 src/highorder.jl diff --git a/src/AutoGrad.jl b/src/AutoGrad.jl index f27c63d..1017703 100644 --- a/src/AutoGrad.jl +++ b/src/AutoGrad.jl @@ -21,6 +21,7 @@ end importall Base # defining getindex, sin, etc. export grad, gradloss, check_grads, gradcheck, gradcheckN, getval export @primitive, @zerograd, recorder, Rec, Grad # the last three are required for the macros to work +export jacobian, hessian, jvp, vjp datapath = joinpath(dirname(@__FILE__),"..","data") include("core.jl") @@ -29,6 +30,7 @@ include("dotfunc.jl") include("gradcheck.jl") include("util.jl") include("interfaces.jl") +include("highorder.jl") include("base/reduce.jl") include("base/number.jl") include("base/float.jl") diff --git a/src/highorder.jl b/src/highorder.jl new file mode 100644 index 0000000..ab26bfc --- /dev/null +++ b/src/highorder.jl @@ -0,0 +1,61 @@ +""" + jacobian(f) + +Computes the Jacobian of the vector valued function `f`. +For `n = length(x)`, `m = length(f(x))`, the Jacobian +matrix `jacobian(f)(x)` has dimension `m x n`. +""" +function jacobian(f) + x -> begin + m = length(f(x)) + hcat([grad(x->f(x)[i])(x) for i=1:m]...)' + end +end + +""" + hessian(f) + +Compute the hessian of the scalar valued function `f`. +It is equivalent to `jacobian(grad(f))`. +""" +hessian(f) = jacobian(grad(f)) + + +""" + vjp(f) + +Returns the vector-Jacobian product operator +for the vector valued function `f`: + + vjp(f)(x, v) == v' * jacobian(f)(x) + +This is more efficient than computing the full Jacobian +if only a few products are needed. See also `jvp`. +""" +function vjp(f) + (x, v) -> begin + grad(x -> v'f(x))(x)' + end +end + + + +# ref https://j-towns.github.io/2017/06/12/A-new-trick.html +""" + jvp(f) + +Returns the Jacobian-vector product operator +for the vector valued function `f`: + + jvp(f)(x, u) == jacobian(f)(x) * u + +This is more efficient than computing the full Jacobian +if only a few products are needed. See also `vjp`. +""" +function jvp(f) + (x, u) -> begin + v = f(x) # any value of v should do + g = v -> vjp(f)(x, v)' + vjp(g)(v, u)' + end +end \ No newline at end of file diff --git a/test/highorder.jl b/test/highorder.jl index 71c52e9..32079e3 100644 --- a/test/highorder.jl +++ b/test/highorder.jl @@ -10,6 +10,16 @@ include("header.jl") g7 = grad(g6); @test g7(1)==-cos(1) g8 = grad(g7); @test g8(1)==sin(1) g9 = grad(g8); @test g9(1)==cos(1) + + A = rand(2, 3) + f(x) = A*x + jacobian(f)(rand(3)) == A + + v = rand(3) + jvp(f)(rand(3), v) == A*v + + u = rand(2) + vjp(f)(rand(3), u) == u'A end nothing From 325b44e4fb148d1d18bd4a9976654c969e073566 Mon Sep 17 00:00:00 2001 From: CarloLucibello Date: Tue, 20 Feb 2018 23:29:51 -0500 Subject: [PATCH 2/3] add hvp, vhp --- src/AutoGrad.jl | 2 +- src/highorder.jl | 39 +++++++++++++++++++++++++++++---------- test/highorder.jl | 10 ++++++++-- 3 files changed, 38 insertions(+), 13 deletions(-) diff --git a/src/AutoGrad.jl b/src/AutoGrad.jl index 1017703..14b4f1f 100644 --- a/src/AutoGrad.jl +++ b/src/AutoGrad.jl @@ -21,7 +21,7 @@ end importall Base # defining getindex, sin, etc. export grad, gradloss, check_grads, gradcheck, gradcheckN, getval export @primitive, @zerograd, recorder, Rec, Grad # the last three are required for the macros to work -export jacobian, hessian, jvp, vjp +export jacobian, hessian, jvp, vjp, hvp, vhp datapath = joinpath(dirname(@__FILE__),"..","data") include("core.jl") diff --git a/src/highorder.jl b/src/highorder.jl index ab26bfc..332b7bc 100644 --- a/src/highorder.jl +++ b/src/highorder.jl @@ -12,15 +12,6 @@ function jacobian(f) end end -""" - hessian(f) - -Compute the hessian of the scalar valued function `f`. -It is equivalent to `jacobian(grad(f))`. -""" -hessian(f) = jacobian(grad(f)) - - """ vjp(f) @@ -58,4 +49,32 @@ function jvp(f) g = v -> vjp(f)(x, v)' vjp(g)(v, u)' end -end \ No newline at end of file +end + +""" + hessian(f) + +Compute the hessian of the scalar valued function `f`. +It is equivalent to `jacobian(grad(f))`. +""" +hessian(f) = jacobian(grad(f)) + +hvp(f) = jvp(grad(f)) + +""" + vhp(f) + +Returns the vector-Hessian product operat. +It is equivalent to `vjp(grad(f))` +""" +vhp(f) = vjp(grad(f)) + + +""" + hvp(f) + +Returns the Hessian product operat. +It is equivalent to `jvp(grad(f))` +""" +hvp(f) = vjp(grad(f))' # can use vjp instead of jvp since the + # hessian is symmetric \ No newline at end of file diff --git a/test/highorder.jl b/test/highorder.jl index 32079e3..a0d7566 100644 --- a/test/highorder.jl +++ b/test/highorder.jl @@ -11,13 +11,19 @@ include("header.jl") g8 = grad(g7); @test g8(1)==sin(1) g9 = grad(g8); @test g9(1)==cos(1) + + A = Symmetric(rand(3, 3)) + f(x) = x'*A*x/2 + hessian(f)(rand(3)) == A + hvp(f)(rand(3), v) == A*v + u = rand(2) + vhp(f)(rand(3), u) == u'A + A = rand(2, 3) f(x) = A*x jacobian(f)(rand(3)) == A - v = rand(3) jvp(f)(rand(3), v) == A*v - u = rand(2) vjp(f)(rand(3), u) == u'A end From 395b095011f48d79fca50f206aaf649acfa8cf7c Mon Sep 17 00:00:00 2001 From: CarloLucibello Date: Wed, 21 Feb 2018 01:19:52 -0500 Subject: [PATCH 3/3] fixes and readme --- README.md | 54 +++++++++++++++++++++++++++++++++++++++++++---- src/highorder.jl | 13 +++++------- test/highorder.jl | 29 +++++++++++++++++-------- 3 files changed, 75 insertions(+), 21 deletions(-) diff --git a/README.md b/README.md index 7e7c124..efa5b94 100644 --- a/README.md +++ b/README.md @@ -28,12 +28,12 @@ for a description of how the code works in detail. ## Installation You can install AutoGrad in Julia using: -``` +```julia julia> Pkg.add("AutoGrad") ``` In order to use it in your code start with: -``` +```julia using AutoGrad ``` @@ -42,7 +42,7 @@ using AutoGrad Here is a linear regression example simplified from [housing.jl](https://github.com/denizyuret/AutoGrad.jl/blob/master/examples/housing.jl): -``` +```julia using AutoGrad function loss(w) @@ -94,7 +94,7 @@ or using the `@primitive` and `@zerograd` macros in [util.jl](https://github.com/denizyuret/AutoGrad.jl/blob/master/src/util.jl) Here is an example: -``` +```julia @primitive hypot(x1::Number,x2::Number),dy,y (dy*x1/y) (dy*x2/y) ``` @@ -111,6 +111,52 @@ example `hypot(x1::Array,x2::Array)` is another hypot method. In AutoGrad.jl each method can independently be defined as a primitive and can have its own specific gradient. +## Jacobian and Higher Order Derivatives +Since `grad` typically returns a differentiable function, it is really easy +to take higher order derivatives of functions with scalar input and scalar output: +```julia +julia> g1 = grad(sin) +(::gradfun) (generic function with 1 method) + +julia> g1(1) == cos(1) +true + +julia> g2 = grad(g1) +(::gradfun) (generic function with 1 method) + +julia> g2(1) == -sin(1) +true +``` +More care has to be taken when functions with vector inputs are considered instead, +since the resulting gradient will be vector valued and AutoGrad is only able to +differentiate scalar valued functions. In this case each output component of the gradient +has to be differentiated separately in order to obtain the Hessian. +AutoGrad provides the utilty function `hessian`, such that + +```julia +julia> A = rand(3, 3); + +julia> f(x) = x'*A*x/2; + +julia> hessian(f)(rand(3)) == (A + A')/2 +true +``` +Similarly the `jacobian` method returns the Jacobian matrix +of a vector valued function: +```julia + +```julia +julia> A = rand(2, 3); + +julia> f(x) = A*x; + +julia> jacobian(f)(rand(3)) == (A + A')/2 +true +``` +Computing Jacobians and Hessians can be heavily resource expensive, +both in time and in memory. In some use cases Jacobian/Hessian product +operators `jvp,vjp,hvp,vhp` can be used as a much faster alternative. + ## Code structure [core.jl](https://github.com/denizyuret/AutoGrad.jl/blob/master/src/core.jl) diff --git a/src/highorder.jl b/src/highorder.jl index 332b7bc..d757a2b 100644 --- a/src/highorder.jl +++ b/src/highorder.jl @@ -59,22 +59,19 @@ It is equivalent to `jacobian(grad(f))`. """ hessian(f) = jacobian(grad(f)) -hvp(f) = jvp(grad(f)) - """ vhp(f) -Returns the vector-Hessian product operat. -It is equivalent to `vjp(grad(f))` +Returns the vector-Hessian product operator. +We have the equivalence `vhp(f)(x, v) == vjp(grad(f))(x, v)`. """ -vhp(f) = vjp(grad(f)) - +vhp(f) = vjp(x->vec(grad(f)(x))) """ hvp(f) Returns the Hessian product operat. -It is equivalent to `jvp(grad(f))` +It is equivalent to `jvp(jacobian(f))` """ -hvp(f) = vjp(grad(f))' # can use vjp instead of jvp since the +hvp(f) = (x,v) -> vhp(f)(x,v)' # can use vjp instead of jvp since the # hessian is symmetric \ No newline at end of file diff --git a/test/highorder.jl b/test/highorder.jl index a0d7566..b5892ae 100644 --- a/test/highorder.jl +++ b/test/highorder.jl @@ -11,21 +11,32 @@ include("header.jl") g8 = grad(g7); @test g8(1)==sin(1) g9 = grad(g8); @test g9(1)==cos(1) - - A = Symmetric(rand(3, 3)) + sym(A) = (A + A')/2 + + A = sym(rand(3, 3)) f(x) = x'*A*x/2 - hessian(f)(rand(3)) == A - hvp(f)(rand(3), v) == A*v - u = rand(2) - vhp(f)(rand(3), u) == u'A + @test hessian(f)(rand(3)) == A + v = rand(3) + @test hvp(f)(rand(3), v) ≈ A*v + u = rand(3) + @test vhp(f)(rand(3), u) ≈ u'A A = rand(2, 3) f(x) = A*x - jacobian(f)(rand(3)) == A + @test jacobian(f)(rand(3)) ≈ A + + A = rand(3, 3) + f(x) = x'*A*x/2 + a = rand(3) + @test jacobian(f)(a) ≈ a' * sym(A) + @test hessian(f)(a) == sym(A) + + A = rand(2, 3) + f(x) = A*x v = rand(3) - jvp(f)(rand(3), v) == A*v + @test jvp(f)(rand(3), v) == A*v u = rand(2) - vjp(f)(rand(3), u) == u'A + @test vjp(f)(rand(3), u) == u'A end nothing