"Dener, Alp" <[email protected]> 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.
