Ability to get directional derivatives without multiplying by some seed vector

[modelica-trac-importer]

Reported by Sam.Lishak on 3 Dec 2015 14:50 UTC In order to compute the Jacobian of a model (A, B, C, D matrices), the only option is to call fmi2GetDirectionalDerivative for each state/input and reassemble the columns, as described in section 2.1.9 of the FMI 2.0 standard.

In situations where the full Jacobian is required, this can be needlessly expensive if the underlying model must recompute the entire Jacobian on every call to fmi2GetDirectionalDerivative. It would be a huge advantage if functionality could be added allowing the computation of the full sparse Jacobian in one call, given a list of value references to the known and unknown variables. This is one of the main factors blocking us from adopting the FMI standard: our simulations make very slow progress as they rely on evaluating the full model Jacobian.

---

Migrated-From: https://trac.fmi-standard.org/ticket/360

[modelica-trac-importer]

Comment by mfriedrich on 4 Dec 2015 10:39 UTC As noted in the FMI standard a good implementation of a FMU export should introduce caching of variables!

Hence, when the FMU can only calculate the full Jacobian the FMU should calculate the full Jacobian at the first call to fmi2GetDirectionalDerivative, cache the full Jacobian and return the requested column. All successive calls to fmi2GetDirectionalDerivative then just return the requested column of the cached Jacobian. Of course the FMU must invalidate the cache on any fmi2Set... calls which have influence on the Jacobian. But if the importing tool needs the full Jacobian (the case here) the cache will perform perfectly.

So this is a performance issue of the FMU and should be fixed by the exporting tool. No need the extend the FMI standard.

[modelica-trac-importer]

Comment by Sam.Lishak on 4 Dec 2015 10:46 UTC Thanks for that, we will raise this with the makers of the exporting tool we're using instead.

[modelica-trac-importer]

Comment by otter on 4 Dec 2015 11:14 UTC The suggestion with a cached Jacobian is not the intended way. The principles behind the FMI interface are sketched below:

There are two problems with a direct Jacobian interface:

• There are several widely used ways how a sparse Jacobian is represented. If the FMI interface would pick any of them, there will be users that have to copy the Jacobian to their used data-structure whenever the Jacobian is computed.

• Worse, computing a Jacobian from variables that are recursively computed from each other is hard. Typically, the number of operations and the C-Code to compute it grows very quickly (roughly O(n!^2), where n is the number of variables). For large systems, this is not a practical approach. There is an own branch in mathematics called "Automatic Differentiation" that has developed methods to compute analytic Jacobians in a much more effectiv way. In the FMI 2.0 specification this is shortly sketched on page 27-28 (see copy below). More details are given in the paper http://www.ep.liu.se/ecp/076/018/ecp12076018.pdf. In most realistic cases, this is the most efficient way to compute Jacobians (and the one where the code-size is minimal as well - grows typically with O(n)).

To summarize: The FMI standard only gives a basic interface for partial derivatives. To efficiently compute a sparse Jacobian for a large system, the calling environment has to deal with it (and this requires some work, see below).

---

From FMI 2.0 page 27-28

If the sparsity of a matrix shall be taken into account, then the matrix can be constructed in the following way:

1. The incidence information of the matrix (whether an element is zero or not zero) is extracted from the xml file from element .
2. A so called graph coloring algorithm is employed to determine the columns of the matrix that can be computed by one call of fmi2GetDirectionalDerivative. Efficient graph coloring algorithms are freely available, such as library ColPack (http://www.cscapes.org/coloringpage/) written in C/C++ (LGPL), or the routines by Coleman, Garbow, Moré: “Software for estimating sparse Jacobian matrices”, ACM Transactions on Mathematical Software - TOMS , vol. 10, no. 3, pp. 346-347, 1984. See e.g. http://www.netlib.org/toms/618.
3. For the columns determined in (2), one call to fmi2DirectionalDerivative is made. After each such call, the elements of the resulting directional derivative vector are copied into their correct locations of the partial derivative matrix.

[modelica-trac-importer]

Comment by Sam.Lishak on 10 Dec 2015 11:20 UTC I’m surprised that there is no FMI standard for a sparse Jacobian alongside the standards made for the rest of the model interface. FMI already requires vendors to conform to the standards for parameters, model structure, states etc. Is imposing an FMI standard on a sparse Jacobian format that different? I would have no problems with adopting whatever sparse representation that were used in FMI, or copying this structure to another format if it were more convenient in that situation.

Regardless, the representation of sparse matrices is not the problem we have at the moment, as the models we work with are not particularly large (<100 independent variables) and are not that sparse anyway - the only reason I mentioned sparsity is that I assumed it would be more convenient for it to be returned in that format for those who deal with large sparse systems.

From the point of view of a user who has generated an FMU (in this case, from Dymola) and wishes to use it in a simulation, I'm not sure how I can use techniques such as AD as you mentioned. It seems like your suggestions are directed more towards the modelling tool developers rather than the users of such tools. It is clear that the way Dymola generates the Jacobian can be improved to match the FMI standard better by not computing the entire Jacobian on every fmi2GetDirectionalDerivative call. I'm not sure whether the reason behind this is because the elements are computed recursively, but if they are and that can't be changed then caching seems like the most appropriate option.

As a commercial end user, it is disappointing that it is not more straight forward to get hold of the state space model. My feeling is that there would be a benefit in a method to get either or all of the A, B, C, D matrices out in one FMI command. Indeed, tools such as PyFMI provide convenience methods to do this for you (get_state_space_representation) so it's clearly not that much of an unusual thing to want. For small dense systems, there's no real advantage in employing graph colouring methods and it seems unnecessarily cumbersome to have to make multiple FMI calls just to get the Jacobian.

In the end, even with an efficient colouring routine we are looking at the FMI being at least 20 times slower for our relatively small dense model example than calling the Dymola model directly. We can accept that the model code can be changed to make the FMI call more efficient, but also reiterate that we believe it would be useful to directly retrieve the sparse data via FMI with a prior known structure in one go.

[modelica-trac-importer]

Comment by andreas.nicolai on 10 Dec 2015 18:26 UTC Hi Sam,

you state that the performance may drop significantly (factor 20) when retrieving the Jacobian through several calls instead of doing a simple memory dump (possible, when internal storage structure is known). However, in my benchmark tests - given a good caching implementation (no model evaluation between derivative query function calls) - the performance difference is rather small. And compared to other costs during the simulation, e.g. model evaluation and rest of the integrator code, the pure transfer speed is really marginal.

Therefore, it would help the discussion if you could provide some profiling results that show the need for such improvements. Otherwise we are talking here about a "feeling, that this could be faster" which may not warrant a standard change/enhancement.

The problem with defining a sparsity storage format is that performance-oriented model developers will choose the format that best suites the task: for example, for our FMU that incorporates and Finite-Volume-Solver we use a modified Ellpack-Itpack format, but our building-simulation FMU we have a CSR format. In the room model we use a band structure with additional borders. So we have three matrix formats in three FMUs, all selected for performance reasons. When we were "forced" by the FMU standard to create matrix memory copies into a standardized format before handing out the derivatives, this would effectively ruin the performance gain you wish to achieve by fixing such a storage format.

[modelica-trac-importer]

Comment by Sam.Lishak on 11 Dec 2015 10:07 UTC Hi Andreas,

Our references to poor performance in computing and retrieving a Jacobian specifically refer to FMUs created in Dymola, where no Jacobian caching has been implemented by the developers yet and the Jacobian is computed in its entirety on every fmi2GetDirectionalDerivative call. We've already raised this with the developers separately.

Despite the original request being submitted on performance grounds, it's now more to do with convenience. We would probably not have raised it in the first place had we not experienced the performance issues with Dymola's implementation, as we would have been able to put up with the current FMI method of building the Jacobian up over multiple calls. However, it does seem odd that FMI is focused on a specific method of computing a Jacobian (having a seed vector only really seems relevant to automatic differentiation), whereas there is reluctance to provide a method that gets the whole matrix out because a sparse matrix form cannot be agreed on. It's a bit frustrating to have the ability to calculate all the elements we want out of the Jacobian in one call, only to have them multiplied through by a vector before they're returned.

As previously mentioned, the representation of sparse matrices is not something we're particularly concerned about, although we do appreciate that it's a lot more important in other situations. Perhaps something more flexible that could be considered is the concept of a "seed matrix"? That way fmi2GetDirectionalDerivative could function the same as before, but would be extended such that if you pass in the identity matrix as your seed, you're returned a dense matrix which is the full Jacobian of the variables you requested.

[modelica-trac-importer]

Modified by andreas.junghanns on 6 Feb 2017 10:46 UTC

[modelica-trac-importer]

Comment by aviel on 13 Feb 2017 17:25 UTC Changed the milestone to 2.1, since it seems to be a new feature request.

[modelica-trac-importer]

Modified by cbertsch on 1 Jun 2018 14:20 UTC

[modelica-trac-importer]

Comment by cbertsch on 11 Jun 2018 10:28 UTC (FCP Alignment Meeting)--> future

[chrbertsch]

The sparsity structure of the Jacobian is known and with coloring the Jacobian matrix can be computed efficiently, see Martin's comments above. This takes good advantage of the sparsity. An efficient implementation of the the partial derivatives functions is a tool issue and should be discussed with the tool vendors instead. Regarding getting out the sparse Jacobian matrices out with one call: there are different sparse matrix formats out there, so one would have to select one which would need transforming/copying in most cases anyway (ok, for less elements.)

So following the argument of @MartinOtter above I suggest to close this issue (and possibly discuss getting out the derivative matrices out at one time in the duplicate tickt #645

[chrbertsch]

Regular FMI Design webmeeting

@jph-tavella and Sam: Please provide the concrete wording and API in a PR together with an example.

[jph-tavella]

I recall here the current definition of the primitive fmi2GetDirectionalDerivative in order to reuse in this comment the correct names of the attributes. fmi2Status fmi2GetDirectionalDerivative(fmi2Component c, const fmi2ValueReference vUnknown_ref[], size_t nUnknown, const fmi2ValueReference vKnown_ref[], size_t nKnown, const fmi2Real dvKnown[], fmi2Real dvUnknown[])

Note: within the current FMI 2.0 version, it is assumed the seed vector dvKnown has the same size as the vector vKnown_ref even if this is not explicitly stated.

The proposal to improve the scope of use of the primitive, may be very simple: 1/ We have first to explicitly write that the size of the seed vector is EITHER the size of the vector vKnown_ref (that is to say nKnown) OR the square of the size of the vector vKnown_ref (that is to say nKnown2). 2/ In case the size of the seed vector is nKnown2, the size of the calculated dvUnknown vector is nUnknown*nKnown. 3/ Here is an example in case the size of the seed vector is nKnown**2 (example for an FMU with 3 inputs u1, u2, u3 and 2 outputs y1, y2): vUnknown_ref={y1,y2}, vKnown_ref={u1,u2,u3} dvKnown={1,0,0,0,1,0,0,0,1} dvUnknown is calculated as {dy1/du1,dy2/du1,dy1/du2,dy2/du2,dy1/du3,dy2/du3} 4/ In case the size of the seed vector is nKnown, the size of the calculated dvUnknown vector is nUnknown. 5/ Here is an example in case the size of the seed vector is nKnown (example for an FMU with 3 inputs u1, u2, u3 and 2 outputs y1, y2): vUnknown_ref={y1,y2}, vKnown_ref={u1,u2,u3} dvKnown={1,0,0} dvUnknown is calculated as {dy1/du1,dy2/du1}

[chrbertsch]

@jph-tavella: With your proposal one would have to extend the arguments of fmi3GetDirectionalDerivative by new arguments ndvUnknownand ndvKnown, right?

And does the dvKnown-matrix in your use cases have to be the the unit matrix, or could it be any matrix (with different seed vectors as columns)?

[jph-tavella]

@chrbertsch: That's right for ndvKnown as the value can be either nKnown or nKnown**2. Not sure for ndvUnknown because, as with the primitive in FMI 2.0, the size of the returned vector can be deduced from nUnknown and nKnown when knowing ndvUnknown. And yes in my case the dvKnown vector is derived from the unit matrix but it can be different if needed. So here is the new proposal for the primitive with an s at the end of its name: fmi3Status fmi3GetDirectionalDerivatives(fmi3Component c, const fmi2ValueReference vUnknown_ref[], size_t nUnknown, const fmi2ValueReference vKnown_ref[], size_t nKnown, const fmi2Real dvKnown[], size_t ndvKnown , fmi2Real dvUnknown[])

[jph-tavella]

Sorry: fmi3Status fmi3GetDirectionalDerivatives(fmi3Component c, const fmi3ValueReference vUnknown_ref[], size_t nUnknown, const fmi3ValueReference vKnown_ref[], size_t nKnown, const fmi3Real dvKnown[], size_t ndvKnown , fmi3Real dvUnknown[])

[chrbertsch]

FMI Design Webmeeting:

Masoud: In case of variables that vectors or matrices, this gets very complicated Markus: We should keep the getDirectionalDerivatives as it is. We could add an getAdjointDerivatives function which might be beneficial in some cases, we could add a function getJacobian function, but do not see the need We should not add a matrix of seed vectors to the API Torsten: in languages like Python or Java, function calls are expensive. Klaus: one could combine the calls to the individidual getDirectionalDerivatives in a C++ function and call this from Pyhton etc. A call to a getJacobianMatrix would not make so much sense if we do not have a sparse format. One can use the sparsity with existing methods. Klaus: this feature could also be added later in FMI 3.1

Agreeement, that we do not want the seed matrix solution.

Poll: Do we want to add a new optional getJacobianMatrix function? Pro: Contra: Christian, Masoud, Klaus, Otto, Markus Abstain: Andreas P., Henrik, Torsten