source: trunk/cppad/core/sparse_hes.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: 13.9 KB
Line 
1# ifndef CPPAD_CORE_SPARSE_HES_HPP
2# define CPPAD_CORE_SPARSE_HES_HPP
3/* --------------------------------------------------------------------------
4CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
5
6CppAD is distributed under multiple licenses. This distribution is under
7the terms of the
8                    Eclipse Public License Version 1.0.
9
10A copy of this license is included in the COPYING file of this distribution.
11Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
12-------------------------------------------------------------------------- */
13
14/*
15$begin sparse_hes$$
16$spell
17        const
18        Taylor
19        rc
20        rcv
21        nr
22        nc
23        hes
24        std
25        cppad
26        colpack
27        cmake
28        Jacobian
29$$
30
31$section Computing Sparse Hessians$$
32
33$head Syntax$$
34$icode%n_sweep% = %f%.sparse_hes(
35        %x%, %w%, %subset%, %pattern%, %coloring%, %work%
36)%$$
37
38$head Purpose$$
39We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ to denote the
40function corresponding to $icode f$$.
41Here $icode n$$ is the $cref/domain/seq_property/Domain/$$ size,
42and $icode m$$ is the $cref/range/seq_property/Range/$$ size, or $icode f$$.
43The syntax above takes advantage of sparsity when computing the Hessian
44$latex \[
45        H(x) = \dpow{2}{x} \sum_{i=0}^{m-1} w_i F_i (x)
46\] $$
47In the sparse case, this should be faster and take less memory than
48$cref Hessian$$.
49The matrix element $latex H_{i,j} (x)$$ is the second partial of
50$latex w^\R{T} F (x)$$ with respect to $latex x_i$$ and $latex x_j$$.
51
52$head SizeVector$$
53The type $icode SizeVector$$ is a $cref SimpleVector$$ class with
54$cref/elements of type/SimpleVector/Elements of Specified Type/$$
55$code size_t$$.
56
57$head BaseVector$$
58The type $icode BaseVector$$ is a $cref SimpleVector$$ class with
59$cref/elements of type/SimpleVector/Elements of Specified Type/$$
60$code size_t$$.
61
62$head f$$
63This object has prototype
64$codei%
65        ADFun<%Base%> %f%
66%$$
67Note that the Taylor coefficients stored in $icode f$$ are affected
68by this operation; see
69$cref/uses forward/sparse_hes/Uses Forward/$$ below.
70
71$head x$$
72This argument has prototype
73$codei%
74        const %BaseVector%& %x%
75%$$
76and its size is $icode n$$.
77It specifies the point at which to evaluate the Hessian
78$latex H(x)$$.
79
80$head w$$
81This argument has prototype
82$codei%
83        const %BaseVector%& %w%
84%$$
85and its size is $icode m$$.
86It specifies the weight for each of the components of $latex F(x)$$;
87i.e. $latex w_i$$ is the weight for $latex F_i (x)$$.
88
89$head subset$$
90This argument has prototype
91$codei%
92        sparse_rcv<%SizeVector%, %BaseVector%>& %subset%
93%$$
94Its row size and column size is $icode n$$; i.e.,
95$icode%subset%.nr() == %n%$$ and $icode%subset%.nc() == %n%$$.
96It specifies which elements of the Hessian are computed.
97$list number$$
98The input value of its value vector
99$icode%subset%.val()%$$ does not matter.
100Upon return it contains the value of the corresponding elements
101of the Hessian.
102$lnext
103All of the row, column pairs in $icode subset$$ must also appear in
104$icode pattern$$; i.e., they must be possibly non-zero.
105$lnext
106The Hessian is symmetric, so one has a choice as to which off diagonal
107elements to put in $icode subset$$.
108It will probably be more efficient if one makes this choice so that
109the there are more entries in each non-zero column of $icode subset$$;
110see $cref/n_sweep/sparse_hes/n_sweep/$$ below.
111$lend
112
113$head pattern$$
114This argument has prototype
115$codei%
116        const sparse_rc<%SizeVector%>& %pattern%
117%$$
118Its row size and column size is $icode n$$; i.e.,
119$icode%pattern%.nr() == %n%$$ and $icode%pattern%.nc() == %n%$$.
120It is a sparsity pattern for the Hessian $latex H(x)$$.
121If the $th i$$ row ($th j$$ column) does not appear in $icode subset$$,
122the $th i$$ row ($th j$$ column) of $icode pattern$$ does not matter
123and need not be computed.
124This argument is not used (and need not satisfy any conditions),
125when $cref/work/sparse_hes/work/$$ is non-empty.
126
127$subhead subset$$
128If the $th i$$ row and $th i$$ column do not appear in $icode subset$$,
129the $th i$$ row and column of $icode pattern$$ do not matter.
130In this case the $th i-th$$ row and column may have no entries in
131$icode pattern$$ even though they are possibly non-zero in $latex H(x)$$.
132(This can be used to reduce the amount of computation required to find
133$icode pattern$$.)
134
135$head coloring$$
136The coloring algorithm determines which rows and columns
137can be computed during the same sweep.
138This field has prototype
139$codei%
140        const std::string& %coloring%
141%$$
142This value only matters when work is empty; i.e.,
143after the $icode work$$ constructor or $icode%work%.clear()%$$.
144
145$subhead cppad.symmetric$$
146This coloring takes advantage of the fact that the Hessian matrix
147is symmetric when find a coloring that requires fewer
148$cref/sweeps/sparse_hes/n_sweep/$$.
149
150$subhead cppad.general$$
151This is the same as the sparse Jacobian
152$cref/cppad/sparse_jac/coloring/cppad/$$ method
153which does not take advantage of symmetry.
154
155$subhead colpack.star$$
156If $cref colpack_prefix$$ was specified on the
157$cref/cmake command/cmake/CMake Command/$$ line,
158you can set $icode coloring$$ to $code colpack.star$$.
159This also takes advantage of the fact that the Hessian matrix is symmetric.
160
161$head work$$
162This argument has prototype
163$codei%
164        sparse_hes_work& %work%
165%$$
166We refer to its initial value,
167and its value after $icode%work%.clear()%$$, as empty.
168If it is empty, information is stored in $icode work$$.
169This can be used to reduce computation when
170a future call is for the same object $icode f$$,
171and the same subset of the Hessian.
172If either of these values change, use $icode%work%.clear()%$$ to
173empty this structure.
174
175$head n_sweep$$
176The return value $icode n_sweep$$ has prototype
177$codei%
178        size_t %n_sweep%
179%$$
180It is the number of first order forward sweeps
181used to compute the requested Hessian values.
182Each first forward sweep is followed by a second order reverse sweep
183so it is also the number of reverse sweeps.
184It is also the number of colors determined by the coloring method
185mentioned above.
186This is proportional to the total computational work,
187not counting the zero order forward sweep,
188or combining multiple columns and rows into a single sweep.
189
190$head Uses Forward$$
191After each call to $cref Forward$$,
192the object $icode f$$ contains the corresponding
193$cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
194After a call to $code sparse_hes$$
195the zero order coefficients correspond to
196$codei%
197        %f%.Forward(0, %x%)
198%$$
199All the other forward mode coefficients are unspecified.
200
201$head Example$$
202$children%
203        example/sparse/sparse_hes.cpp
204%$$
205The files $cref sparse_hes.cpp$$
206is an example and test of $code sparse_hes$$.
207It returns $code true$$, if it succeeds, and $code false$$ otherwise.
208
209$head Subset Hessian$$
210The routine
211$cref sparse_sub_hes.cpp$$
212is an example and test that compute a subset of a sparse Hessian.
213It returns $code true$$, for success, and $code false$$ otherwise.
214
215$end
216*/
217# include <cppad/core/cppad_assert.hpp>
218# include <cppad/local/sparse_internal.hpp>
219# include <cppad/local/color_general.hpp>
220# include <cppad/local/color_symmetric.hpp>
221
222/*!
223\file sparse_hes.hpp
224Sparse Hessian calculation routines.
225*/
226namespace CppAD { // BEGIN_CPPAD_NAMESPACE
227
228/*!
229Class used to hold information used by Sparse Hessian routine in this file,
230so it does not need to be recomputed every time.
231*/
232class sparse_hes_work {
233        public:
234                /// row and column indicies for return values
235                /// (some may be reflected by symmetric coloring algorithms)
236                CppAD::vector<size_t> row;
237                CppAD::vector<size_t> col;
238                /// indices that sort the row and col arrays by color
239                CppAD::vector<size_t> order;
240                /// results of the coloring algorithm
241                CppAD::vector<size_t> color;
242
243                /// constructor
244                sparse_hes_work(void)
245                { }
246                /// inform CppAD that this information needs to be recomputed
247                void clear(void)
248                {
249                        row.clear();
250                        col.clear();
251                        order.clear();
252                        color.clear();
253                }
254};
255// ----------------------------------------------------------------------------
256/*!
257Calculate sparse Hessians using forward mode
258
259\tparam Base
260the base type for the recording that is stored in the ADFun object.
261
262\tparam SizeVector
263a simple vector class with elements of type size_t.
264
265\tparam BaseVector
266a simple vector class with elements of type Base.
267
268\param x
269a vector of length n, the number of independent variables in f
270(this ADFun object).
271
272\param w
273a vector of length m, the number of dependent variables in f
274(this ADFun object).
275
276\param subset
277specifices the subset of the sparsity pattern where the Hessian is evaluated.
278subset.nr() == n,
279subset.nc() == n.
280
281\param pattern
282is a sparsity pattern for the Hessian of w^T * f;
283pattern.nr() == n,
284pattern.nc() == n,
285where m is number of dependent variables in f.
286
287\param coloring
288determines which coloring algorithm is used.
289This must be cppad.symmetric, cppad.general, or colpack.star.
290
291\param work
292this structure must be empty, or contain the information stored
293by a previous call to sparse_hes.
294The previous call must be for the same ADFun object f
295and the same subset.
296
297\return
298This is the number of first order forward
299(and second order reverse) sweeps used to compute thhe Hessian.
300*/
301template <class Base>
302template <class SizeVector, class BaseVector>
303size_t ADFun<Base>::sparse_hes(
304        const BaseVector&                    x        ,
305        const BaseVector&                    w        ,
306        sparse_rcv<SizeVector , BaseVector>& subset   ,
307        const sparse_rc<SizeVector>&         pattern  ,
308        const std::string&                   coloring ,
309        sparse_hes_work&                     work     )
310{       size_t n = Domain();
311        //
312        CPPAD_ASSERT_KNOWN(
313                subset.nr() == n,
314                "sparse_hes: subset.nr() not equal domain dimension for f"
315        );
316        CPPAD_ASSERT_KNOWN(
317                subset.nc() == n,
318                "sparse_hes: subset.nc() not equal domain dimension for f"
319        );
320        CPPAD_ASSERT_KNOWN(
321                size_t( x.size() ) == n,
322                "sparse_hes: x.size() not equal domain dimension for f"
323        );
324        CPPAD_ASSERT_KNOWN(
325                size_t( w.size() ) == Range(),
326                "sparse_hes: w.size() not equal range dimension for f"
327        );
328        //
329        // work information
330        vector<size_t>& row(work.row);
331        vector<size_t>& col(work.col);
332        vector<size_t>& color(work.color);
333        vector<size_t>& order(work.order);
334        //
335        // subset information
336        const SizeVector& subset_row( subset.row() );
337        const SizeVector& subset_col( subset.col() );
338        //
339        // point at which we are evaluationg the Hessian
340        Forward(0, x);
341        //
342        // number of elements in the subset
343        size_t K = subset.nnz();
344        //
345        // check for case were there is nothing to do
346        // (except for call to Forward(0, x)
347        if( K == 0 )
348                return 0;
349        //
350# ifndef NDEBUG
351        if( color.size() != 0 )
352        {       CPPAD_ASSERT_KNOWN(
353                        color.size() == n,
354                        "sparse_hes: work is non-empty and conditions have changed"
355                );
356                CPPAD_ASSERT_KNOWN(
357                        row.size() == K,
358                        "sparse_hes: work is non-empty and conditions have changed"
359                );
360                CPPAD_ASSERT_KNOWN(
361                        col.size() == K,
362                        "sparse_hes: work is non-empty and conditions have changed"
363                );
364                //
365                for(size_t k = 0; k < K; k++)
366                {       bool ok = row[k] == subset_row[k] && col[k] == subset_col[k];
367                        ok     |= row[k] == subset_col[k] && col[k] == subset_row[k];
368                        CPPAD_ASSERT_KNOWN(
369                                ok,
370                                "sparse_hes: work is non-empty and conditions have changed"
371                        );
372                }
373        }
374# endif
375        //
376        // check for case where input work is empty
377        if( color.size() == 0 )
378        {       // compute work color and order vectors
379                CPPAD_ASSERT_KNOWN(
380                        pattern.nr() == n,
381                        "sparse_hes: pattern.nr() not equal domain dimension for f"
382                );
383                CPPAD_ASSERT_KNOWN(
384                        pattern.nc() == n,
385                        "sparse_hes: pattern.nc() not equal domain dimension for f"
386                );
387                //
388                // initialize work row, col to be same as subset row, col
389                row.resize(K);
390                col.resize(K);
391                // cannot assign vectors becasue may be of different types
392                // (SizeVector and CppAD::vector<size_t>)
393                for(size_t k = 0; k < K; k++)
394                {       row[k] = subset_row[k];
395                        col[k] = subset_col[k];
396                }
397                //
398                // convert pattern to an internal version of its transpose
399                vector<size_t> internal_index(n);
400                for(size_t j = 0; j < n; j++)
401                        internal_index[j] = j;
402                bool transpose   = true;
403                bool zero_empty  = false;
404                bool input_empty = true;
405                local::sparse_list internal_pattern;
406                internal_pattern.resize(n, n);
407                local::set_internal_sparsity(zero_empty, input_empty,
408                        transpose, internal_index, internal_pattern, pattern
409                );
410                //
411                // execute coloring algorithm
412                // (we are using transpose becasue coloring groups rows, not columns)
413                color.resize(n);
414                if( coloring == "cppad.general" )
415                        local::color_general_cppad(internal_pattern, col, row, color);
416                else if( coloring == "cppad.symmetric" )
417                        local::color_general_cppad(internal_pattern, col, row, color);
418                else if( coloring == "colpack.star" )
419                {
420# if CPPAD_HAS_COLPACK
421                        local::color_general_colpack(internal_pattern, col, row, color);
422# else
423                        CPPAD_ASSERT_KNOWN(
424                                false,
425                                "sparse_hes: coloring = colpack.star "
426                                "and colpack_prefix not in cmake command line."
427                        );
428# endif
429                }
430                else CPPAD_ASSERT_KNOWN(
431                        false,
432                        "sparse_hes: coloring is not valid."
433                );
434                //
435                // put sorting indices in color order
436                SizeVector key(K);
437                order.resize(K);
438                for(size_t k = 0; k < K; k++)
439                        key[k] = color[ col[k] ];
440                index_sort(key, order);
441        }
442        // Base versions of zero and one
443        Base one(1.0);
444        Base zero(0.0);
445        //
446        size_t n_color = 1;
447        for(size_t j = 0; j < n; j++) if( color[j] < n )
448                n_color = std::max(n_color, color[j] + 1);
449        //
450        // initialize the return Hessian values as zero
451        for(size_t k = 0; k < K; k++)
452                subset.set(k, zero);
453        //
454        // direction vector for calls to first order forward
455        BaseVector dx(n);
456        //
457        // return values for calls to second order reverse
458        BaseVector ddw(2 * n);
459        //
460        // loop over colors
461        size_t k = 0;
462        for(size_t ell = 0; ell < n_color; ell++)
463        {       CPPAD_ASSERT_UNKNOWN( color[ col[ order[k] ] ] == ell );
464                //
465                // combine all columns with this color
466                for(size_t j = 0; j < n; j++)
467                {       dx[j] = zero;
468                        if( color[j] == ell )
469                                dx[j] = one;
470                }
471                // call forward mode for all these rows at once
472                Forward(1, dx);
473                //
474                // evaluate derivative of w^T * F'(x) * dx
475                ddw = Reverse(2, w);
476                //
477                // set the corresponding components of the result
478                while( k < K && color[ col[order[k]] ] == ell )
479                {       size_t index = row[ order[k] ] * 2 + 1;
480                        subset.set(order[k], ddw[index] );
481                        k++;
482                }
483        }
484        return n_color;
485}
486
487} // END_CPPAD_NAMESPACE
488
489# endif
Note: See TracBrowser for help on using the repository browser.