pylops / pyproximal Goto Github PK

Now that the pylops FISTA works on matrices of data, it makes sense that the regularization parameter can be a vector of values (one per vector of data in the data matrix).

Implementing this feature will require that the proximals work at the matrix level rather than the vector level.

Edit: I actually think Python would already take care of making proximals works with matrices.

In the blind deconvolution method, I would like to update the point spread funciton after updating the data term with ADMM using LinearizedADMM.
Does Pyploximal have a function similar to deconvblind that returns psf in Matlab?
https://jp.mathworks.com/help/images/ref/deconvblind.html?lang=en

For reference: openjournals/joss-reviews#6326

Here are some comments I have regarding the paper itself:

• In the statement of need, it would be helpful to cite a few well-known packages in order to highlight the key distinguishing features of this package. This is already done in the repo README and can probably be reused.
• One DOI is missing -- see here for the suggested DOI

This is part of the review at openjournals/joss-reviews#6326

The installation documentation mentions

git clone https://github.com/mrava87/pyproximal.git

While GitHub redirects links it may be helpful to update this to reflect the new ownership.

Don't know what's going on behind the scenes but skimage TV is a fair bit quicker than your implementation. Could be looped in with something like the following:

from pyproximal.ProxOperator import _check_tau
from skimage.restoration import denoise_tv_chambolle, denoise_tv_bregman

class TV:
    def __init__(self, sigma, dims, isotropic=True, max_iter=1000):
        self.sigma = sigma
        self.isotropic = isotropic
        self.max_iter = max_iter
        self.count = 0
        self.dims = dims

    def _increment_count(func):
        """Increment counter"""

        def wrapped(self, *args, **kwargs):
            self.count += 1
            return func(self, *args, **kwargs)

        return wrapped

    @_increment_count
    @_check_tau
    def prox(self, x, tau):
        x = x.reshape(self.dims)
        if self.isotropic:
            return denoise_tv_chambolle(
                x, weight=tau * self.sigma, max_num_iter=self.max_iter
            ).ravel()
        else:
            return denoise_tv_bregman(
                x,
                weight=tau * self.sigma,
                isotropic=self.isotropic,
                max_num_iter=self.max_iter,
            ).ravel()

Hi,

Great work on this package. I’ve actually started building around it here https://github.com/jameschapman19/scikit-prox using your operators to solve regularised versions of some scikit-learn models. I like how you’ve chosen the proximal operators as the level of abstraction it has made it really easy to work with.

My question is whether you have or would consider combined proximal operators by Dykstra or otherwise. It’s been done in https://github.com/neurospin/pylearn-parsimony.

The nuclear norm is the soft-thresholding operator acting on the singular values. This essentially means that even the larger singular values are penalized just as hard as the smaller ones, which leads to the phenomenon shrinking bias which has been studied in e.g. [1]. Ideally, the smaller singular values (likely stemming from noise) should be penalized proportionally harder, and this can be done in two ways:

Weighted nuclear norm (WNN)

The weighted nuclear norm adds a penalty w_i for the i:th singular value, where {w_i} is increasing. [2]

Concave function acting on the singular value

There are many works that add instead apply a concave map f to the i:th singular σ_i -> f(σ_i). These include:

• SCAD [3]
• Log [4]
• MCP [5]
• ETP [6]
• Geman [7]
• f_mu (special case of #51)

Proposition for enhancement

Implement the above penalties and their proximal operators.

References

[1] Carlsson et al. "An unbiased approach to low rank recovery", https://arxiv.org/abs/1909.13363
[2] Gu et al. "Weighted Nuclear Norm Minimization with Application to Image Denoising" In the Proceedings of the Conference on Computer Vision and Pattern Recognition (CVPR), 2014
[3] Fan and Li "Variable selection via nonconcave penalized likelihood and its oracle properties" Journal of the American Statistical Association, 96(456):1348–1360, 2001
[4] Friedman "Fast sparse regression and classification", International Journal of Forecasting, 28(3):722 – 738, 2012.
[5] Zhang et al. "Nearly unbiased variable selection under minimax concave penalty" The Annals of Statistics, 38(2):894–942, 2010
[6] Gao et al. "A feasible nonconvex relaxation approach to feature selection", In Proceedings of the Conference on Artificial Intelligence (AAAI), 2011
[7] Geman and Yang "Nonlinear image recovery with half-quadratic regularization", IEEE Transactions on Image Processing, 4:932 – 946, 08 1995

There are now an increasing amount of penalties being merged. Currently on main:

• Box
• Simplex
• Intersection
• AffineSet
• Quadratic
• Euclidean
• EuclideanBall
• L0Ball
• L1
• L1Ball
• L2
• L2Convolve
• L21
• L21_plus_L1
• Nonlinear
• Huber
• Nuclear
• NuclearBall
• Orthogonal
• VStack
• SCAD
• Log
• ETP
• Geman
• QuadraticEnvelopeCard

I think it would be beneficial to categorize/group these penalties. Today we only have this flat structure, but new users could get confused of what penalties can be applied to what problems. There are many ways of structuring these, but as a first step I suggest

• First level: Separate between vector penalties and matrix-only penalties.
• Second level (perhaps not yet): Separate between convex and non-convex penalties

Any thoughts? They should also be sorted in an alphabetical order in the documentation.

Hi, and thank you for a nice project.

I actually started implement a framework similar to yours, but I am lucky that I found pyproximal, especially considering the similarities in most of the design choices you have made. No need to write things twice, so I am instead hoping I can contribute a bit.

What I was mostly discouraged by was how slow a simple ADMM lasso implementation was using your framework. When compared to my own implementation it was 15x slower, but I could quickly make it faster using two things:

1.) Change scipy.linalg.solve -> np.linalg.solve for dense matrices in LinearOperator (in the twin pyLops library), and
2.) Cache the Op1 operator in L2.prox()

I understand that you might want to keep the scipy.linalg.solve for backwards compatability, but I urge you to make it configurable in some way, since it is really slow. I would be happy to create PRs if you are interested.

Again, thank you for a nice framework.

This is part of the review at openjournals/joss-reviews#6326.

On this page, the line

For this example we will use the well-known Rosenbrock function:

should be followed by a mathematical description of the Rosenbrock function, but only a generic function is presented.

Motivation

PyProximal solvers are currently only stopped when the maximum number of iterations is reached.

It is desirable to introduce stopping criteria to prematurely stop a solver based on some condition. A first attempt at introducing a stopping criterion to ProximalPoint and ProximalGradient based on the change in objective function between two consecutive iterations can be found in #176

Definition of done

The same criterion should be implemented to all solvers

• More of such criteria should be implemented (e.g., operating only on one term of an objective function instead of the overall objective function)

Hello everyone,
given that you have contributed some code to PyProximal, I would like to give you the opportunity to be included as co-author of this JOSS paper. I am writing it mostly to ensure PyProximal can be properly cited by researchers leveraging it in their published work going forward.

I kindly ask you to thumb up if you would like to be included in the author list. If I do not see a thumb up I will assume that you are not interested.

Moreover, to make the paper stronger, I would like to add a couple of more use-cases of PyProximal: ideally they should be backed up by peer-reviewed publications (and/or openly available codes). It would be great if you could provide me with a use-case from your team or other people you know in your field.

To access the paper please head over to https://github.com/mrava87/pyproximal/actions/runs/7271942518 and click on the paper text under Artifacts (the PDF will start downloading). The paper .md file is currently sitting in this branch: https://github.com/mrava87/pyproximal/tree/joss. I am happy to push it to the main repo in a branch with the same name if you want to act directly on it, or just take any of your comments/additions and add it to the paper myself.

Contributors:
@marcusvaltonen
@NickLuiken
@olivierleblanc
@eurunuela
@shnaqvi

There is already support for BilinearOperator in pyproximal/pyproximal/utils/bilinear.py and the PALM optimizer; however,
they do not scale to second-order methods such as Levenberg-Marquardt (LM) and Variable Projection (VarPro). In recent years, such methods have been become popular alternatives to splitting methods such as ADMM for certain types of problems (e.g. image analysis and computer vision).

Examples

As an example the nuclear norm penalty has a variational formulation

which was utilized in [1] using a bilinear framework (no second order method involved, though).

This also scales to the weighted nuclear norm, and it has been successfully implemented for vision tasks using second-order methods in e.g. [2, 7, 8]. Other non-convex penalties have been studied as well, e.g. [3, 4, 5, 6].

Suggestion on enhancement

I suggest adding support for Levenberg-Marquardt (LM) and Variable Projection (VarPro). As a second step I suggest adding some of the mentioned regularizers.

References

[1] Cabral et al. "Unifying Nuclear Norm and Bilinear Factorization Approaches for Low-Rank Matrix Decomposition", In the Proceedings of the IEEE International Conference on Computer Vision (ICCV), 2013
[2] Iglesias et al. "Accurate Optimization of Weighted Nuclear Norm for Non-Rigid Structure from Motion", In the Proceedings of the European Conference on Computer Vision (ECCV), 2020
[3] Valtonen Ornhag et al. "Bilinear parameterization for differentiable rank-regularization", In the Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPR), 2020
[4] Valtonen Ornhag et al. "Differentiable Fixed-Rank Regularisation using Bilinear Parameterisation", In the Proceedings of the British Machine Vision Conference (BMVC), 2019
[5] Valtonen Ornhag et al. "Bilinear Parameterization for Non-Separable Singular Value Penalties", In the Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2021
[6] Xu et al. "A unified convex surrogate for the schatten-p norm". In Proceedings of the Conference on Artificial Intelligence (AAAI), 2017
[7] Kumar "Non-rigid structure from motion: Prior-free factorization method revisited" In IEEE Winter Conference on Applications of Computer Vision (WACV), 2020
[8] McDonald et al. "Spectral k-support norm regularization", In Advances in Neural Information Processing Systems (NIPS), 2014

Is it possible to apply two convolution filters to the variable X, compare the elements of the two arrays, and choose the larger value?　And can pyproximal　solve the problem?

Currently the ADMM solver accommodates only for the special case of A and B identities.

We should extend the solver to more general cases and have A=None and B=None as defaults (which means identities)

I have a large-scale problem where I use the L1-regularized least square for compressed sensing by Stephen Boyd . However, the computation takes something like 170 seconds, which drove me to look for something more efficient. I found that also ISTA, FISTA, and Twist are able to fulfill the same task, but I did not find an implementation of the L1-regularized least square in the proximal library.

Q1: Is the L1-regularized least square already exist in the library?
Q2: Are ISTA, FISTA, and Twist capable of solving large-scale problems (A matrix shape (149600, 144))?
Q2: is the ADMM solver capable of solving such problems? If yes, is it more efficient?

In [1] a general low-rank inducing framework was introduced and showed many benefits compared to e.g. the nuclear norm when applied to common computer vision tasks. This was later generalized in [2].

Suggestion on enhancement

Implement [2] and treat the special cases R_g, R_mu and R_i=k, separately. (See papers for details).

References

[1] Larsson and Olsson, "Convex Low Rank Approximation", International Journal of Computer Vision (IJCV), 2016
[2] Valtonen Ornhag and Olsson, "A unified optimization framework for low-rank inducing penalties", In the Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2020

Add pydata-sphinx-theme in environment-dev.yml

As discussed in #143, we should implement the Proximal Average Proximal Gradient solver to handle multiple proximals.

See https://proceedings.neurips.cc/paper/2013/hash/49182f81e6a13cf5eaa496d51fea6406-Abstract.html

The deprecation of Python 3.6 and 3.7 (PR #46) and bumping of scipy version to >=1.8.0 is a big step. I would highly recommend doing a release at this point.

For example, I tried building on Python 3.10.2 with scipy==1.8.0 locally (Ubuntu 20.04) with pyproximal==0.2.0, but it failed since the scipy interface changed. Specifically it was line 8 in pyproximal/proximal/VStack.py that failed:

from scipy.sparse.linalg.interface import _get_dtype

Furthermore, when testing my repository in a Docker container (specifically python:3.9.10-slim) I was for this reason forced to use the main branch; however, these standard docker containers do not ship with git installed, and therefore I had to customize the container to first install it before being able to run it.

Maybe not a big deal, but if everything works of the shelf then people are encouraged to use the library more.

Once Pylops v1.18.1 has been released, add the following test to test_norms:

@pytest.mark.parametrize("par", [(par1), (par2)])
def test_L2_dense(par):
    """L2 norm of Op*x with dense Op and proximal/dual proximal
    """
    b = np.zeros(par['nx'], dtype=par['dtype'])
    d = np.random.normal(0., 1., par['nx']).astype(par['dtype'])
    l2 = L2(Op=MatrixMult(np.diag(d), dtype=par['dtype']),
            b=b, sigma=par['sigma'], densesolver='numpy')

    # norm
    x = np.random.normal(0., 1., par['nx']).astype(par['dtype'])
    assert l2(x) == (par['sigma'] / 2.) * np.linalg.norm(d * x) ** 2

    # prox: since Op is a Diagonal operator the denominator becomes
    # 1 + sigma*tau*d[i] for every i
    tau = 2.
    den = 1. + par['sigma'] * tau * d ** 2
    assert_array_almost_equal(l2.prox(x, tau), x / den, decimal=4)

This needs checked:

See https://regularize.wordpress.com/2017/02/24/moreaus-identity-for-monotone-operators/amp/