source: trunk/cppad/core/sparse_hessian.hpp @ 3904

Last change on this file since 3904 was 3904, checked in by bradbell, 2 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 4933d2c95c3d6736ab87a3cad45a8a43051c8d35
end hash code: 596554f576900948514ab4028728d583ada5333d

commit 596554f576900948514ab4028728d583ada5333d
Author: Brad Bell <bradbell@…>
Date: Sat Mar 25 05:59:25 2017 -0700

Document and test that it sparse_hes is more efficient when subset has
larger columns (not nrows).

commit bd3974dc485fd1c9b4d0f8e374661fa361da3dc0
Author: Brad Bell <bradbell@…>
Date: Sat Mar 25 04:49:09 2017 -0700

  1. Change sparse_hes to group colums instead of rows (uses reverse mode).
  2. Advance to cppad-20170325.

commit 22d86e19bdd9848e146165ce4f7279086abdaea1
Author: Brad Bell <bradbell@…>
Date: Sat Mar 25 02:53:58 2017 -0700

sparse_hes.hpp: Improve discussion of pattern subset.

File size: 23.7 KB
Line 
1# ifndef CPPAD_CORE_SPARSE_HESSIAN_HPP
2# define CPPAD_CORE_SPARSE_HESSIAN_HPP
3
4/* --------------------------------------------------------------------------
5CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
6
7CppAD is distributed under multiple licenses. This distribution is under
8the terms of the
9                    Eclipse Public License Version 1.0.
10
11A copy of this license is included in the COPYING file of this distribution.
12Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
13-------------------------------------------------------------------------- */
14
15/*
16$begin sparse_hessian$$
17$spell
18        jacobian
19        recomputed
20        CppAD
21        valarray
22        std
23        Bool
24        hes
25        const
26        Taylor
27        cppad
28        cmake
29        colpack
30$$
31
32$section Sparse Hessian$$
33
34$head Syntax$$
35$icode%hes% = %f%.SparseHessian(%x%, %w%)
36%hes% = %f%.SparseHessian(%x%, %w%, %p%)
37%n_sweep% = %f%.SparseHessian(%x%, %w%, %p%, %row%, %col%, %hes%, %work%)
38%$$
39
40$head Purpose$$
41We use $latex n$$ for the $cref/domain/seq_property/Domain/$$ size,
42and $latex m$$ for the $cref/range/seq_property/Range/$$ size of $icode f$$.
43We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the
44$cref/AD function/glossary/AD Function/$$
45corresponding to $icode f$$.
46The syntax above sets $icode hes$$ to the Hessian
47$latex \[
48        H(x) = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x)
49\] $$
50This routine takes advantage of the sparsity of the Hessian
51in order to reduce the amount of computation necessary.
52If $icode row$$ and $icode col$$ are present, it also takes
53advantage of the reduced set of elements of the Hessian that
54need to be computed.
55One can use speed tests (e.g. $cref speed_test$$)
56to verify that results are computed faster
57than when using the routine $cref Hessian$$.
58
59$head f$$
60The object $icode f$$ has prototype
61$codei%
62        ADFun<%Base%> %f%
63%$$
64Note that the $cref ADFun$$ object $icode f$$ is not $code const$$
65(see $cref/Uses Forward/sparse_hessian/Uses Forward/$$ below).
66
67$head x$$
68The argument $icode x$$ has prototype
69$codei%
70        const %VectorBase%& %x%
71%$$
72(see $cref/VectorBase/sparse_hessian/VectorBase/$$ below)
73and its size
74must be equal to $icode n$$, the dimension of the
75$cref/domain/seq_property/Domain/$$ space for $icode f$$.
76It specifies
77that point at which to evaluate the Hessian.
78
79$head w$$
80The argument $icode w$$ has prototype
81$codei%
82        const %VectorBase%& %w%
83%$$
84and size $latex m$$.
85It specifies the value of $latex w_i$$ in the expression
86for $icode hes$$.
87The more components of $latex w$$ that are identically zero,
88the more sparse the resulting Hessian may be (and hence the more efficient
89the calculation of $icode hes$$ may be).
90
91$head p$$
92The argument $icode p$$ is optional and has prototype
93$codei%
94        const %VectorSet%& %p%
95%$$
96(see $cref/VectorSet/sparse_hessian/VectorSet/$$ below)
97If it has elements of type $code bool$$,
98its size is $latex n * n$$.
99If it has elements of type $code std::set<size_t>$$,
100its size is $latex n$$ and all its set elements are between
101zero and $latex n - 1$$.
102It specifies a
103$cref/sparsity pattern/glossary/Sparsity Pattern/$$
104for the Hessian $latex H(x)$$.
105
106$subhead Purpose$$
107If this sparsity pattern does not change between calls to
108$codei SparseHessian$$, it should be faster to calculate $icode p$$ once and
109pass this argument to $codei SparseHessian$$.
110If you specify $icode p$$, CppAD will use the same
111type of sparsity representation
112(vectors of $code bool$$ or vectors of $code std::set<size_t>$$)
113for its internal calculations.
114Otherwise, the representation
115for the internal calculations is unspecified.
116
117$subhead work$$
118If you specify $icode work$$ in the calling sequence,
119it is not necessary to keep the sparsity pattern; see the heading
120$cref/p/sparse_hessian/work/p/$$ under the $icode work$$ description.
121
122$subhead Column Subset$$
123If the arguments $icode row$$ and $icode col$$ are present,
124and $cref/color_method/sparse_hessian/work/color_method/$$ is
125$code cppad.general$$ or $code cppad.symmetric$$,
126it is not necessary to compute the entire sparsity pattern.
127Only the following subset of column values will matter:
128$codei%
129        { %col%[%k%] : %k% = 0 , %...% , %K%-1 }
130%$$.
131
132
133$head row, col$$
134The arguments $icode row$$ and $icode col$$ are optional and have prototype
135$codei%
136        const %VectorSize%& %row%
137        const %VectorSize%& %col%
138%$$
139(see $cref/VectorSize/sparse_hessian/VectorSize/$$ below).
140They specify which rows and columns of $latex H (x)$$ are
141returned and in what order.
142We use $latex K$$ to denote the value $icode%hes%.size()%$$
143which must also equal the size of $icode row$$ and $icode col$$.
144Furthermore,
145for $latex k = 0 , \ldots , K-1$$, it must hold that
146$latex row[k] < n$$ and $latex col[k] < n$$.
147In addition,
148all of the $latex (row[k], col[k])$$ pairs must correspond to a true value
149in the sparsity pattern $icode p$$.
150
151$head hes$$
152The result $icode hes$$ has prototype
153$codei%
154        %VectorBase% %hes%
155%$$
156In the case where $icode row$$ and $icode col$$ are not present,
157the size of $icode hes$$ is $latex n * n$$ and
158its size is $latex n * n$$.
159In this case, for $latex i = 0 , \ldots , n - 1 $$
160and $latex ell = 0 , \ldots , n - 1$$
161$latex \[
162        hes [ j * n + \ell ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } ( x )
163\] $$
164$pre
165
166$$
167In the case where the arguments $icode row$$ and $icode col$$ are present,
168we use $latex K$$ to denote the size of $icode hes$$.
169The input value of its elements does not matter.
170Upon return, for $latex k = 0 , \ldots , K - 1$$,
171$latex \[
172        hes [ k ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } (x)
173        \; , \;
174        \; {\rm where} \;
175        j = row[k]
176        \; {\rm and } \;
177        \ell = col[k]
178\] $$
179
180$head work$$
181If this argument is present, it has prototype
182$codei%
183        sparse_hessian_work& %work%
184%$$
185This object can only be used with the routines $code SparseHessian$$.
186During its the first use, information is stored in $icode work$$.
187This is used to reduce the work done by future calls to $code SparseHessian$$
188with the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$.
189If a future call is made where any of these values have changed,
190you must first call $icode%work%.clear()%$$
191to inform CppAD that this information needs to be recomputed.
192
193$subhead color_method$$
194The coloring algorithm determines which rows and columns
195can be computed during the same sweep.
196This field has prototype
197$codei%
198        std::string %work%.color_method
199%$$
200This value only matters on the first call to $code sparse_hessian$$ that
201follows the $icode work$$ constructor or a call to
202$icode%work%.clear()%$$.
203$codei%
204
205"cppad.symmetric"
206%$$
207This is the default coloring method (after a constructor or $code clear()$$).
208It takes advantage of the fact that the Hessian matrix
209is symmetric to find a coloring that requires fewer
210$cref/sweeps/sparse_hessian/n_sweep/$$.
211$codei%
212
213"cppad.general"
214%$$
215This is the same as the $code "cppad"$$ method for the
216$cref/sparse_jacobian/sparse_jacobian/work/color_method/$$ calculation.
217$codei%
218
219"colpack.star"
220%$$
221This method requires that
222$cref colpack_prefix$$ was specified on the
223$cref/cmake command/cmake/CMake Command/$$ line.
224It also takes advantage of the fact that the Hessian matrix is symmetric.
225
226$subhead p$$
227If $icode work$$ is present, and it is not the first call after
228its construction or a clear,
229the sparsity pattern $icode p$$ is not used.
230This enables one to free the sparsity pattern
231and still compute corresponding sparse Hessians.
232
233$head n_sweep$$
234The return value $icode n_sweep$$ has prototype
235$codei%
236        size_t %n_sweep%
237%$$
238It is the number of first order forward sweeps
239used to compute the requested Hessian values.
240Each first forward sweep is followed by a second order reverse sweep
241so it is also the number of reverse sweeps.
242This is proportional to the total work that $code SparseHessian$$ does,
243not counting the zero order forward sweep,
244or the work to combine multiple columns into a single
245forward-reverse sweep pair.
246
247$head VectorBase$$
248The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with
249$cref/elements of type/SimpleVector/Elements of Specified Type/$$
250$icode Base$$.
251The routine $cref CheckSimpleVector$$ will generate an error message
252if this is not the case.
253
254$head VectorSet$$
255The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with
256$cref/elements of type/SimpleVector/Elements of Specified Type/$$
257$code bool$$ or $code std::set<size_t>$$;
258see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
259of the difference.
260The routine $cref CheckSimpleVector$$ will generate an error message
261if this is not the case.
262
263$subhead Restrictions$$
264If $icode VectorSet$$ has elements of $code std::set<size_t>$$,
265then $icode%p%[%i%]%$$ must return a reference (not a copy) to the
266corresponding set.
267According to section 26.3.2.3 of the 1998 C++ standard,
268$code std::valarray< std::set<size_t> >$$ does not satisfy
269this condition.
270
271$head VectorSize$$
272The type $icode VectorSize$$ must be a $cref SimpleVector$$ class with
273$cref/elements of type/SimpleVector/Elements of Specified Type/$$
274$code size_t$$.
275The routine $cref CheckSimpleVector$$ will generate an error message
276if this is not the case.
277
278$head Uses Forward$$
279After each call to $cref Forward$$,
280the object $icode f$$ contains the corresponding
281$cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
282After a call to any of the sparse Hessian routines,
283the zero order Taylor coefficients correspond to
284$icode%f%.Forward(0, %x%)%$$
285and the other coefficients are unspecified.
286
287$children%
288        example/sparse/sparse_hessian.cpp%
289        example/sparse/sub_sparse_hes.cpp%
290        example/sparse/sparse_sub_hes.cpp
291%$$
292
293$head Example$$
294The routine
295$cref sparse_hessian.cpp$$
296is examples and tests of $code sparse_hessian$$.
297It return $code true$$, if it succeeds and $code false$$ otherwise.
298
299$head Subset Hessian$$
300The routine
301$cref sub_sparse_hes.cpp$$
302is an example and test that compute a sparse Hessian
303for a subset of the variables.
304It returns $code true$$, for success, and $code false$$ otherwise.
305
306$end
307-----------------------------------------------------------------------------
308*/
309# include <cppad/local/std_set.hpp>
310# include <cppad/local/color_general.hpp>
311# include <cppad/local/color_symmetric.hpp>
312
313namespace CppAD { // BEGIN_CPPAD_NAMESPACE
314/*!
315\file sparse_hessian.hpp
316Sparse Hessian driver routine and helper functions.
317*/
318// ===========================================================================
319/*!
320class used by SparseHessian to hold information
321so it does not need to be recomputed.
322*/
323class sparse_hessian_work {
324        public:
325                /// Coloring method: "cppad", or "colpack"
326                /// (this field is set by user)
327                std::string color_method;
328                /// row and column indicies for return values
329                /// (some may be reflected by star coloring algorithm)
330                CppAD::vector<size_t> row;
331                CppAD::vector<size_t> col;
332                /// indices that sort the user row and col arrays by color
333                CppAD::vector<size_t> order;
334                /// results of the coloring algorithm
335                CppAD::vector<size_t> color;
336
337                /// constructor
338                sparse_hessian_work(void) : color_method("cppad.symmetric")
339                { }
340                /// inform CppAD that this information needs to be recomputed
341                void clear(void)
342                {       color_method = "cppad.symmetric";
343                        row.clear();
344                        col.clear();
345                        order.clear();
346                        color.clear();
347                }
348};
349// ===========================================================================
350/*!
351Private helper function that does computation for all Sparse Hessian cases.
352
353\tparam Base
354is the base type for the recording that is stored in this ADFun<Base object.
355
356\tparam VectorBase
357is a simple vector class with elements of type \a Base.
358
359\tparam VectorSet
360is a simple vector class with elements of type
361\c bool or \c std::set<size_t>.
362
363\tparam VectorSize
364is sparse_pack or sparse_list.
365
366\param x [in]
367is a vector specifing the point at which to compute the Hessian.
368
369\param w [in]
370is the weighting vector that defines a scalar valued function by
371a weighted sum of the components of the vector valued function
372$latex F(x)$$.
373
374\param sparsity [in]
375is the sparsity pattern for the Hessian that we are calculating.
376
377\param user_row [in]
378is the vector of row indices for the returned Hessian values.
379
380\param user_col [in]
381is the vector of columns indices for the returned Hessian values.
382It must have the same size as user_row.
383
384\param hes [out]
385is the vector of Hessian values.
386It must have the same size as user_row.
387The return value <code>hes[k]</code> is the second partial of
388\f$ w^{\rm T} F(x)\f$ with respect to the
389<code>row[k]</code> and <code>col[k]</code> component of \f$ x\f$.
390
391\param work
392This structure contains information that is computed by \c SparseHessianCompute.
393If the sparsity pattern, \c row vector, or \c col vectors
394are not the same between calls to \c SparseHessianCompute,
395\c work.clear() must be called to reinitialize \c work.
396
397\return
398Is the number of first order forward sweeps used to compute the
399requested Hessian values.
400(This is also equal to the number of second order reverse sweeps.)
401The total work, not counting the zero order
402forward sweep, or the time to combine computations, is proportional to this
403return value.
404*/
405template<class Base>
406template <class VectorBase, class VectorSet, class VectorSize>
407size_t ADFun<Base>::SparseHessianCompute(
408        const VectorBase&           x           ,
409        const VectorBase&           w           ,
410              VectorSet&            sparsity    ,
411        const VectorSize&           user_row    ,
412        const VectorSize&           user_col    ,
413              VectorBase&           hes         ,
414              sparse_hessian_work&  work        )
415{
416        using   CppAD::vectorBool;
417        size_t i, k, ell;
418
419        CppAD::vector<size_t>& row(work.row);
420        CppAD::vector<size_t>& col(work.col);
421        CppAD::vector<size_t>& color(work.color);
422        CppAD::vector<size_t>& order(work.order);
423
424        size_t n = Domain();
425
426        // some values
427        const Base zero(0);
428        const Base one(1);
429
430        // check VectorBase is Simple Vector class with Base type elements
431        CheckSimpleVector<Base, VectorBase>();
432
433        // number of components of Hessian that are required
434        size_t K = hes.size();
435        CPPAD_ASSERT_UNKNOWN( size_t( user_row.size() ) == K );
436        CPPAD_ASSERT_UNKNOWN( size_t( user_col.size() ) == K );
437
438        CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n );
439        CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n );
440        CPPAD_ASSERT_UNKNOWN( row.size() == 0   || row.size() == K );
441        CPPAD_ASSERT_UNKNOWN( col.size() == 0   || col.size() == K );
442
443
444        // Point at which we are evaluating the Hessian
445        Forward(0, x);
446
447        // check for case where nothing (except Forward above) to do
448        if( K == 0 )
449                return 0;
450
451        // Rows of the Hessian (i below) correspond to the forward mode index
452        // and columns (j below) correspond to the reverse mode index.
453        if( color.size() == 0 )
454        {
455                CPPAD_ASSERT_UNKNOWN( sparsity.n_set() ==  n );
456                CPPAD_ASSERT_UNKNOWN( sparsity.end() ==  n );
457
458                // copy user rwo and col to work space
459                row.resize(K);
460                col.resize(K);
461                for(k = 0; k < K; k++)
462                {       row[k] = user_row[k];
463                        col[k] = user_col[k];
464                }
465
466                // execute coloring algorithm
467                color.resize(n);
468                if( work.color_method == "cppad.general" )
469                        local::color_general_cppad(sparsity, row, col, color);
470                else if( work.color_method == "cppad.symmetric" )
471                        local::color_symmetric_cppad(sparsity, row, col, color);
472                else if( work.color_method == "colpack.star" )
473                {
474# if CPPAD_HAS_COLPACK
475                        local::color_symmetric_colpack(sparsity, row, col, color);
476# else
477                        CPPAD_ASSERT_KNOWN(
478                                false,
479                                "SparseHessian: work.color_method = colpack.star "
480                                "and colpack_prefix missing from cmake command line."
481                        );
482# endif
483                }
484                else
485                {       CPPAD_ASSERT_KNOWN(
486                                false,
487                                "SparseHessian: work.color_method is not valid."
488                        );
489                }
490
491                // put sorting indices in color order
492                VectorSize key(K);
493                order.resize(K);
494                for(k = 0; k < K; k++)
495                        key[k] = color[ row[k] ];
496                index_sort(key, order);
497
498        }
499        size_t n_color = 1;
500        for(ell = 0; ell < n; ell++) if( color[ell] < n )
501                n_color = std::max(n_color, color[ell] + 1);
502
503        // direction vector for calls to forward (rows of the Hessian)
504        VectorBase u(n);
505
506        // location for return values from reverse (columns of the Hessian)
507        VectorBase ddw(2 * n);
508
509        // initialize the return value
510        for(k = 0; k < K; k++)
511                hes[k] = zero;
512
513        // loop over colors
514        k = 0;
515        for(ell = 0; ell < n_color; ell++)
516        {       CPPAD_ASSERT_UNKNOWN( color[ row[ order[k] ] ] == ell );
517
518                // combine all rows with this color
519                for(i = 0; i < n; i++)
520                {       u[i] = zero;
521                        if( color[i] == ell )
522                                u[i] = one;
523                }
524                // call forward mode for all these rows at once
525                Forward(1, u);
526
527                // evaluate derivative of w^T * F'(x) * u
528                ddw = Reverse(2, w);
529
530                // set the corresponding components of the result
531                while( k < K && color[ row[ order[k] ] ] == ell )
532                {       hes[ order[k] ] = ddw[ col[ order[k] ] * 2 + 1 ];
533                        k++;
534                }
535        }
536        return n_color;
537}
538// ===========================================================================
539// Public Member Functions
540// ===========================================================================
541/*!
542Compute user specified subset of a sparse Hessian.
543
544The C++ source code corresponding to this operation is
545\verbatim
546        SparceHessian(x, w, p, row, col, hes, work)
547\endverbatim
548
549\tparam Base
550is the base type for the recording that is stored in this ADFun<Base object.
551
552\tparam VectorBase
553is a simple vector class with elements of type \a Base.
554
555\tparam VectorSet
556is a simple vector class with elements of type
557\c bool or \c std::set<size_t>.
558
559\tparam VectorSize
560is a simple vector class with elements of type \c size_t.
561
562\param x [in]
563is a vector specifing the point at which to compute the Hessian.
564
565\param w [in]
566is the weighting vector that defines a scalar valued function by
567a weighted sum of the components of the vector valued function
568$latex F(x)$$.
569
570\param p [in]
571is the sparsity pattern for the Hessian that we are calculating.
572
573\param row [in]
574is the vector of row indices for the returned Hessian values.
575
576\param col [in]
577is the vector of columns indices for the returned Hessian values.
578It must have the same size are r.
579
580\param hes [out]
581is the vector of Hessian values.
582It must have the same size are r.
583The return value <code>hes[k]</code> is the second partial of
584\f$ w^{\rm T} F(x)\f$ with respect to the
585<code>row[k]</code> and <code>col[k]</code> component of \f$ x\f$.
586
587\param work
588This structure contains information that is computed by \c SparseHessianCompute.
589If the sparsity pattern, \c row vector, or \c col vectors
590are not the same between calls to \c SparseHessian,
591\c work.clear() must be called to reinitialize \c work.
592
593\return
594Is the number of first order forward sweeps used to compute the
595requested Hessian values.
596(This is also equal to the number of second order reverse sweeps.)
597The total work, not counting the zero order
598forward sweep, or the time to combine computations, is proportional to this
599return value.
600*/
601template<class Base>
602template <class VectorBase, class VectorSet, class VectorSize>
603size_t ADFun<Base>::SparseHessian(
604        const VectorBase&     x    ,
605        const VectorBase&     w    ,
606        const VectorSet&      p    ,
607        const VectorSize&     row  ,
608        const VectorSize&     col  ,
609        VectorBase&           hes  ,
610        sparse_hessian_work&  work )
611{
612        size_t n    = Domain();
613        size_t K    = hes.size();
614# ifndef NDEBUG
615        size_t k;
616        CPPAD_ASSERT_KNOWN(
617                size_t(x.size()) == n ,
618                "SparseHessian: size of x not equal domain dimension for f."
619        );
620        CPPAD_ASSERT_KNOWN(
621                size_t(row.size()) == K && size_t(col.size()) == K ,
622                "SparseHessian: either r or c does not have the same size as ehs."
623        );
624        CPPAD_ASSERT_KNOWN(
625                work.color.size() == 0 || work.color.size() == n,
626                "SparseHessian: invalid value in work."
627        );
628        for(k = 0; k < K; k++)
629        {       CPPAD_ASSERT_KNOWN(
630                        row[k] < n,
631                        "SparseHessian: invalid value in r."
632                );
633                CPPAD_ASSERT_KNOWN(
634                        col[k] < n,
635                        "SparseHessian: invalid value in c."
636                );
637        }
638        if( work.color.size() != 0 )
639                for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN(
640                        work.color[j] <= n,
641                        "SparseHessian: invalid value in work."
642        );
643# endif
644        // check for case where there is nothing to compute
645        size_t n_sweep = 0;
646        if( K == 0 )
647                return n_sweep;
648
649        typedef typename VectorSet::value_type Set_type;
650        typedef typename local::internal_sparsity<Set_type>::pattern_type Pattern_type;
651        Pattern_type s;
652        if( work.color.size() == 0 )
653        {       bool transpose = false;
654                const char* error_msg = "SparseHessian: sparsity pattern"
655                " does not have proper row or column dimension";
656                sparsity_user2internal(s, p, n, n, transpose, error_msg);
657        }
658        n_sweep = SparseHessianCompute(x, w, s, row, col, hes, work);
659        return n_sweep;
660}
661/*!
662Compute a sparse Hessian.
663
664The C++ source code coresponding to this operation is
665\verbatim
666        hes = SparseHessian(x, w, p)
667\endverbatim
668
669
670\tparam Base
671is the base type for the recording that is stored in this
672ADFun<Base object.
673
674\tparam VectorBase
675is a simple vector class with elements of the \a Base.
676
677\tparam VectorSet
678is a simple vector class with elements of type
679\c bool or \c std::set<size_t>.
680
681\param x [in]
682is a vector specifing the point at which to compute the Hessian.
683
684\param w [in]
685The Hessian is computed for a weighted sum of the components
686of the function corresponding to this ADFun<Base> object.
687The argument \a w specifies the weights for each component.
688It must have size equal to the range dimension for this ADFun<Base> object.
689
690\param p [in]
691is a sparsity pattern for the Hessian.
692
693\return
694Will be a vector of size \c n * n containing the Hessian of
695at the point specified by \a x
696(where \c n is the domain dimension for this ADFun<Base> object).
697*/
698template <class Base>
699template <class VectorBase, class VectorSet>
700VectorBase ADFun<Base>::SparseHessian(
701        const VectorBase& x, const VectorBase& w, const VectorSet& p
702)
703{       size_t i, j, k;
704
705        size_t n = Domain();
706        VectorBase hes(n * n);
707
708        CPPAD_ASSERT_KNOWN(
709                size_t(x.size()) == n,
710                "SparseHessian: size of x not equal domain size for f."
711        );
712
713        typedef typename VectorSet::value_type Set_type;
714        typedef typename local::internal_sparsity<Set_type>::pattern_type Pattern_type;
715
716        // initialize the return value as zero
717        Base zero(0);
718        for(i = 0; i < n; i++)
719                for(j = 0; j < n; j++)
720                        hes[i * n + j] = zero;
721
722        // arguments to SparseHessianCompute
723        Pattern_type          s;
724        CppAD::vector<size_t> row;
725        CppAD::vector<size_t> col;
726        sparse_hessian_work   work;
727        bool transpose = false;
728        const char* error_msg = "SparseHessian: sparsity pattern"
729        " does not have proper row or column dimension";
730        sparsity_user2internal(s, p, n, n, transpose, error_msg);
731        k = 0;
732        for(i = 0; i < n; i++)
733        {       typename Pattern_type::const_iterator itr(s, i);
734                j = *itr;
735                while( j != s.end() )
736                {       row.push_back(i);
737                        col.push_back(j);
738                        k++;
739                        j = *(++itr);
740                }
741        }
742        size_t K = k;
743        VectorBase H(K);
744
745        // now we have folded this into the following case
746        SparseHessianCompute(x, w, s, row, col, H, work);
747
748        // now set the non-zero return values
749        for(k = 0; k < K; k++)
750                hes[ row[k] * n + col[k] ] = H[k];
751
752        return hes;
753}
754/*!
755Compute a sparse Hessian
756
757The C++ source code coresponding to this operation is
758\verbatim
759        hes = SparseHessian(x, w)
760\endverbatim
761
762
763\tparam Base
764is the base type for the recording that is stored in this
765ADFun<Base object.
766
767\tparam VectorBase
768is a simple vector class with elements of the \a Base.
769
770\param x [in]
771is a vector specifing the point at which to compute the Hessian.
772
773\param w [in]
774The Hessian is computed for a weighted sum of the components
775of the function corresponding to this ADFun<Base> object.
776The argument \a w specifies the weights for each component.
777It must have size equal to the range dimension for this ADFun<Base> object.
778
779\return
780Will be a vector of size \c n * n containing the Hessian of
781at the point specified by \a x
782(where \c n is the domain dimension for this ADFun<Base> object).
783*/
784template <class Base>
785template <class VectorBase>
786VectorBase ADFun<Base>::SparseHessian(const VectorBase &x, const VectorBase &w)
787{       size_t i, j, k;
788        typedef CppAD::vectorBool VectorBool;
789
790        size_t m = Range();
791        size_t n = Domain();
792
793        // determine the sparsity pattern p for Hessian of w^T F
794        VectorBool r(n * n);
795        for(j = 0; j < n; j++)
796        {       for(k = 0; k < n; k++)
797                        r[j * n + k] = false;
798                r[j * n + j] = true;
799        }
800        ForSparseJac(n, r);
801        //
802        VectorBool s(m);
803        for(i = 0; i < m; i++)
804                s[i] = w[i] != 0;
805        VectorBool p = RevSparseHes(n, s);
806
807        // compute sparse Hessian
808        return SparseHessian(x, w, p);
809}
810
811} // END_CPPAD_NAMESPACE
812# endif
Note: See TracBrowser for help on using the repository browser.