2
$\begingroup$

I need to calculate the Jacobian matrix of a subroutine F(U). Both F and U are of size N(=O($10^5$)). Using Tapenade, I differentiated the routine in tangent mode. I cannot calculate the full Jacobian directly because of the large memory requirement.

  1. Is there a package that I can use to parse the differentiated FORTRAN routines to get the sparsity pattern? preferably in python.

  2. I plan to naively construct the Jacobian matrix using N tangent calculations, using sparsity pattern to reduce the memory usage. Is there a better way of doing this?

I would like to stick with Tapenade because of other issues.

Improve this question
$\endgroup$
5
  • $\begingroup$ Would it work to construct the sparsity pattern on the fly, by reading where the zeros are in the individual calculations? E.g., you could transform each individual derivatives column to a compressed column format, then write it to the Jacobian stored in compressed columns. That way no precomputed sparsity pattern is required. It's probably easier to change storage format afterwards. $\endgroup$
    – Kirill
    Aug 13, 2015 at 3:44
  • $\begingroup$ Yes, this is something I plan to do, but this is not the cleanest way though. If I don't know the number of non zero elements before hand, I will have to declare an array which is big enough to store them. Also, a zero value of an element does not imply independence(just that it is zero for those set of input parameters). EDIT: zero value and independence is an issue if I want to reuse the pattern. $\endgroup$
    – 0b1100001
    Aug 13, 2015 at 3:57
  • $\begingroup$ I was thinking that's a much more trivial programming task than the general pre-computed solution, if you were to do it yourself. I meant you would dynamically construct the Jacobian each time from scratch: given that you're already computing $O(N^2)$ elements, the overhead of going over them one by one and putting everything together would be negligible. Is the ability to reuse the pattern so important? Constructing it is fairly cheap. $\endgroup$
    – Kirill
    Aug 13, 2015 at 4:02
  • $\begingroup$ That's right. I should probably do this. I was being a little ambitious here, how about inserting the sparsity pattern in the code as constant and compiling at run time and let compiler optimize and remove all the multiplications by zero. $\endgroup$
    – 0b1100001
    Aug 13, 2015 at 4:11
  • $\begingroup$ So I tried this, and it is really slow. Takes around 5 hours to build the Jacobian. This weekend, I am planning of parsing the tapenade generated code using python and finding some approximate dependency. Essentially there are lots of FOR loop inside dFdU.DX. If i can decrease there range, I could get some performance improvement. $\endgroup$
    – 0b1100001
    Aug 14, 2015 at 17:39

4 Answers 4

Reset to default
1
$\begingroup$

I need to calculate the Jacobian matrix of a subroutine F(U). Both F and U are of size $N(=O(10^5))$. Using Tapenade, I differentiated the routine in tangent mode. I cannot calculate the full Jacobian directly because of the large memory requirement.

There are AD (automatic differentiation) software packages that can do this calculation; I know DAEPACK can, because it came out of my PhD advisor's group. It appears that Tapenade does not have sparsity analysis capabilities.

  1. Is there a package that I can use to parse the differentiated FORTRAN routines to get the sparsity pattern? preferably in python.

My intuition is that it probably doesn't exist, because it would require a lot of static analysis. A cursory Google search backs that up.

  1. I plan to naively construct the Jacobian matrix using N tangent calculations, using sparsity pattern to reduce the memory usage. Is there a better way of doing this?

Probably. It really depends on what you want to do with the Jacobian matrix. If you mostly want to solve linear systems that are relatively well-conditioned or are amenable to preconditioning, I'd suggest instead considering matrix-free approaches, which would allow you to make use of the tangents you calculate. A Jacobian-free Newton-Krylov approach might even be viable, which would make AD unnecessary.

However, if you really need the entire Jacobian matrix, I'd suggest sparse AD instead.

Improve this answer
$\endgroup$
3
  • $\begingroup$ I really want to stick with Tapenade, because the code is quite old and fragile and getting the adjoint to work took some time. Even with Tapenade I had to manually change some parts to make it work. I will probably look into Jacobian free Newton Krylov approach. $\endgroup$
    – 0b1100001
    Aug 14, 2015 at 17:36
  • $\begingroup$ one minor comment: tangent calculations via AD are not subject to finite differencing errors. The derivatives (i.e., the Jacobian entries) can be expected to have the same accuracy as the original function (F). $\endgroup$
    – GoHokies
    Aug 27, 2015 at 12:13
  • $\begingroup$ Good point. I think I misinterpreted the OP when it said "I plan to construct the Jacobian matrix using N tangent calculations", not realizing they were still planning on using AD. I will update my answer. $\endgroup$
    – Geoff Oxberry
    Aug 27, 2015 at 18:00
3
$\begingroup$

If you were using Julia, then SparsityDetection.jl solves (1) by using non-standard analysis with concolic fuzzing to directly determine the sparsity pattern from the function's AST, and SparseDiffTools.jl performs matrix coloring and colored automatic differentiation to do a sparse calculation of the Jacobian to solve (2).

There are other systems that can do this, just not in a high level language. TAF is a good one in Fortran, DAEPACK as well, and ADOL-C can do this in C++. These all detect sparsity and use matrix coloring techniques. I do not know of a solution in Python.

Improve this answer
$\endgroup$
0
$\begingroup$

Regarding item 2 of your question, assuming you have the sparsity pattern, it is often possible to compute the Jacobian with far fewer than N evaluations of the function. The idea is that because the Jacobian is sparse, you can evaluate more than a single column with a single function evaluation. Take a look at this paper for the details:

Software for Estimating Sparse Jacobian Matrices

Improve this answer
$\endgroup$
0
$\begingroup$

A cursory Google search revealed this resource.

It is a set of FORTRAN subroutines for determining the Jacobian sparsity pattern of a general multivariate function. What you will want to look at is subroutine RP01A:

C This subroutine automatically generates the sparsity pattern
C of a sparse Jacobian matrix with an arbitrary or band structure.

There is a list of numerical examples (albeit for very small N) that you can use as a starting point for your own problem. I am not sure that their algorithm can handle something like $N = 10^5$, but I think it's worth giving it a shot.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.