Rationales behind the creation of PyOp

[ghisvail]

Dear Ryan and Daniel,

Could you explain in your README, the rationales behind the creation of yet another linear operator package ? To do so, you might want to answer these specific questions:

• Acknowledge the existence of "competing" solutions of the linear operator design, such as `scipy,linoporpyoperators``,
• Detail how your package is different from these existing solutions,
• Explain what prevents your from contributing to these packages,

Best regards,

• Ghis

[argriffing]

Hey guys (@pv, @jakevdp, @ghisvail, @larsmans, @dwhensley, @rdodesigns, @ogrisel, @MechCoder, @pchanial, ...) it seems really obvious that the state of python linear operators development should be improved somehow. Ideally @pchanial's original scipy PRs many years ago would have been merged and this discussion would be moot. But even if there is not a lot of spare time or energy to improve this corner of scipy, would it be worth identifying a place on the internet to put this discussion instead of having it scattered across various mailing lists and github repos and personal emails? I personally don't have experience with this kind of social coordination but maybe some of you guys do?

[rdodesigns]

So the tl;dr is we would be totally stoked to work towards a common concept and place our imaging related stuff on that base, whatever the base is. I can setup up a mailing list if that works for everyone (@argriffing I am not sure if this is best either but we could try it?)

I'll write up what @ghisvail requested tonight, but I wanted state why we have this particular package. We needed something more than the scipy operators (which did not have simple composition with respect to transposing an operator) and we were focused on image reconstruction for medical imaging. At the time, we unfortunately did not know about PyOperators (it looks great), so our image reconstruction code is written around PyOp.

For us, the important features were

• block composition (for 3D image data sets),
• composition with transposition (for gradient descent type algorithms),
• ease of operating on nD arrays (the data),
• validation through the adjoint test (required for some mathematical guarantees to hold) and shape matching (making sure the output is the correct shape), and
• the ability to substitute numpy.ndarrays with functional linear operators transparently (for modularity/testing).

For example, we could easily define functions on images themselves without having to deal with the reshaping of the vector often or dealing with 1D array inputs (which are assumed to be column vectors).

h, w = (10, 10)
op_shape = (h*w, h*w) ## Say the domain and codomain are the same size

## matvectorized pulls out each column of the input matrix, reshapes it into 
## the shape of the image, and then flattens the result to a column in the 
## resulting matrix. Order is how the input column should be reshaped.
@matvectorized((h, w), order = 'F')
def forward(img):
    ## Operation on an img, where img is already in it's correct shape (2d/3d/nd)

## The `-1` is the same as for the reshape command.
@matvectorized((h, -1), order='F')
def adjoint(co_img):
    ## Operation on the codomain image, where co_img is already in it's 
    ## correct shape (2d/3d/nd)

A = LinearOperator(op_shape, forward, adjoint)

## Test the operator with a bunch of random input.
for _ in range(1000):
    adjointTest(A)

## Turn the operator to work in the Fourier domain instead of the image domain, with 
## automatic fft, fftshift, ifftshift, ifft (in that left to right order).
B = fftwrap(A, shift='all')

### Gradient Descent
y = # some initial data
x = # initial guess to the solution

alpha = # some step control parameter.

## Stepping in the projected Landweber algorithm, with projection function p.
new_x = p(x - B.T.dot(B.dot(x) - y)))

With the aforementioned features a lot of our operators became only a handful of lines that could be tested on the images themselves, instead of the vectorized versions. The shape checking on the output kept us honest.

The stuff in PyOp do not strongly depend on the underlying representation of the linear operator though, so the functions our lab uses could be easily modified to use the structure of other code out there. We just wanted something that the students in the lab could use/modify easily; for us this meant a flat (single class) functional/FauxO style but that design choice was mostly due to prior experience.

[ghisvail]

Hi all,

I did not mean to start such a long thread here. I was only genuinely interested to know what PyOp was bringing to the current linear operator offering.

I too happen to work in the field of medical imaging (MRI reconstruction), and the prospect of wrapping complicated operations into simple LinearOperator-like interface is particularly important to raise the level of abstraction in our code and build smarter tools. Judging by the feature set you guys wanted to implement, I felt that it was no different than what is achieved in linop, hence me asking for clarifications.

Regarding what @argriffing raised, I feel that scipy, linop and pyoperators can definitely coexist since they solve different things. Scipy lead the way to provide a common interface, but the heaviness of the scipy stack can be a hindrance for whatever application, which only needs this interface. This is what linop actually provides, a separate package with a scipy-like LinearOperator interface, some derived classes for common mathematical operators and support for composition. Then, from my personal experience, someone who wants smarter operators should move to pyoperators, which features smarter compositions, introspection of memory usage, parallel computing and more...

EDIT: fixed some typos