codegen - Maple Programming Help

Home : Support : Online Help : Programming : codegen Package : codegen/HESSIAN

codegen

 HESSIAN
 compute the HESSIAN matrix of a Maple procedure

 Calling Sequence HESSIAN(F) HESSIAN(F, X) HESSIAN(F, X, ...)

Parameters

 F - Maple procedure X - list of symbols (parameters of F)

Description

 • The first argument F is a Maple procedure which computes a function of x1,x2,...,xn.  The HESSIAN command outputs a new procedure H, which when executed at given values for x1,x2,...,xn, returns a matrix of the second partial derivatives of H w.r.t. x1,...,xn at the given values.  For example, given

 F := proc(x, y) local t; t := exp(-x); y*t + t end proc

 The output of H := HESSIAN(F); is the procedure

 H := proc(x, y) local grd1, grd2, df, grd, df1, t, dfr0; t := exp(-x); df1 := y + 1; grd1 := - t*df1; grd2 := t; df := array(1 .. 4); dfr0 := array(1 .. 4); df[3] := 1; df[2] := - df[3]*t; df[1] := - df[3]*df1; dfr0[4] := 1; dfr0[1] := dfr0[4]; grd := array(1 .. 2, 1 .. 2); grd[1, 1] := - df[1]*exp(-x); grd[1, 2] := df[2]; grd[2, 1] := - dfr0[1]*exp(-x); grd[2, 2] := 0; return grd end proc

 The H procedure can be optimized by optimize(H). When H is called with inputs $1.0,1.0$, it outputs the matrix

$\left[\begin{array}{cc}0.7357588824& -0.3678794412\\ -0.3678794412& 0\end{array}\right]$

 • The code in H is constructed by applying the GRADIENT command to F twice.  The GRADIENT command uses automatic differentiation. This often leads to a more efficient computation than symbolic differentiation, that is, what you would obtain from using linalg[hessian]. See codegen[GRADIENT for further details on automatic differentiation. The remaining arguments to HESSIAN are optional, they are described below.
 • By default, HESSIAN computes the partial derivatives of F w.r.t. all the parameters present in F.  The optional argument X, a list of symbols, may be used to specify which parameters to take the derivative w.r.t.
 • Two algorithms are supported, the so-called forward and reverse modes. By default, HESSIAN tries to use the reverse mode since it usually leads to a more efficient code.  If it is unable to use the reverse mode, the forward mode is used.  The user may specify which algorithm is to be used by giving the optional argument mode=forward or mode=reverse.
 • The matrix of partial derivatives is, by default, returned as an array. The optional argument result_type=list, result_type=array, or result_type=seq specifies that the matrix of derivatives returned by H is to be a Maple list, array, and sequence respectively.
 • The command with(codegen,HESSIAN) allows the use of the abbreviated form of this command.

Examples

 > $\mathrm{with}\left(\mathrm{codegen}\right):$
 > F := proc(x,y) local t; t := x*y; x+t-y*t; end proc;
 ${F}{≔}{\mathbf{proc}}\left({x}{,}{y}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{local}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{t}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{t}{≔}{x}{*}{y}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{x}{+}{t}{-}{y}{*}{t}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (1)
 > $F\left(x,y\right)$
 ${-}{x}{}{{y}}^{{2}}{+}{x}{}{y}{+}{x}$ (2)
 > $H≔\mathrm{HESSIAN}\left(F\right)$
 ${H}{≔}{\mathbf{proc}}\left({x}{,}{y}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{local}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df}}{,}{\mathrm{df1}}{,}{\mathrm{dfr0}}{,}{\mathrm{grd}}{,}{\mathrm{grd1}}{,}{\mathrm{grd2}}{,}{\mathrm{t1}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df1}}{≔}{−}{y}{+}{1}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t1}}{≔}{\mathrm{df1}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd1}}{≔}{\mathrm{t1}}{*}{y}{+}{1}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd2}}{≔}{\mathrm{t1}}{*}{x}{-}{x}{*}{y}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df}}{≔}{\mathrm{array}}{}\left({1}{..}{4}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr0}}{≔}{\mathrm{array}}{}\left({1}{..}{4}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df}}{[}{3}{]}{≔}{1}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df}}{[}{2}{]}{≔}{\mathrm{df}}{[}{3}{]}{*}{y}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df}}{[}{1}{]}{≔}{\mathrm{df}}{[}{2}{]}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr0}}{[}{4}{]}{≔}{1}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr0}}{[}{2}{]}{≔}{\mathrm{dfr0}}{[}{4}{]}{*}{x}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr0}}{[}{1}{]}{≔}{\mathrm{dfr0}}{[}{2}{]}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd}}{≔}{\mathrm{array}}{}\left({1}{..}{2}{,}{1}{..}{2}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd}}{[}{1}{,}{1}{]}{≔}{0}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd}}{[}{1}{,}{2}{]}{≔}{\mathrm{t1}}{*}{\mathrm{df}}{[}{3}{]}{-}{\mathrm{df}}{[}{1}{]}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd}}{[}{2}{,}{1}{]}{≔}{\mathrm{dfr0}}{[}{4}{]}{*}\left({\mathrm{t1}}{-}{y}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd}}{[}{2}{,}{2}{]}{≔}{−}{\mathrm{dfr0}}{[}{4}{]}{*}{x}{-}{\mathrm{dfr0}}{[}{1}{]}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{return}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (3)
 > $\mathrm{optimize}\left(H\right)$
 ${\mathbf{proc}}\left({x}{,}{y}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{local}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df1}}{,}{\mathrm{grd}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df1}}{≔}{−}{y}{+}{1}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd}}{≔}{\mathrm{array}}{}\left({1}{..}{2}{,}{1}{..}{2}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd}}{[}{1}{,}{1}{]}{≔}{0}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd}}{[}{1}{,}{2}{]}{≔}{\mathrm{df1}}{-}{y}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd}}{[}{2}{,}{1}{]}{≔}{\mathrm{grd}}{[}{1}{,}{2}{]}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd}}{[}{2}{,}{2}{]}{≔}{−}{2}{*}{x}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{grd}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (4)
 > $\mathrm{H1}≔H\left(x,y\right)$
 ${\mathrm{H1}}{≔}{\mathrm{grd}}$ (5)
 > $\mathrm{eval}\left(\mathrm{H1}\right)$
 $\left[\begin{array}{cc}{0}& {-}{2}{}{y}{+}{1}\\ {-}{2}{}{y}{+}{1}& {-}{2}{}{x}\end{array}\right]$ (6)
 > $\mathrm{HESSIAN}\left(F,\mathrm{mode}=\mathrm{forward},\mathrm{result_type}=\mathrm{seq}\right)$
 ${\mathbf{proc}}\left({x}{,}{y}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{local}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dt1}}{,}{\mathrm{t1}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dt1}}{≔}{\mathrm{array}}{}\left({1}{..}{2}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dt1}}{[}{1}{]}{≔}{0}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dt1}}{[}{2}{]}{≔}{−}{1}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t1}}{≔}{−}{y}{+}{1}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{return}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{y}{*}{\mathrm{dt1}}{[}{1}{]}{,}{y}{*}{\mathrm{dt1}}{[}{2}{]}{+}{\mathrm{t1}}{,}{x}{*}{\mathrm{dt1}}{[}{1}{]}{+}{\mathrm{t1}}{-}{y}{,}{x}{*}{\mathrm{dt1}}{[}{2}{]}{-}{x}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (7)

This example we compute the Hessian w.r.t. the first two parameters phi and omega only.  Since the torus program returns a vector of values, the result is of dimension 3.

 > torus  := proc(phi,omega,R,r) local x,y,z;      x := cos(phi)*(R+r*cos(omega));      y := sin(phi)*(R+r*cos(omega));      z := r*sin(omega);      [x,y,z] end proc:
 > $\mathrm{torus}≔\mathrm{optimize}\left(\mathrm{torus}\right)$
 ${\mathrm{torus}}{≔}{\mathbf{proc}}\left({\mathrm{φ}}{,}{\mathrm{ω}}{,}{R}{,}{r}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{local}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t1}}{,}{\mathrm{t2}}{,}{\mathrm{t4}}{,}{\mathrm{t5}}{,}{\mathrm{t6}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t1}}{≔}{\mathrm{cos}}{}\left({\mathrm{φ}}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t2}}{≔}{\mathrm{cos}}{}\left({\mathrm{ω}}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t4}}{≔}{r}{*}{\mathrm{t2}}{+}{R}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t5}}{≔}{\mathrm{sin}}{}\left({\mathrm{φ}}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t6}}{≔}{\mathrm{sin}}{}\left({\mathrm{ω}}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}\left[{\mathrm{t1}}{*}{\mathrm{t4}}{,}{\mathrm{t5}}{*}{\mathrm{t4}}{,}{r}{*}{\mathrm{t6}}\right]\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (8)
 > $\mathrm{optimize}\left(\mathrm{HESSIAN}\left(\mathrm{torus},\mathrm{result_type}=\mathrm{list}\right)\right)$
 ${\mathbf{proc}}\left({\mathrm{φ}}{,}{\mathrm{ω}}{,}{R}{,}{r}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{local}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df}}{,}{\mathrm{df2}}{,}{\mathrm{dfr0}}{,}{\mathrm{dfr3}}{,}{\mathrm{dfr4}}{,}{\mathrm{t1}}{,}{\mathrm{t2}}{,}{\mathrm{t30}}{,}{\mathrm{t4}}{,}{\mathrm{t42}}{,}{\mathrm{t5}}{,}{\mathrm{t6}}{,}{\mathrm{t7}}{,}{\mathrm{t9}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t1}}{≔}{\mathrm{cos}}{}\left({\mathrm{φ}}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t2}}{≔}{\mathrm{cos}}{}\left({\mathrm{ω}}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t4}}{≔}{\mathrm{t2}}{*}{r}{+}{R}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t5}}{≔}{\mathrm{sin}}{}\left({\mathrm{φ}}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t6}}{≔}{\mathrm{sin}}{}\left({\mathrm{ω}}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df2}}{≔}{\mathrm{t1}}{*}{r}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df}}{≔}{\mathrm{array}}{}\left({1}{..}{8}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr0}}{≔}{\mathrm{array}}{}\left({1}{..}{8}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr3}}{≔}{\mathrm{array}}{}\left({1}{..}{8}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr4}}{≔}{\mathrm{array}}{}\left({1}{..}{8}\right){;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df}}{[}{5}{]}{≔}{−}{\mathrm{t4}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df}}{[}{4}{]}{≔}{−}{\mathrm{t5}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df}}{[}{3}{]}{≔}{\mathrm{df}}{[}{4}{]}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t7}}{≔}{\mathrm{df}}{[}{3}{]}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{df}}{[}{2}{]}{≔}{\mathrm{t7}}{*}{r}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr0}}{[}{7}{]}{≔}{−}{\mathrm{t6}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr0}}{[}{6}{]}{≔}{−}{\mathrm{df2}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr0}}{[}{1}{]}{≔}{\mathrm{dfr0}}{[}{7}{]}{*}{r}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr3}}{[}{2}{]}{≔}{\mathrm{df2}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr4}}{[}{8}{]}{≔}{\mathrm{dfr0}}{[}{7}{]}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr4}}{[}{6}{]}{≔}{−}{\mathrm{t5}}{*}{r}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t9}}{≔}{\mathrm{dfr4}}{[}{8}{]}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{dfr4}}{[}{5}{]}{≔}{\mathrm{t9}}{*}{r}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t30}}{≔}{\mathrm{t2}}{*}{\mathrm{t1}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathrm{t42}}{≔}\left[{0}{,}{0}{,}{0}{,}{0}\right]{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}\left[\left[\left[{\mathrm{df}}{[}{5}{]}{*}{\mathrm{t1}}{,}{−}{\mathrm{df}}{[}{2}{]}{*}{\mathrm{t6}}{,}{\mathrm{df}}{[}{3}{]}{,}{\mathrm{t7}}{*}{\mathrm{t2}}\right]{,}\left[{−}{\mathrm{dfr0}}{[}{1}{]}{*}{\mathrm{t5}}{,}{\mathrm{dfr0}}{[}{6}{]}{*}{\mathrm{t2}}{,}{0}{,}{\mathrm{dfr4}}{[}{8}{]}{*}{\mathrm{t1}}\right]{,}\left[{−}{\mathrm{t5}}{,}{0}{,}{0}{,}{0}\right]{,}\left[{−}{\mathrm{t5}}{*}{\mathrm{t2}}{,}{−}{\mathrm{t1}}{*}{\mathrm{t6}}{,}{0}{,}{0}\right]\right]{,}\left[\left[{−}{\mathrm{t5}}{*}{\mathrm{t4}}{,}{−}{\mathrm{dfr3}}{[}{2}{]}{*}{\mathrm{t6}}{,}{\mathrm{t1}}{,}{\mathrm{t30}}\right]{,}\left[{\mathrm{dfr4}}{[}{5}{]}{*}{\mathrm{t1}}{,}{\mathrm{dfr4}}{[}{6}{]}{*}{\mathrm{t2}}{,}{0}{,}{\mathrm{t9}}{*}{\mathrm{t5}}\right]{,}\left[{\mathrm{t1}}{,}{0}{,}{0}{,}{0}\right]{,}\left[{\mathrm{t30}}{,}{−}{\mathrm{t5}}{*}{\mathrm{t6}}{,}{0}{,}{0}\right]\right]{,}\left[{\mathrm{t42}}{,}\left[{0}{,}{−}{r}{*}{\mathrm{t6}}{,}{0}{,}{\mathrm{t2}}\right]{,}{\mathrm{t42}}{,}\left[{0}{,}{\mathrm{t2}}{,}{0}{,}{0}\right]\right]\right]\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (9)