Tutorial: Non-Rigid Structure from Motion

[marcusvaltonen]

I thought this tutorial could expand the repertoire of pyproximal, with some fun and visual applications from computer vision.

In the future, when time permits, I plan to implement the bilinear methods #50 which fits well in this context (and I could expand the tutorial then as well).

[marcusvaltonen]

I did try to wrap the block diagonal operations for Frobenius norm proximal operator, but it was very slow. Perhaps you could shine some light on that?

import pylops

from pyproximal.proximal.L2 import L2


class BlockDiagFrobenius(L2):
    r"""Block-diagonal Frobenius wrapped."""
    def __init__(self, dim, blkdiag, B=None, sigma=1., niter=10, x0=None, warm=True,
                 densesolver=None):
        self.dim = dim
        R = pylops.BlockDiag(blkdiag)
        Op = pylops.Kronecker(pylops.Identity(dim[1]), R)
        if x0:
            x0 = x0.T.ravel()
        super().__init__(Op, B.T.ravel(), None, sigma, None, None, niter, x0, warm, densesolver)

    def prox(self, x, tau):
        # Make sure that the array is column-wise stacked
        xx = x.reshape(self.dim).T.ravel()
        y = super().prox(xx, tau)
        return y.reshape(self.dim, order='F').ravel()

[mrava87]

Fantastic!

Let me merge this and I can take a look at your question in coming days.

Usually Knonecker is what I call the last resort, let me explain why. When you have a method that works say on 1D arrays and you want to apply it to all columns of a 2d array, Kronecker allows you to do that but all it really does under the hood is a for loop https://github.com/PyLops/pylops/blob/31048ec1e05a1819b53ed59144d23c673cf69b7f/pylops/basicoperators/Kronecker.py#L72 (see how matmat is implemented). This is nice for quick testing to see if something makes sense before embarking on more coding, but the efficient approach would be to take the operator you are interested and make it broacast in nd. Most of our operators actually do that (see for example Diagonal).

One intermediate step is to enable multiprocess in BlockDiag operator with nproc. That should speed things up with a few caveats, the heavy stuff should be in the operators that make up the blockdiag whislt the dim[1] should be rather small. This is easy to try as it is all ready for you. Another thing which may make sense in the opposite scenario is to add multiprocess functionalities to matmat, I never thought about it simply because I never had a use case but one is here so I have no problem against adding this functionality.

Disclaimer: we are very close to releasing PyLops2. I adds a lot of new interesting features and improves some early poor design (see PyLops/pylops#438). This means of course some breaking changes, and the next step is to make pyproximal PyLops2 compliant. This also means we dont accept new features in PyLops1 branch (only bug fixes).