bSDMM algorithm is broken when sources resize

[fred3m]

In the bSDMM algorithm, the L matrix is copied, not referenced in https://github.com/pmelchior/proxmin/blob/master/proxmin/algorithms.py#L480-L486. So for algorithms that resize X and Ls, the new _L is not resized in proxmin, and results in a dimensional error.

proxmin should allow a user to pass a list of matix adapters to L, preventing this behavior.

[pmelchior]

I don't think the explanation for this bug is correct. bSDMM only holds a reference to the matrix through MatrixAdapter, which itself only holds a reference: https://github.com/pmelchior/proxmin/blob/master/proxmin/utils.py#L35

I could see that we get the spectral norm wrong after resizing, but you report a dimensionality problem, which I don't understand.

[fred3m]

Try it:

import numpy as np
import scarlet
import scarlet.constraint as sc

# Load the sample images
#data = np.load("../data/test_sim/data.npz")
data = np.load("/Users/fred/projects/scarlet/data/test_sim/data.npz")
images = data["images"]

# Load the thruth detection catalog
from astropy.table import Table as ApTable
catalog = ApTable.read("/Users/fred/projects/scarlet/data/test_sim/true_catalog.fits")
bg_rms = np.array([20]*len(images))

# Initialize a constraint that is used for all of the sources
constraints = (
    sc.SimpleConstraint(), # sed sum to unity, all elements of SED and morph are non-negative
    sc.MonotonicityConstraint(use_nearest=False), # prox_g monotonicity
    sc.DirectSymmetryConstraint() # prox_f perfect symmetry
)

# Initialize the Sources
sources = [scarlet.source.PointSource(
    (src['y'],src['x']), # center coordinates in `images`
    images, # data cube (bands, Ny, Nx)
    (15,15), # initial shape of the bounding box
    constraints=constraints,
) for src in catalog]

blend = scarlet.Blend(sources).set_data(images, bg_rms=bg_rms)
blend.fit(11)

output:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-14-bd13a15d594b> in <module>
     36 
     37 blend = scarlet.Blend(sources).set_data(images, bg_rms=bg_rms)
---> 38 blend.fit(11)

~/projects/scarlet/scarlet/blend.py in fit(self, steps, e_rel)
    187             res = proxmin.algorithms.bsdmm(X, self._prox_f, self._steps_f, proxs_g, steps_g=steps_g,
    188                 Ls=self._Ls, update_order=update_order, steps_g_update=steps_g_update, max_iter=steps,
--> 189                 e_rel=self._e_rel, e_abs=self._e_abs)
    190 
    191         converged, errors = res

~/projects/proxmin/proxmin/algorithms.py in bsdmm(X, proxs_f, steps_f_cb, proxs_g, steps_g, Ls, update_order, steps_g_update, max_iter, e_rel, e_abs, traceback)
    528 
    529             # update the variables
--> 530             LX[j], R[j], S[j] = utils.update_variables(X[j], Z[j], U[j], proxs_f_j, steps_f[j], proxs_g[j], steps_g_[j], _L[j])
    531 
    532             # convergence criteria, adapted from Boyd 2011, Sec 3.3.1

~/projects/proxmin/proxmin/utils.py in update_variables(X, Z, U, prox_f, step_f, prox_g, step_g, L)
    347     else:
    348         M = len(prox_g)
--> 349         dX = np.sum([step_f/step_g[i] * L[i].T.dot(L[i].dot(X) - Z[i] + U[i]) for i in range(M)], axis=0)
    350         X[:] = prox_f(X - dX, step_f)
    351         LX = [None] * M

~/projects/proxmin/proxmin/utils.py in <listcomp>(.0)
    347     else:
    348         M = len(prox_g)
--> 349         dX = np.sum([step_f/step_g[i] * L[i].T.dot(L[i].dot(X) - Z[i] + U[i]) for i in range(M)], axis=0)
    350         X[:] = prox_f(X - dX, step_f)
    351         LX = [None] * M

~/projects/proxmin/proxmin/utils.py in dot(self, X)
     66         # dot product
     67         if self.axis == 1:
---> 68             return self.L.dot(X.reshape(-1)).reshape(X.shape[0], -1)
     69         raise NotImplementedError("MatrixAdapter.dot() is not useful with axis=0.\n"
     70                                   "Use regular matrix dot product instead!")

~/miniconda3/envs/scarlet/lib/python3.7/site-packages/scipy/sparse/base.py in dot(self, other)
    362 
    363         """
--> 364         return self * other
    365 
    366     def power(self, n, dtype=None):

~/miniconda3/envs/scarlet/lib/python3.7/site-packages/scipy/sparse/base.py in __mul__(self, other)
    498             # dense row or column vector
    499             if other.shape != (N,) and other.shape != (N, 1):
--> 500                 raise ValueError('dimension mismatch')
    501 
    502             result = self._mul_vector(np.ravel(other))

ValueError: dimension mismatch

[fred3m]

This is what broke the docs the day that you merged the resize branch in September. I tried it out and the error occurs when some of the sources are resized from (15,15) to (25, 25).

[fred3m]

I just double checked, and blend._Ls has shape 625x625 for a source that was resized, but _L in promin has shape 125 (15x15). So proxmin doesn't see the new _L. I think this is because _Ls is a dynamic property of the Blend class but when proxmin is initialized it points to the L matrix that it was initialized with.

[pmelchior]

Good point. It is very well possible, and now even likely, that I didn't test resizing with a L-type constraint. I wonder if we should keep dragging L-type constraints around in scarlet, which is the only application where we need dynamic resizing. Do we still use any in scarlet?

[fred3m]

Isn't Remy using it? I don't think that we use it anymore but I'm not sure if anyone else does. Maybe we should start a mailing list so that we have a better way to keep track of users (outside of LSST/DESC) and solicit feedback when we are implementing a breaking change.

[pmelchior]

If we want to use wavelet-sparsity, then yes, otherwise no.

I'm conflicted about the mailing list, I don't think they work well to address the underlying issue. But I'd still be willing to move fast and break things at this point.

Two alternatives:

• No L matrices in scarlet (if resizing is used)
• A mechanism to get the resized matrix

I prefer the latter, and it should be reasonably straightforward to solve it in the same way we resize some other operators (like your monotonicity prox).

In either case, I opt to solve this issue on the scarlet side, not in proxmin.

[fred3m]

I agree that we should use the later option and it should be straight forward to implement in the scarlet side.