JuliaImageRecon / RegularizedLeastSquares.jl

MIT License
20 stars 9 forks source link

Provide the Option to supply only the Normal Operator [breaking] #63

Closed JakobAsslaender closed 9 months ago

JakobAsslaender commented 9 months ago

Hi @tknopp,

this is a first draft of what I mentioned to you. You are absolutely right that I could have hacked the old version by supplying an identity operator as A. But this really feels like a hack to me, in particular, as it requires intimate knowledge of each algorithm and whether it uses A for anything but the back-projection. The current version allows you to choose whether you want to supply A or AHA or both where this makes sense and makes A compulsory for algorithms that require it.

While I was at it, I also tried to tidy up and homogenize the interface a bit. In particular, I

I left the createLinearSolver interface unchanged, i.e., it softens the breaking aspect of this change but also limits the power of this interface. I never use it as I rather use the constructors directly, so I don't have any strong feelings about it. Not sure if having two interfaces helps or hurts in terms of clarity.

Things I still plan to do:

Thoughts? I am aware that I took a rather radical re-write approach to the interface and am sure I have overlooked a thing or two in terms of functionality you had in mind. I would appreciate any feedback you might have! Also happy to have Zoom call for a quicker convergence to mutual understanding.

Thanks for considering!

nHackel commented 9 months ago

Hello! Tobi, Mirco and I talked a bit about the PR. Overall I think all the implemented and planned changes are very good, we just have two small-ish points.

The first is in regards to the kwargs. We stumbled over a similar issue with the dangerous "sink" behavior in our MPIReco.jl package (the "fix" there actually started our effort in updating RegLS lately). Having just an explicit list of keywords is definitely safer, but it would cause a few issues in MPIReco and MRIReco. However, I think I found a solution where we can have "both":

The direct constructors only use an explicit list of keyword arguments and instead the createLinearSolver function receives any number of keyword arguments and just selects the fitting ones for the given solver. This requires a little bit of introspection, where we look at the method table for the constructors with methods(Kaczmarz) for example and then with Base.kwarg_decl.(...) we can look at the list of keywords.

This way we can still have the safety when directly using a solver and the flexibility if the utility function is used. createLinearSolver could even give (debug or optional) logging messages about unused parameters. I can add this functionality later on top of this PR.

The second comment is about the constructor and if A should be an argument or a keyword argument. Since A is an "understandable" and often times mandatory argument for all solvers, even if it is just used to create the normal operator in the default case, I feel like it would better to have it be an argument instead of a keyword.

If we then define the function as follows: FISTA(A=nothing; AHA = ..., reg = ...) then you can still call FISTA(AHA = foo, ...) and only supply the normal operator. I haven't tried it out yet, but I think that should cover your use case.

Best regards, Niklas

JakobAsslaender commented 9 months ago

Hi @nHackel et al.,

thanks for your fast and positive response!

Regarding kwargs... and createLinearSolver: That sounds like a good solution! What is your main goal with that function? Is it mostly an interface for MRIReco.jl and MPIReco.jl? I was wondering how public-facing this function should be. I find it more intuitive and transparent to use the individual constructors as the public-facing interface. To this end, I was going to suggest modifying the documentation: here is an example of how to create a TV regularization object, bake it into an ADMM regularizer and call solve!. And then just put the API information with all the regularizers solvers. One step further, we could remove the export of createLinearSolver and make it a function that is only intended for package developers who need additional functionality.

Regarding A: I hear you. I rewrote the interface to take A as a non-key argument. The nice thing about this is that we can dispatch the two versions, making error messages more transparent. I implemented two constructors, e.g., FISTA(A; AHA = A'*A, rho = 0.95, ...) and FISTA(; AHA, rho = 0.95, ..., i.e., AHA is a required keyword if A is not supplied.

Let me know what you think about my suggestion regarding the documentation. Happy to modify it accordingly.

-ja

JakobAsslaender commented 9 months ago

Fixes Issue #54 .

nHackel commented 9 months ago

Hello @JakobAsslaender, I really like the two constructor version for A and AHA. That's a great idea!

Regarding createLinearSolver: I think our main goal here is to have an "easy" constructor for a solver, that can take a hodgepodge of parameters and return a "best-effort" solver. For MPIReco/MRIReco the hodgepodge is a result of the reco frameworks supporting the "union" of all solver parameters in a way and not the solver specific sets. This is partially a result of the least-squares details being hidden for the high-level users.

I think another use case we could consider is if a user is not entirely sure which solver is the correct one for their given least-squares problem. We could add a createLinearSolver version where the solverType is not specified and instead the function checks which solver are actually applicable for the given A, AHA, x and reg via the isapplicable function.

If a user is sure which solver they want to use they can directly construct it and I think we can reflect that like you said by directly building the solver in the documentation. And maybe have an additional small section for the createLinearSolver in the docs.

I'll try to take a look at the code today or tomorrow. Did you get around to all the changes you want to do?

JakobAsslaender commented 9 months ago

I think it would be great to have the function createLinearSolver to pick an algorithm based on the regularization!

I think I am done with the changes I had in mind. From my end, this branch can be merged if no one has any objections. We can wait for your additions, but I guess this could also be a separate branch and if we decide to make a release, this change would be breaking and your changes shouldn't be, so I would be fine with releasing the code as is and making a separate one for you changes if you don't want to rush it.