From 491d2aa0b756fd188c1f381c862e317e696a234e Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Wed, 14 Dec 2016 09:32:04 +0100 Subject: [PATCH] Return lower order results (#94) * return value for gradient * update tests * add docstrings and export gradient and hessian * remove unnecessary qualifiers after exporting gradient [ci skip] * fix docstring, write something about all in docs * make examples a header --- docs/src/demos.md | 2 +- docs/src/man/automatic_differentiation.md | 27 +- src/ContMechTensors.jl | 1 + src/automatic_differentiation.jl | 301 +++++++++++++++------- test/test_ad.jl | 52 ++-- 5 files changed, 268 insertions(+), 115 deletions(-) diff --git a/docs/src/demos.md b/docs/src/demos.md index 6b9e35e..91ca467 100644 --- a/docs/src/demos.md +++ b/docs/src/demos.md @@ -105,7 +105,7 @@ julia> F = one(Tensor{2,3}) + rand(Tensor{2,3}); julia> C = tdot(F); -julia> S_AD = 2 * ContMechTensors.gradient(C -> Ψ(C, μ, Kb), C) +julia> S_AD = 2 * gradient(C -> Ψ(C, μ, Kb), C) 3×3 ContMechTensors.SymmetricTensor{2,3,Float64,6}: 4.30534e11 -2.30282e11 -8.52861e10 -2.30282e11 4.38793e11 -2.64481e11 diff --git a/docs/src/man/automatic_differentiation.md b/docs/src/man/automatic_differentiation.md index dc6aba7..38c35fe 100644 --- a/docs/src/man/automatic_differentiation.md +++ b/docs/src/man/automatic_differentiation.md @@ -7,24 +7,37 @@ end # Automatic Differentiation +```@index +Pages = ["automatic_differentiation.md"] +``` + `ContMechTensors` supports forward mode automatic differentiation (AD) of tensorial functions to compute first order derivatives (gradients) and second order derivatives (Hessians). It does this by exploiting the `Dual` number defined in `ForwardDiff.jl`. While `ForwardDiff.jl` can itself be used to differentiate tensor functions it is a bit awkward because `ForwardDiff.jl` is written to work with standard Julia `Array`s. One therefore has to send the input argument as an `Array` to `ForwardDiff.jl`, convert it to a `Tensor` and then convert the output `Array` to a `Tensor` again. This can also be inefficient since these `Array`s are allocated on the heap so one needs to preallocate which can be annoying. Instead, it is simpler to use `ContMechTensors` own AD API to do the differentiation. This does not require any conversions and everything will be stack allocated so there is no need to preallocate. -The API for AD in `ContMechTensors` is `ContMechTensors.gradient(f, A)` and `ContMechTensors.hessian(f, A)` where `f` is a function and `A` is a first or second order tensor. For `gradient` the function can return a scalar, vector (in case the input is a vector) or a second order tensor. For `hessian` the function should return a scalar. +The API for AD in `ContMechTensors` is `gradient(f, A)` and `hessian(f, A)` where `f` is a function and `A` is a first or second order tensor. For `gradient` the function can return a scalar, vector (in case the input is a vector) or a second order tensor. For `hessian` the function should return a scalar. + +When evaluating the function with dual numbers, the value (value and gradient in the case of hessian) is obtained automatically, along with the gradient. To obtain the lower order results `gradient` and `hessian` accepts a third arguement, a `Symbol`. Note that the symbol is only used to dispatch to the correct function, and thus it can be any symbol. In the examples the symbol `:all` is used to obtain all the lower order derivatives and values. + +```@docs +gradient +hessian +``` + +## Examples We here give a few examples of differentiating various functions and compare with the analytical solution. -## Norm of a vector +### Norm of a vector $f(\mathbf{x}) = |\mathbf{x}| \quad \Rightarrow \quad \partial f / \partial \mathbf{x} = \mathbf{x} / |\mathbf{x}|$ ```jldoctest julia> x = rand(Vec{2}); -julia> ContMechTensors.gradient(norm, x) +julia> gradient(norm, x) 2-element ContMechTensors.Tensor{1,2,Float64,2}: 0.61036 0.792124 @@ -35,14 +48,14 @@ julia> x / norm(x) 0.792124 ``` -## Determinant of a second order symmetric tensor +### Determinant of a second order symmetric tensor $f(\mathbf{A}) = \det \mathbf{A} \quad \Rightarrow \quad \partial f / \partial \mathbf{A} = \mathbf{A}^{-T} \det \mathbf{A}$ ```jldoctest julia> A = rand(SymmetricTensor{2,2}); -julia> ContMechTensors.gradient(det, A) +julia> gradient(det, A) 2×2 ContMechTensors.SymmetricTensor{2,2,Float64,3}: 0.566237 -0.766797 -0.766797 0.590845 @@ -53,7 +66,7 @@ julia> inv(A)' * det(A) -0.766797 0.590845 ``` -## Hessian of a quadratic potential +### Hessian of a quadratic potential $\psi(\mathbf{e}) = 1/2 \mathbf{e} : \mathsf{E} : \mathbf{e} \quad \Rightarrow \quad \partial \psi / (\partial \mathbf{e} \otimes \partial \mathbf{e}) = \mathsf{E}^\text{sym}$ @@ -64,7 +77,7 @@ julia> E = rand(Tensor{4,2}); julia> ψ(ϵ) = 1/2 * ϵ ⊡ E ⊡ ϵ; -julia> E_sym = ContMechTensors.hessian(ψ, rand(Tensor{2,2})); +julia> E_sym = hessian(ψ, rand(Tensor{2,2})); julia> norm(majorsymmetric(E) - E_sym) 0.0 diff --git a/src/ContMechTensors.jl b/src/ContMechTensors.jl index 562225a..3707463 100644 --- a/src/ContMechTensors.jl +++ b/src/ContMechTensors.jl @@ -10,6 +10,7 @@ export AbstractTensor, SymmetricTensor, Tensor, Vec, FourthOrderTensor, SecondOr export otimes, ⊗, ⊡, dcontract, dev, vol, symmetric, skew, minorsymmetric, majorsymmetric export minortranspose, majortranspose, isminorsymmetric, ismajorsymmetric export tdot, dotdot +export hessian#, gradient @deprecate extract_components(tensor) Array(tensor) diff --git a/src/automatic_differentiation.jl b/src/automatic_differentiation.jl index 06a6ee7..d75488e 100644 --- a/src/automatic_differentiation.jl +++ b/src/automatic_differentiation.jl @@ -1,14 +1,11 @@ -using ForwardDiff -using ContMechTensors -import ForwardDiff: Dual, partials - +import ForwardDiff: Dual, partials, value ###################### # Extraction methods # ###################### -# Extractions are supposed to unpack the partials in a tensor -# and put it into a tensor of higher order. +# Extractions are supposed to unpack the value and the partials +# The partials should be put into a tensor of higher order. # The extraction methods need to know the input type to the function # that generated the result. The reason for this is that there is no # difference in the output type (except the number of partials) for @@ -16,130 +13,186 @@ import ForwardDiff: Dual, partials # For dim = 1 the type is exactly the same. # Scalar -@inline function _extract{N, T}(v::Dual{N, T}, ::Vec{N}) - Vec{N, T}(partials(v).values) +@inline function _extract_value{N, T}(v::Dual{N, T}, ::Vec{N}) + return value(v) +end +@inline function _extract_gradient{N, T}(v::Dual{N, T}, ::Vec{N}) + return Vec{N, T}(partials(v).values) end -@inline function _extract{N, T, dim, T2}(v::Dual{N, T}, ::SymmetricTensor{2, dim, T2, N}) - SymmetricTensor{2, dim, T}(partials(v).values) +@inline function _extract_value{N, T, dim, T2}(v::Dual{N, T}, ::SymmetricTensor{2, dim, T2, N}) + return value(v) +end +@inline function _extract_gradient{N, T, dim, T2}(v::Dual{N, T}, ::SymmetricTensor{2, dim, T2, N}) + return SymmetricTensor{2, dim, T}(partials(v).values) end -@inline function _extract{N, T, dim, T2}(v::Dual{N, T}, ::Tensor{2, dim, T2, N}) - Tensor{2, dim, T}(partials(v).values) +@inline function _extract_value{N, T, dim, T2}(v::Dual{N, T}, ::Tensor{2, dim, T2, N}) + return value(v) +end +@inline function _extract_gradient{N, T, dim, T2}(v::Dual{N, T}, ::Tensor{2, dim, T2, N}) + return Tensor{2, dim, T}(partials(v).values) end # Vec -@inline function _extract{D <: Dual}(v::Vec{1, D}, ::Any) +@inline function _extract_value{D <: Dual}(v::Vec{1, D}, ::Any) + @inbounds begin + v1 = value(v[1]) + f = Vec{1}((v1,)) + end + return f +end +@inline function _extract_gradient{D <: Dual}(v::Vec{1, D}, ::Any) @inbounds begin p1 = partials(v[1]) - d = Tensor{2, 1}((p1[1],)) + ∇f = Tensor{2, 1}((p1[1],)) end - return d + return ∇f end -@inline function _extract{D <: Dual}(v::Vec{2, D}, ::Any) +@inline function _extract_value{D <: Dual}(v::Vec{2, D}, ::Any) @inbounds begin - p1 = partials(v[1]) - p2 = partials(v[2]) - d = Tensor{2, 2}((p1[1], p2[1], p1[2], p2[2])) + v1, v2 = value(v[1]), value(v[2]) + f = Vec{2}((v1, v2)) end - return d + return f +end +@inline function _extract_gradient{D <: Dual}(v::Vec{2, D}, ::Any) + @inbounds begin + p1, p2 = partials(v[1]), partials(v[2]) + ∇f = Tensor{2, 2}((p1[1], p2[1], p1[2], p2[2])) + end + return ∇f end -@inline function _extract{D <: Dual}(v::Vec{3, D}, ::Any) +@inline function _extract_value{D <: Dual}(v::Vec{3, D}, ::Any) @inbounds begin - p1 = partials(v[1]) - p2 = partials(v[2]) - p3 = partials(v[3]) - d = Tensor{2, 3}((p1[1], p2[1], p3[1], p1[2], p2[2], p3[2], p1[3], p2[3], p3[3])) + v1, v2, v3 = value(v[1]), value(v[2]), value(v[3]) + f = Vec{3}((v1, v2, v3)) + end + return f +end +@inline function _extract_gradient{D <: Dual}(v::Vec{3, D}, ::Any) + @inbounds begin + p1, p2, p3 = partials(v[1]), partials(v[2]), partials(v[3]) + ∇f = Tensor{2, 3}((p1[1], p2[1], p3[1], p1[2], p2[2], p3[2], p1[3], p2[3], p3[3])) end - return d + return ∇f end # Second order tensor -@inline function _extract{D <: Dual}(v::Tensor{2, 1, D}, ::Any) +@inline function _extract_value{D <: Dual}(v::Tensor{2, 1, D}, ::Any) + @inbounds begin + v1 = value(v[1,1]) + f = Tensor{2, 1}((v1,)) + end + return f +end +@inline function _extract_gradient{D <: Dual}(v::Tensor{2, 1, D}, ::Any) @inbounds begin p1 = partials(v[1,1]) - d = Tensor{4, 1}((p1[1],)) + ∇f = Tensor{4, 1}((p1[1],)) end - return d + return ∇f end -@inline function _extract{D <: Dual}(v::SymmetricTensor{2, 1, D}, ::Any) +@inline function _extract_value{D <: Dual}(v::SymmetricTensor{2, 1, D}, ::Any) + @inbounds begin + v1 = value(v[1,1]) + f = SymmetricTensor{2, 1}((v1,)) + end + return f +end +@inline function _extract_gradient{D <: Dual}(v::SymmetricTensor{2, 1, D}, ::Any) @inbounds begin p1 = partials(v[1,1]) - d = SymmetricTensor{4, 1}((p1[1],)) + ∇f = SymmetricTensor{4, 1}((p1[1],)) end - return d + return ∇f end -@inline function _extract{D <: Dual}(v::Tensor{2, 2, D}, ::Any) +@inline function _extract_value{D <: Dual}(v::Tensor{2, 2, D}, ::Any) @inbounds begin - p1 = partials(v[1,1]) - p2 = partials(v[2,1]) - p3 = partials(v[1,2]) - p4 = partials(v[2,2]) - d = Tensor{4, 2}((p1[1], p2[1], p3[1], p4[1], - p1[2], p2[2], p3[2], p4[2], - p1[3], p2[3], p3[3], p4[3], - p1[4], p2[4], p3[4], p4[4])) + v1, v2, v3, v4 = value(v[1,1]), value(v[2,1]), value(v[1,2]), value(v[2,2]) + f = Tensor{2, 2}((v1, v2, v3, v4)) end - return d + return f +end +@inline function _extract_gradient{D <: Dual}(v::Tensor{2, 2, D}, ::Any) + @inbounds begin + p1, p2, p3, p4 = partials(v[1,1]), partials(v[2,1]), partials(v[1,2]), partials(v[2,2]) + ∇f = Tensor{4, 2}((p1[1], p2[1], p3[1], p4[1], + p1[2], p2[2], p3[2], p4[2], + p1[3], p2[3], p3[3], p4[3], + p1[4], p2[4], p3[4], p4[4])) + end + return ∇f end -@inline function _extract{D <: Dual}(v::SymmetricTensor{2, 2, D}, ::Any) +@inline function _extract_value{D <: Dual}(v::SymmetricTensor{2, 2, D}, ::Any) @inbounds begin - p1 = partials(v[1,1]) - p2 = partials(v[2,1]) - p3 = partials(v[2,2]) - d = SymmetricTensor{4, 2}((p1[1], p2[1], p3[1], - p1[2], p2[2], p3[2], - p1[3], p2[3], p3[3])) + v1, v2, v3 = value(v[1,1]), value(v[2,1]), value(v[2,2]) + f = SymmetricTensor{2, 2}((v1, v2, v3)) + end + return f +end +@inline function _extract_gradient{D <: Dual}(v::SymmetricTensor{2, 2, D}, ::Any) + @inbounds begin + p1, p2, p3 = partials(v[1,1]), partials(v[2,1]), partials(v[2,2]) + ∇f = SymmetricTensor{4, 2}((p1[1], p2[1], p3[1], + p1[2], p2[2], p3[2], + p1[3], p2[3], p3[3])) end - return d + return ∇f end -@inline function _extract{D <: Dual}(v::Tensor{2, 3, D}, ::Any) +@inline function _extract_value{D <: Dual}(v::Tensor{2, 3, D}, ::Any) @inbounds begin - p1 = partials(v[1,1]) - p2 = partials(v[2,1]) - p3 = partials(v[3,1]) - p4 = partials(v[1,2]) - p5 = partials(v[2,2]) - p6 = partials(v[3,2]) - p7 = partials(v[1,3]) - p8 = partials(v[2,3]) - p9 = partials(v[3,3]) - d = Tensor{4, 3}((p1[1], p2[1], p3[1], p4[1], p5[1], p6[1], p7[1], p8[1], p9[1], - p1[2], p2[2], p3[2], p4[2], p5[2], p6[2], p7[2], p8[2], p9[2], # ### # - p1[3], p2[3], p3[3], p4[3], p5[3], p6[3], p7[3], p8[3], p9[3], # # # # - p1[4], p2[4], p3[4], p4[4], p5[4], p6[4], p7[4], p8[4], p9[4], ### ### ### - p1[5], p2[5], p3[5], p4[5], p5[5], p6[5], p7[5], p8[5], p9[5], - p1[6], p2[6], p3[6], p4[6], p5[6], p6[6], p7[6], p8[6], p9[6], - p1[7], p2[7], p3[7], p4[7], p5[7], p6[7], p7[7], p8[7], p9[7], - p1[8], p2[8], p3[8], p4[8], p5[8], p6[8], p7[8], p8[8], p9[8], - p1[9], p2[9], p3[9], p4[9], p5[9], p6[9], p7[9], p8[9], p9[9])) - end - return d -end - -@inline function _extract{D <: Dual}(v::SymmetricTensor{2, 3, D}, ::Any) + v1, v2, v3 = value(v[1,1]), value(v[2,1]), value(v[3,1]) + v4, v5, v6 = value(v[1,2]), value(v[2,2]), value(v[3,2]) + v7, v8, v9 = value(v[1,3]), value(v[2,3]), value(v[3,3]) + f = Tensor{2,3}((v1, v2, v3, v4, v5, v6, v7, v8, v9)) + end + return f +end +@inline function _extract_gradient{D <: Dual}(v::Tensor{2, 3, D}, ::Any) @inbounds begin - p1 = partials(v[1,1]) - p2 = partials(v[2,1]) - p3 = partials(v[3,1]) - p4 = partials(v[2,2]) - p5 = partials(v[3,2]) - p6 = partials(v[3,3]) + p1, p2, p3 = partials(v[1,1]), partials(v[2,1]), partials(v[3,1]) + p4, p5, p6 = partials(v[1,2]), partials(v[2,2]), partials(v[3,2]) + p7, p8, p9 = partials(v[1,3]), partials(v[2,3]), partials(v[3,3]) + ∇f = Tensor{4, 3}((p1[1], p2[1], p3[1], p4[1], p5[1], p6[1], p7[1], p8[1], p9[1], + p1[2], p2[2], p3[2], p4[2], p5[2], p6[2], p7[2], p8[2], p9[2], # ### # + p1[3], p2[3], p3[3], p4[3], p5[3], p6[3], p7[3], p8[3], p9[3], # # # # + p1[4], p2[4], p3[4], p4[4], p5[4], p6[4], p7[4], p8[4], p9[4], ### ### ### + p1[5], p2[5], p3[5], p4[5], p5[5], p6[5], p7[5], p8[5], p9[5], + p1[6], p2[6], p3[6], p4[6], p5[6], p6[6], p7[6], p8[6], p9[6], + p1[7], p2[7], p3[7], p4[7], p5[7], p6[7], p7[7], p8[7], p9[7], + p1[8], p2[8], p3[8], p4[8], p5[8], p6[8], p7[8], p8[8], p9[8], + p1[9], p2[9], p3[9], p4[9], p5[9], p6[9], p7[9], p8[9], p9[9])) + end + return ∇f +end - d = SymmetricTensor{4, 3}((p1[1], p2[1], p3[1], p4[1], p5[1], p6[1], - p1[2], p2[2], p3[2], p4[2], p5[2], p6[2], - p1[3], p2[3], p3[3], p4[3], p5[3], p6[3], - p1[4], p2[4], p3[4], p4[4], p5[4], p6[4], - p1[5], p2[5], p3[5], p4[5], p5[5], p6[5], - p1[6], p2[6], p3[6], p4[6], p5[6], p6[6])) +@inline function _extract_value{D <: Dual}(v::SymmetricTensor{2, 3, D}, ::Any) + @inbounds begin + v1, v2, v3 = value(v[1,1]), value(v[2,1]), value(v[3,1]) + v4, v5, v6 = value(v[2,2]), value(v[3,2]), value(v[3,3]) + f = SymmetricTensor{2, 3}((v1, v2, v3, v4, v5, v6)) end - return d + return f +end +@inline function _extract_gradient{D <: Dual}(v::SymmetricTensor{2, 3, D}, ::Any) + @inbounds begin + p1, p2, p3 = partials(v[1,1]), partials(v[2,1]), partials(v[3,1]) + p4, p5, p6 = partials(v[2,2]), partials(v[3,2]), partials(v[3,3]) + ∇f = SymmetricTensor{4, 3}((p1[1], p2[1], p3[1], p4[1], p5[1], p6[1], + p1[2], p2[2], p3[2], p4[2], p5[2], p6[2], + p1[3], p2[3], p3[3], p4[3], p5[3], p6[3], + p1[4], p2[4], p3[4], p4[4], p5[4], p6[4], + p1[5], p2[5], p3[5], p4[5], p5[5], p6[5], + p1[6], p2[6], p3[6], p4[6], p5[6], p6[6])) + end + return ∇f end ################## @@ -147,7 +200,7 @@ end ################## # Loaders are supposed to take a tensor of real values and convert it -# into a tensor of dual values where the seeds are correctly defind. +# into a tensor of dual values where the seeds are correctly defined. @inline function _load{T}(v::Vec{1, T}) @inbounds v_dual = Vec{1}((Dual(v[1], one(T)),)) @@ -234,13 +287,81 @@ end return v_dual end -function gradient{F}(f::F, v::Union{SecondOrderTensor, Vec}) +""" +Computes the gradient of the input function. If the (pseudo)-keyword `all` +is given, the value of the function is also returned as a second output argument + +```julia +gradient(f::Function, v::Union{SecondOrderTensor, Vec}) +gradient(f::Function, v::Union{SecondOrderTensor, Vec}, :all) +``` + +**Example:** + +```jldoctest +julia> A = rand(SymmetricTensor{2, 2}); + +julia> ∇f = gradient(norm, A) +2×2 ContMechTensors.SymmetricTensor{2,2,Float64,3}: + 0.434906 0.56442 + 0.56442 0.416793 + +julia> ∇f, f = gradient(norm, A, :all); +``` +""" +function Base.gradient{F}(f::F, v::Union{SecondOrderTensor, Vec}) v_dual = _load(v) res = f(v_dual) - return _extract(res, v) + return _extract_gradient(res, v) end +function Base.gradient{F}(f::F, v::Union{SecondOrderTensor, Vec}, ::Symbol) + v_dual = _load(v) + res = f(v_dual) + return _extract_gradient(res, v), _extract_value(res, v) +end + +""" +Computes the hessian of the input function. If the (pseudo)-keyword `all` +is given, the lower order results (gradient and value) of the function is +also returned as a second and third output argument. + +```julia +hessian(f::Function, v::Union{SecondOrderTensor, Vec}) +hessian(f::Function, v::Union{SecondOrderTensor, Vec}, :all) +``` + +**Example:** + +```jldoctest +julia> A = rand(SymmetricTensor{2, 2}); +julia> ∇∇f = hessian(norm, A) +2×2×2×2 ContMechTensors.SymmetricTensor{4,2,Float64,9}: +[:, :, 1, 1] = + 0.596851 -0.180684 + -0.180684 -0.133425 + +[:, :, 2, 1] = + -0.180684 0.133546 + 0.133546 -0.173159 + +[:, :, 1, 2] = + -0.180684 0.133546 + 0.133546 -0.173159 + +[:, :, 2, 2] = + -0.133425 -0.173159 + -0.173159 0.608207 + +julia> ∇∇f, ∇f, f = hessian(norm, A, :all); +``` +""" function hessian{F}(f::F, v::Union{SecondOrderTensor, Vec}) gradf = y -> gradient(f, y) return gradient(gradf, v) end + +function hessian{F}(f::F, v::Union{SecondOrderTensor, Vec}, ::Symbol) + gradf = y -> gradient(f, y) + return gradient(gradf, v), gradient(f, v, :all)... +end diff --git a/test/test_ad.jl b/test/test_ad.jl index f1ab371..fadbd2d 100644 --- a/test/test_ad.jl +++ b/test/test_ad.jl @@ -49,32 +49,50 @@ for dim in 1:3 II_sym = one(SymmetricTensor{4, dim, T}) # Gradient of scalars - @test ∇(norm, v) ≈ v / norm(v) - @test ∇(norm, A) ≈ A / norm(A) - @test ∇(norm, A_sym) ≈ A_sym / norm(A_sym) - @test ∇(v -> 3*v, v) ≈ diagm(Tensor{2, dim}, 3.0) + @test ∇(norm, v) ≈ ∇(norm, v, :all)[1] ≈ v / norm(v) + @test ∇(norm, v, :all)[2] ≈ norm(v) + @test ∇(norm, A) ≈ ∇(norm, A, :all)[1] ≈ A / norm(A) + @test ∇(norm, A, :all)[2] ≈ norm(A) + @test ∇(norm, A_sym) ≈ ∇(norm, A_sym, :all)[1] ≈ A_sym / norm(A_sym) + @test ∇(norm, A_sym, :all)[2] ≈ norm(A_sym) + @test ∇(v -> 3*v, v) ≈ ∇(v -> 3*v, v, :all)[1] ≈ diagm(Tensor{2, dim}, 3.0) + @test ∇(v -> 3*v, v, :all)[2] ≈ 3*v # https://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivatives_of_the_invariants_of_a_second-order_tensor I1, DI1 = A -> trace(A), A -> one(A) I2, DI2 = A -> 1/2 * (trace(A)^2 - trace(A⋅A)), A -> I1(A) * one(A) - A' I3, DI3 = A -> det(A), A -> det(A) * inv(A)' - @test ∇(I1, A) ≈ DI1(A) - @test ∇(I2, A) ≈ DI2(A) - @test ∇(I3, A) ≈ DI3(A) - @test ∇(I1, A_sym) ≈ DI1(A_sym) - @test ∇(I2, A_sym) ≈ DI2(A_sym) - @test ∇(I3, A_sym) ≈ DI3(A_sym) + @test ∇(I1, A) ≈ ∇(I1, A, :all)[1] ≈ DI1(A) + @test ∇(I1, A, :all)[2] ≈ I1(A) + @test ∇(I2, A) ≈ ∇(I2, A, :all)[1] ≈ DI2(A) + @test ∇(I2, A, :all)[2] ≈ I2(A) + @test ∇(I3, A) ≈ ∇(I3, A, :all)[1] ≈ DI3(A) + @test ∇(I3, A, :all)[2] ≈ I3(A) + @test ∇(I1, A_sym) ≈ ∇(I1, A_sym, :all)[1] ≈ DI1(A_sym) + @test ∇(I1, A_sym, :all)[2] ≈ I1(A_sym) + @test ∇(I2, A_sym) ≈ ∇(I2, A_sym, :all)[1] ≈ DI2(A_sym) + @test ∇(I2, A_sym, :all)[2] ≈ I2(A_sym) + @test ∇(I3, A_sym) ≈ ∇(I3, A_sym, :all)[1] ≈ DI3(A_sym) + @test ∇(I3, A_sym, :all)[2] ≈ I3(A_sym) # Gradient of second order tensors - @test ∇(identity, A) ≈ II - @test ∇(identity, A_sym) ≈ II_sym - @test ∇(transpose, A) ⊡ B ≈ B' - @test ∇(transpose, A_sym) ⊡ B_sym ≈ B_sym' - @test ∇(inv, A) ⊡ B ≈ - inv(A) ⋅ B ⋅ inv(A) - @test ∇(inv, A_sym) ⊡ B_sym ≈ - inv(A_sym) ⋅ B_sym ⋅ inv(A_sym) + @test ∇(identity, A) ≈ ∇(identity, A, :all)[1] ≈ II + @test ∇(identity, A, :all)[2] ≈ A + @test ∇(identity, A_sym) ≈ ∇(identity, A_sym, :all)[1] ≈ II_sym + @test ∇(identity, A_sym, :all)[2] ≈ A_sym + @test ∇(transpose, A) ⊡ B ≈ ∇(transpose, A, :all)[1] ⊡ B ≈ B' + @test ∇(transpose, A, :all)[2] ≈ A' + @test ∇(transpose, A_sym) ⊡ B_sym ≈ ∇(transpose, A_sym, :all)[1] ⊡ B_sym ≈ B_sym' + @test ∇(transpose, A_sym, :all)[2] ≈ A_sym' + @test ∇(inv, A) ⊡ B ≈ ∇(inv, A, :all)[1] ⊡ B ≈ - inv(A) ⋅ B ⋅ inv(A) + @test ∇(inv, A, :all)[2] ≈ inv(A) + @test ∇(inv, A_sym) ⊡ B_sym ≈ ∇(inv, A_sym, :all)[1] ⊡ B_sym ≈ - inv(A_sym) ⋅ B_sym ⋅ inv(A_sym) + @test ∇(inv, A_sym, :all)[2] ≈ inv(A_sym) # Hessians of scalars - @test Δ(norm, A).data ≈ vec(ForwardDiff.hessian(x -> sqrt(sumabs2(x)), A.data)) + @test Δ(norm, A).data ≈ Δ(norm, A, :all)[1].data ≈ vec(ForwardDiff.hessian(x -> sqrt(sumabs2(x)), A.data)) + @test Δ(norm, A, :all)[2].data ≈ vec(ForwardDiff.gradient(x -> sqrt(sumabs2(x)), A.data)) + @test Δ(norm, A, :all)[3] ≈ norm(A) end end end # testset