Merge pull request #34 from JuliaAlgebra/st/hessian

Add Hessian

.travis.yml

@@ -4,8 +4,8 @@ os:
   - linux
   # - osx
 julia:
-  - 0.7
   - 1.0
+  - 1.1
   - nightly
 notifications:
   email: false

Project.toml

@@ -2,7 +2,7 @@ name = "StaticPolynomials"
 uuid = "62e018b1-6e46-5407-a5a7-97d4fbcae734"
 license = "MIT"
 repo = "https://github.com/JuliaAlgebra/StaticPolynomials.jl.git"
-version = "1.1.1"
+version = "1.2.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/Project.toml

@@ -3,4 +3,4 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 StaticPolynomials = "62e018b1-6e46-5407-a5a7-97d4fbcae734"
 
 [compat]
-Documenter = "^0.19.6"
+Documenter = "^0.22.0"

docs/make.jl

@@ -1,7 +1,6 @@
 using Documenter, StaticPolynomials
 
 makedocs(
-    format = :html,
     sitename = "StaticPolynomials.jl",
     pages = [
         "Introduction" => "index.md",
@@ -11,9 +10,4 @@ makedocs(
 
 deploydocs(
     repo   = "github.com/JuliaAlgebra/StaticPolynomials.jl.git",
-    target = "build",
-    julia = "0.7",
-    osname = "linux",
-    deps   = nothing,
-    make   = nothing
 )

docs/src/reference.md

@@ -16,6 +16,10 @@ gradient
 gradient!
 evaluate_and_gradient
 evaluate_and_gradient!
+hessian(f::Polynomial, xp...)
+hessian!(U, f::Polynomial, xp...)
+gradient_and_hessian
+gradient_and_hessian!
 differentiate_parameters(::Polynomial, x, p)
 differentiate_parameters!(u, ::Polynomial, x, p)
 ```
@@ -38,6 +42,8 @@ jacobian
 jacobian!
 evaluate_and_jacobian
 evaluate_and_jacobian!
+hessian!(U, F::PolynomialSystem, x::AbstractVector)
+jacobian_and_hessian!
 differentiate_parameters(::PolynomialSystem, x, p)
 differentiate_parameters!(U, ::PolynomialSystem, x, p)
 ```

src/evaluate_codegen.jl

@@ -6,38 +6,38 @@ This assumes that E is in reverse lexicographic order.
 """
 function generate_evaluate(E, ::Type{T}, access_input) where T
     exprs = []
-    last_expr = generate_evaluate!(exprs, E, T, size(E, 1), 1, 1:size(E, 2), access_input)
+    coeffs = map(i -> :(c[$i]), 1:size(E, 2))
+    last_expr = generate_evaluate!(exprs, E, T, size(E, 1), 1, coeffs, access_input)
     Expr(:block, exprs..., last_expr)
 end
 
-function generate_evaluate!(exprs, E, ::Type{T}, nvar, nterm, coeffperm, access_input) where T
+function generate_evaluate!(exprs, E, ::Type{T}, nvar, nterm, coeffs, access_input) where T
     m, n = size(E)
 
     if n == 1
-        return first(monomial_product(T, E[:,1], :(c[$(coeffperm[nterm])]), access_input=access_input))
+        return first(monomial_product(T, E[:,1], coeffs[nterm], access_input=access_input))
     end
 
     if m == 1
-        coeffs = [:(c[$(coeffperm[j])]) for j=nterm:nterm+n-1]
-        return evalpoly(T, E[1,:], coeffs, access_input(nvar))
+        return evalpoly(T, E[1,:], coeffs[nterm:nterm+n-1], access_input(nvar))
     end
 
     # Recursive
-    coeffs = []
+    rec_coeffs = []
     sub_nterm = nterm
     degrees, submatrices = degrees_submatrices(E)
     for E_d in submatrices
-        coeff_subexpr = generate_evaluate!(exprs, E_d, T, nvar - 1, sub_nterm, coeffperm, access_input)
+        coeff_subexpr = generate_evaluate!(exprs, E_d, T, nvar - 1, sub_nterm, coeffs, access_input)
         # If the returned value is just a symbol propagate it
         if isa(coeff_subexpr, Symbol)
-            push!(coeffs, coeff_subexpr)
+            push!(rec_coeffs, coeff_subexpr)
         else
             @gensym c
             push!(exprs, :($c = $coeff_subexpr))
-            push!(coeffs, c)
+            push!(rec_coeffs, c)
         end
         sub_nterm += size(E_d, 2)
     end
 
-    return evalpoly(T, degrees, coeffs, access_input(nvar))
+    return evalpoly(T, degrees, rec_coeffs, access_input(nvar))
 end

src/evaluation.jl

@@ -1,6 +1,7 @@
 export evaluate, gradient, gradient!,
     evaluate_and_gradient, evaluate_and_gradient!,
-    differentiate_parameters, differentiate_parameters!
+    differentiate_parameters, differentiate_parameters!,
+    hessian!, hessian, gradient_and_hessian, gradient_and_hessian!
 
 import Base: @propagate_inbounds
 
@@ -178,6 +179,133 @@ end
     val
 end
 
+
+
+#############
+## Hessian ##
+#############
+
+function _gradient_hessian_impl(f::Type{Polynomial{T, E, P}}; return_grad=true) where {T, E, P}
+    if P == Nothing
+        access_input = x_
+    else
+        n = size(P, 1)
+        access_input(i) = i ≤ n ? :(p[$i]) : :(x[$(i-n)])
+    end
+
+    exps = exponents(E)
+    params = exponents(P)
+
+    n, m = size(E)
+
+    code = []
+    ∇ = [Symbol("∇", i) for i=1:n]
+    vals = [Symbol("val", i) for i=1:n]
+
+
+    # for each variable we simulate the derivative
+    for i in 1:n
+        coeffsᵢ = Union{Symbol, Expr}[]
+        expsᵢ = Vector{eltype(exps)}()
+        paramsᵢ = P === nothing ? nothing : Vector{eltype(exps)}()
+        mᵢ = 0
+
+        # For each term make the symbolic derivative
+        for j in 1:m
+            if exps[i,j] == 0
+                continue
+            end
+            for k in 1:n
+                if i == k
+                    push!(expsᵢ, exps[k, j] - 1)
+                    if exps[i,j] > 1
+                        push!(coeffsᵢ, :($(exps[i,j]) * c[$j]))
+                    else
+                        push!(coeffsᵢ, :(c[$j]))
+                    end
+                else
+                    push!(expsᵢ, exps[k, j])
+                end
+            end
+            if params !== nothing
+                for k in 1:size(params, 1)
+                    push!(paramsᵢ, params[k, j])
+                end
+            end
+
+            mᵢ += 1
+        end
+        if mᵢ == 0
+            zero_tuple = Expr(:tuple, (:(zero($T)) for _=1:n)...)
+            val_grad_codeᵢ = :((zero($T), SVector($zero_tuple)))
+        else
+            Eᵢ = reshape(expsᵢ, (n, mᵢ))
+            Pᵢ = P === nothing ? nothing : reshape(paramsᵢ, (size(paramsᵢ, 1), mᵢ))
+            val_grad_codeᵢ = generate_gradient(Eᵢ, Pᵢ, T, access_input, coeffsᵢ)
+        end
+
+        push!(code, :($(Expr(:tuple, vals[i], ∇[i])) = begin $(val_grad_codeᵢ) end))
+    end
+
+
+    push!(code, :(grad = SVector($(Expr(:tuple, vals...)))))
+    push!(code, :(hess = assemble_matrix(SVector($(Expr(:tuple, (∇[i] for i=1:n)...))))))
+    quote
+        @boundscheck length(x) ≥ $(size(E)[1])
+        c = coefficients(f)
+        $(code...)
+        grad, hess
+    end
+end
+
+@generated function gradient_and_hessian(f::Polynomial{T,E,Nothing}, x::AbstractVector) where {T, E}
+    _gradient_hessian_impl(f)
+end
+@generated function gradient_and_hessian(f::Polynomial, x::AbstractVector, p)
+    _gradient_hessian_impl(f)
+end
+
+@doc """
+     gradient_and_hessian(F::PolynomialSystem, x, [p])
+
+Compute the gradient and Hessian of `f` at `x`.
+""" gradient_and_hessian(F::Polynomial, x)
+
+"""
+    gradient_and_hessian!(u, U, f::Polynomial, x, [p])
+
+Compute the gradient and Hessian of `f` at `x` and store them in `u` and `U`.
+"""
+function gradient_and_hessian!(u, U, f::Polynomial, xp...)
+    grad, hess = gradient_and_hessian(f, xp...)
+    u .= grad
+    U .= hess
+    nothing
+end
+
+"""
+    hessian(f::Polynomial, x, [p])
+
+Compute the hessian of `f` at `x`.
+"""
+function hessian(f::Polynomial, xp...)
+    last(gradient_and_hessian(f, xp...))
+end
+
+"""
+    hessian!(U, f::Polynomial, x, [p])
+
+Compute the hessian of `f` at `x` and store it in `U`.
+"""
+function hessian!(U, f::Polynomial, xp...)
+    U .= hessian(f, xp...)
+    U
+end
+
+
+################
+## Parameters ##
+################
 function _differentiate_parameters_impl(f::Type{Polynomial{T, E, P}}) where {T, E, P}
     @assert P != Nothing