Kedem Method for Derivatives of Implicit Functions

@(@\newcommand{\B}[1]{{\bf #1}} \newcommand{\R}[1]{{\rm #1}}@)@Kedem Method for Derivatives of Implicit Functions
Syntax

implicit_kedem kedem_ad(L_fun, F_fun, solve)

rk = kedem_ad.Forward(k, xk, store)

Purpose
Given a function @(@ L : \B{R}^{n \times m} \rightarrow \B{R}^m @)@, we define the implicit function @(@ Y : \B{R}^n \rightarrow \B{R}^m @)@ by @(@ L[ x , Y(x) ] = 0 @)@. The partial @(@ L_y [ x , Y(x) ] @)@ must be invertible for all @(@ x @)@ used by kedem_ad . We define the reduced function @(@ R(x) = F[ x , Y(x) ] @)@ where @(@ F : \B{R}^{n \times m} \rightarrow \B{R}^p @)@. The object kedem_ad can be used to compute derivatives of the reduced function @(@ R(x) @)@. Note that in the special case where @(@ F(x, y) = y @)@, @(@ R(x) = Y(x) @)@.

L_fun
This argument has prototype



     const CppAD::ADFun<double>& L_fun

and is the CppAD function object corresponding to @(@ L(x, y) @)@.

F_fun
This argument has prototype



     const CppAD::ADFun<double>& F_fun

and is the CppAD function object corresponding to @(@ F(x, y) @)@.

solve
This argument has prototype



     const Solve& solve

The type Solver must support the default constructor and the assignment operator. It must also support the following operations:

solve.function
In the syntax



     y = solve.function(x)

the argument x and the return value y have prototypes



     const VECTOR(double)&  x

     VECTOR(double)         y

where x.size() == n and y.size() == m . The return value satisfies the relation @(@ L(x, y) = 0 @)@.

solve.derivative
In the syntax



     b = solve.derivative(x, y)

the arguments have prototypes



     const VECTOR(double)&  x

     const VECTOR(double)&  y

The return value has prototype



     VECTOR(double) b

This returns the value of @(@ L_y (x, y) @)@ for subsequent calls to solve.linear . Only the elements of @(@ L_y (x, y) @)@ that depend on @(@ (x, y) @)@ need be included in the vector @(@ b @)@.

solve.linear
In the syntax



     u = solve.linear(b, v)

the arguments b , v and the return value u have prototypes



     const VECTOR(double)& b

     const VECTOR(double)& v

     const VECTOR(double)  u

where both vectors have size m . The argument b has prototype



     VECTOR(double)& b

The return value satisfies the equation @[@ u = L_y (x, y)^{-1} v @]@ where @(@ L_y (x, y) @)@ corresponds to b .

k
This is the order for this forward mode calculation and has prototype



     size_t k

Note that computing the k-th order, uses the results of the previous lower order calculations which are internally stored in the function objects.

xk
The argument xk has prototypes



     const VECTOR(double)& xk

its size is m and it is the k-th order Taylor coefficient for x . The results of the calculation are stored internally and used during higher order forward mode calculations. If k > 0 , there must have been a previous call to forward for order k - 1 .

store
This argument has prototype



     bool store

If it is true, the results of this forward mode calculation is store so that a future forward mode calculation of order k+1 is possible. This requires an extra pass over the operation sequence for L_fun . Hence store should be false when a higher order calculation will not be needed.

rk
The return value rk has prototype



     VECTOR(double) rk

its size is m and it is the k-th order Taylor coefficient for R(x) .

Method
Fix @(@ x : \B{R} \rightarrow \B{R}^n @)@ and define @(@ y : \B{R} \rightarrow \B{R}^m @)@ by @(@ x(t) = Y( x(t) ) @)@. The k-th order Taylor coefficients are @(@ x^k = x^{(k)} (0) / k ! @)@ and @(@ y^k = y^{(k)} (0) / k ! @)@.

Zero Order
The zero order coefficient @(@ y^0 @)@ is found using



     y0 = solve.function(x0)

Higher Orders
It follows from the definitions that @[@ 0 = L_x ( x^0 , y^0 ) x^{(1)} (t) + L_y ( x^0 , y^0 ) y^{(1)} (t) @]@ Taking more derivatives we see that, for @(@ k \geq 1 @)@, @[@ 0 = H_k ( x^0 , \cdots , x^k , y^0 , \cdots , y^{k-1} ) + L_y ( x^0 , y^0 ) y^k @]@ We use k-th order forward mode on L_fun , with @(@ y^k = 0 @)@, to determine @[@ z^k = H_k ( x^0 , \cdots , x^k , y^0 , \cdots , y^{k-1} ) @]@ We then solve the equation @[@ 0 = z^k + L_y ( x^0 , y^0 ) y^k @]@ to determine the proper value for @(@ y^k @)@.

Example
The routine test_circle_implicit_kedem is a simple example and test using this class. The routine test_control_reduced_objective is a control problem example and test using this class.

Input File: src/implicit_kedem.hpp