LinAlgError when setting spectra

[b-biswas]

Hello,

I am trying to initialize sources with scarlet.initialization.init_all_sources but I get a numpy.linalg.LinAlgError from this line of the function scarlet.initialization.set_spectra_to_match.

The code I used for this run is uploaded here
and I get the following output:

File "run_scarlet.py", line 195, in <module>
    scarlet_current_predictions = predict_with_scarlet(
  File "run_scarlet.py", line 102, in predict_with_scarlet
    sources, skipped = scarlet.initialization.init_all_sources(
  File "/pbs/throng/lsst/users/bbiswas/miniconda3/envs/madness/lib/python3.8/site-packages/scarlet-1.0.1+g7eb3fe0-py3.8-linux-x86_64.egg/scarlet/initialization.py", line 399, in init_all_sources
    set_spectra_to_match(sources, observations)
  File "/pbs/throng/lsst/users/bbiswas/miniconda3/envs/madness/lib/python3.8/site-packages/scarlet-1.0.1+g7eb3fe0-py3.8-linux-x86_64.egg/scarlet/initialization.py", line 593, in set_spectra_to_match
    covar = np.linalg.inv(mw @ m.T)
  File "<__array_function__ internals>", line 180, in inv
  File "/pbs/throng/lsst/users/bbiswas/miniconda3/envs/madness/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 552, in inv
    ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
  File "/pbs/throng/lsst/users/bbiswas/miniconda3/envs/madness/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 89, in _raise_linalgerror_singular
    raise LinAlgError("Singular matrix")
numpy.linalg.LinAlgError: Singular matrix

Here are a few things I tried to resolve this issue:

• changing the min_snr.
• model only as a single ExtendedSource.
• set_spectra = False

As expected set_spectra = False manages to avoid this error but quite surprisingly the error also disappears when setting shifting=True although this parameter should not have an impact since the sources are well within the boundary of the images.

I am currently using the scarlet btk-v1 release and I noticed in the BTK example runs of scarlet, a try catch block was introduced to catch this exception.

Does this error come from extremely faint sources?

[pmelchior]

Thanks for reporting. The code that fails here is a linear solver to get the best-fit amplitude across all bands, aka the SEDs, for all sources given their initial morphology. In absence of zero or nan weights, there's two ways of getting that:

1. The source is so faint that there's no pixel above the noise level, so the morphology is zero.
2. The morphologies are identical for at least two sources. That happens when the center of two sources lies in the same pixel. Then the morphology initialization yields the same result for both.

In either case the best-fit solution is ill-defined and numerically unstable. Hence, the exception when trying to invert the covariance matrix of the best-fit.

It's interesting to see the btk had to catch this failure, so this isn't completely rare in simulated cases, whereas we have not seen it to be overly problematic with actual detection catalogs, where it's essentially impossible to have two detections so close together. I would say that the solutions of a deblender in such cases will be pretty useless (and that's probably true for any deblender: if the data can't tell the difference, you get either the initial guess or the prior), but exceptions should be avoided or caught. So, I suggest that we fix this problem on our end.

The fix will be to test for both of the cases listed above:

1. If the morphology is zero, it should be replaced by a point source image (this should happen anyway, but we might need to check that it really does).
2. If two or more morphologies are identical (including their position in the frame), we should merge them into one for the linear SED solver. Once we have the SED for the merged morphology, we can split it again into two distinct SEDs, while keeping the morphologies the same. A 50-50 split will do nothing to break the degeneracies, so we should probably tilt one component to the red and the other to the blue, so that each can develop differently during the fit.

Something like

sed1 = sed * np.arange(1, C+1)/C
sed2 = sed * (1 - np.arange(1, C+1)/C)

should work.

[b-biswas]

I took a closer look at the cases where the exception occurs and I confirm that the exception occurs in the fields that have two galaxies too close to one another (centers lie in the same pixel).

[pmelchior]

Thanks for confirming. We should be able to patch that soon