Rather than accumulating callbacks, couldn't you keep an array of operators and devectorize the computation?
A potential drawback would be a loss of extensibility... You could also cache the realized operator after the first evaluation (memoization). —Pierre-Yves On Saturday, May 24, 2014 9:30:43 PM UTC+2, Dominique Orban wrote: > > I just put together a basic linear operator package for Julia: > https://github.com/dpo/linop.jl > > The central idea is that operators act like matrices, but it would be too > expensive to combine them before they are actually applied to a vector. > Another benefit is that you can manipulate them as if they were matrices > and write expressions such as y = (A' * A + B) * x without worrying that A' > * A + B will actually be formed (it won't). > > There are basic tests, but more should be added, and lots is still missing > (e.g., block operators, parallel operators, ...). Entire algorithms could > be wrapped as operators so that conceptually, you could write things like: x > = LSQR(A) * b to apply the LSQR algorithm to solve the least-squares > problem with A and b. > > Feedback most welcome. > > Cheers. >
