[petsc-dev] MATLMVM

[Jed Brown]

"Dener, Alp" <adener at anl.gov> writes:

> Sure. Symmetric Broyden is defined as the convex linear combination of BFGS and DFP updates such that SymBrdn = (1-phi)BFGS + phi*DFP with phi in range of [0, 1]. It is possible, mathematically, to produce both BFGS and DFP applications out of the single Symmetric Broyden object. However, in the two extreme cases where phi = 0 and phi = 1, the pure BFGS and pure DFP updates can be implemented much more efficiently than the Symmetric Broyden formula. Specifically, the BFGS inverse action and the DFP forward action can be unrolled into a two-loop algorithm that significantly reduces dot products and and AXPBY operations compared to the Symmetric Broyden formula with the matching phi scalar. Furthermore, BFGS and DFP implementations require about 25% fewer storage vectors than Symmetric Broyden because they don’t need to carry around the extra information required for the other extreme. Similar benefits apply to the SR-1 method, which is actually also part of the Symmetric Broyden family, except with a dynamically calculated phi value (based on dot products with update vectors) that lie outside the range of [0, 1]. Implementing it on its own allows for explicitly cancelling out some terms in the update, reducing it to a rank-1 formula instead of the rank-2 formula in Symmetric Broyden, and apply it far more efficiently.
>
> It’s possible for all of these to be offered inside Symmetric Broyden with checks on a subtype or the phi value that direct the operations into the correct special case, but that doesn’t reduce the amount of code that needs to be written to maintain the efficiency. And if you wanted to really have minimal code that does everything out of the Symmetric Broyden update formula just to have minimal amount of code, then the special cases would be performing a lot of unnecessary operations that could (and should) have been eliminated.

Thanks.  By analogy, VecAXPBY_Seq has a number of special cases, rather
than depend on the user to choose the best interface themselves.

I think we'd be better off with one matrix type and let the
implementation handle the special cases.  But I understand how that goes
more naturally with making the MatLMVM a strict matrix instead of a
container for optimization, and I don't know if changing any of it
should be a priority right now.