[CppAD] Sparse Hessians and Jacobians

Wed May 9 20:09:54 EDT 2012

Personally, I would at least want to be able to pass reference to allocated memory.  The CppAD definition of a simple vector can be restrictive.  For example, most of my mathematical computation is done with the Eigen numerical library, and it would be good if I could just pass a reference to the start of an Eigen vector.  Now, I was able to write an Eigen plug-in that makes Eigen vectors qualify as CppAD simple vectors.  But that may not always be the case for everyone.

Although passing references for output variables would be different from how CppAD works now, it is not all that different from many C libraries (e.g., GSL, Intel CBLAS, etc).  Returned values might make more sense if they are returned as expressions, as in Eigen.

Nevertheless, in general I think it would be useful to minimize memory allocation within CppAD functions, and let the user re-use blocks of memory that are allocated only once.

I'm looking forward to testing whatever you come up with.  Thanks again for all of your work with CppAD.

On May 9, 2012, at 2:51 PM, <bradbell at seanet.com>
 <bradbell at seanet.com> wrote:

>> If we only want to support one format, this triplet format can work.  It
>> is intuitive, and there are routines available to convert to and from CSC
>> and CSR formats.
>> 
>> I do have a question about how the Jacobian is returned.  If I understand
>> the code correctly, CppAD allocates new memory for a new Jacobian/Hessian
>> every time the function is called.  If I am repeatedly computing gradients
>> and Hessians, this is an awful lot of dynamic memory allocation.  Wouldn't
>> it be more efficient for the user to allocate memory for the output, and
>> then pass in a pointer or reference for that array?
>> 
>> In that case, the call would be something like:
>> 
>> f.SparseJacobian(x, p, r, c, *out)
> 
> Yes, this would be more efficient, but it would be a marked change from
> the other cases where the value is returned. Perhaps one should have the
> option to pass in the allocated memory as you suggest above.
> 
> Perhaps all the CppAD routines should have an interface that uses input
> and output iterators (for arbitrary containers), instead just simple
> vectors. Harder for the beginning user, but more versatile, more
> efficient, and expected by the C++ expert.
> 
>> 
>> Also, will the SparseJacobian / SparseMatrix functions still need to
>> allocate a full m-by-n matrix internally?  The value of a sparse matrix
>> algorithm is not just that the matrix is return in a sparse format, but
>> the algorithm itself exploits the sparseness in terms of its memory usage.
> 
> The sparse Jacobians and Hessians only allocate space corresponding to one
> row (forward mode) or one column (reverse mode). Then (using the sparsity)
> multiple rows or columns are computed as once and then separated.
> 
> On the other hand, the one row or one column is for every intermediate
> variable in the operation sequence; see size_var
>      http://www.coin-or.org/CppAD/Doc/seq_property.htm#size_var
> not just the independent and dependent variables.
> 
>> 
>> 
>>> 
>>> 
>>>  1. Sparse Hessians and Jacobians (Brad Bell)
>>> 
>>> 
>>> ----------------------------------------------------------------------
>>> 
>>> Message: 1
>>> Date: Tue, 08 May 2012 07:45:05 -0700
>>> From: Brad Bell <bradbell at seanet.com>
>>> To: CppAD at list.coin-or.org
>>> Subject: [CppAD] Sparse Hessians and Jacobians
>>> Message-ID: <4FA93171.5040207 at seanet.com>
>>> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>>> 
>>> I am considering adding sparse Jacobian and Hessians returns to CppAD as
>>> per the sparsity discussion under
>>>    http://list.coin-or.org/pipermail/cppad/2012q1/thread.html
>>> 
>>> Reading over the discussion, it seems to me that the format used by
>>> ADOL-C is reasonable; i.e. three vectors rind[], cind[], and values[]
>>> all of the same length could be used to represent Sparse Hessians and
>>> Jacobians (where the corresponding element types are size_t, size_t, and
>>> Base).
>>> 
>>> For example, the sparse Jacobian driver
>>>    http://www.coin-or.org/CppAD/Doc/sparse_jacobian.xml
>>> could be extended to have the syntax
>>>    jac = f.SparseJacobian(x, p, r, c)
>>> where
>>>    x is the point where the Jacobian is being evaluated
>>>    p is the CppAD sparsity pattern (for the entire Jacobian)
>>>    r is the row indices
>>>    c is the column indices
>>>    jac is the Jacobian values at the corresponding row and column
>>> indices.
>>> 
>>> Note that the sparsity pattern p can be computed by either
>>>    http://www.coin-or.org/CppAD/Doc/forsparsejac.xml
>>>    http://www.coin-or.org/CppAD/Doc/revsparsejac.xml
>>> Then the users can decide which non-zero indices they want back and in
>>> what order. For example, if the Jacobian is know to be symmetric, one
>>> only needs to ask for the upper or lower triangle.
>>> 
>>> A similar approach be used for the sparse Hessians.
>>>    http://www.coin-or.org/CppAD/Doc/sparse_hessian.xml
>>> 
>>> 
>>> 
>>> ------------------------------
>>> 
>>> _______________________________________________
>>> CppAD mailing list
>>> CppAD at list.coin-or.org
>>> http://list.coin-or.org/mailman/listinfo/cppad
>>> 
>>> 
>>> End of CppAD Digest, Vol 56, Issue 1
>>> ************************************
>> 
>> Michael Braun
>> MIT Sloan School of Management
>> braunm at mit.edu
>> --------------------------------------------------
>> View my research at
>> http://braunm.scripts.mit.edu/
>> 
>> _______________________________________________
>> CppAD mailing list
>> CppAD at list.coin-or.org
>> http://list.coin-or.org/mailman/listinfo/cppad
>> 
> 
> 

-------------------------------------------
Michael Braun
Associate Professor of Management Science (Marketing Group)
MIT Sloan School of Management
100 Main St.., E62-535
Cambridge, MA 02139
braunm at mit.edu
617-253-3436

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://list.coin-or.org/pipermail/cppad/attachments/20120510/f5c4fa63/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: smime.p7s
Type: application/pkcs7-signature
Size: 1844 bytes
Desc: not available
URL: <http://list.coin-or.org/pipermail/cppad/attachments/20120510/f5c4fa63/attachment-0001.p7s>