Repository: gbrammer/eazy-py

Issue:
Eazy returns "LinAlgError: Singular matrix"
I am running eazy on a mock catalogue and it returns the error "LinAlgError: Singular matrix", when running standard_output.  This occurs *both* with and without `get_err = False`.  No errors occur in earlier fittings for example in `fit_catalog`.
The error appears to be located in the function `safe_invert` from utils.py. The function appears to fail when the array has values close to 1. 
For example the matrix that it breaks on is:

`arr = [[1.,         0.99998837, 0.99992621, 0.99999995],
 [0.99998837, 1.,         0.9999728,  0.99998979],
 [0.99992621, 0.9999728,  1.  ,       0.99992987],
 [0.99999995, 0.99998979, 0.99992987, 1.        ]]`.

The photometry of the source which causes this issue is:
```
fnu = [ 2.02334821e-02  2.94720661e-02  2.23611714e-05  9.08530176e-01
  1.86806238e+00  1.81287825e+00 -2.84472341e-03 -1.93866696e-02
  7.91999046e-03  2.16457676e-02  1.47525286e-02  2.76298486e-02]

efnu = [1.04980990e-02 6.03997521e-02 5.02377301e-02 5.58287837e-02
 9.90000000e+01 9.90000000e+01 5.50845750e-02 1.99999996e-02
 2.00001560e-02 3.16986032e-02 5.02379462e-02 1.82402372e-01]
```


![Image](https://github.com/user-attachments/assets/5531723b-070f-4fb3-bcf8-0a301f79801a)

Hints from issue comments:
@nallen98: Translate file:
```
f_VIS F410
e_VIS E410
f_Y F411
e_Y E411
f_J F412
e_J E412
f_H F413
e_H E413
f_irac1 F18
e_irac1 E18
f_irac2 F19
e_irac2 E19
f_uCFHT F72
e_uCFHT E72
f_gHSC F314
e_gHSC E314
f_rHSC F315
e_rHSC E315
f_iHSC F316
e_iHSC E316
f_zHSC F317
e_zHSC E317
f_uLSST F419
e_uLSST E419
f_gLSST F420
e_gLSST E420
f_rLSST F421
e_rLSST E421
f_iLSST F422
e_iLSST E422
f_zLSST F423
e_zLSST E423
f_yLSST F424
e_yLSST E424
```@nallen98: ```python
params = {}
params['CATALOG_FILE'] = cat_fname

params['MAIN_OUTPUT_FILE'] = f'/data/natalie/Euclid_Cats/Eazy_v2025_Cats/v_{version}/Euclid_{field}_Uniform_v2025.eazypy'
params['MW_EBV'] = 0.0
params['CAT_HAS_EXTCORR'] = True
params['Z_STEP'] = 0.01
params['Z_MIN'] = 0.0
params['Z_MAX'] = 15
params['APPLY_PRIOR'] = False
params['PRIOR_ABZP'] = 25 
params['TEMPLATES_FILE'] = 'templates/sfhz/blue_sfhz_13.param'
params['FIX_ZSPEC'] = False
params['IGM_SCALE_TAU'] = 1.0
params['FILTERS_RES'] = 'FILTER.RES.latest'


self = eazy.photoz.PhotoZ(param_file=None, translate_file="zphot.translate", zeropoint_file=None, 
                          params=params, load_prior=True, load_products=False)

self.fit_catalog(self.idx, n_proc=60)

self.fit_phoenix_stars(sonora=True, stars=None, just_dwarfs=False)

zout, hdu = self.standard_output(
    simple=False, 
    rf_pad_width=0.5,
    rf_max_err=2, 
    prior=True,
    beta_prior=True, 
    absmag_filters=[], 
    extra_rf_filters=[],
    n_proc=0,
    par_skip = len(self.idx)
)

```@nallen98: New file:

[Test_cdfs_objid-108709.csv](https://github.com/user-attachments/files/18869003/Test_cdfs_objid-108709.csv)@gbrammer: Thanks for sending this.  I can't reproduce the error using the provided catalog, which could be due to differences in machine architecture and floating point precision.  It's definitely worrisome that all of the values in the ``arr`` matrix differ from unity only by something similar to machine epsilon, so it may or may not be singular depending on your machine, the c compiler used to build numpy, etc.

In any case, I pushed an update in the PR mentioned above that should resolve this and that is now merged into release version `0.8.4`.  Can you try updating and checking if the problem persists?

Can you also please file an Issue on [eazy-photoz](https://github.com/gbrammer/eazy-photoz/issues) to downsample the LSST filter bandpasses?  They seem to cause the `PhotoZ` initialization to take forever, I think unnecessarily.@gbrammer: @nallen98 , can you please close this issue if the fix above solved it?  Thanks!

Solution patch:
diff --git a/eazy/photoz.py b/eazy/photoz.py
index ab97c33..cd0f809 100644
--- a/eazy/photoz.py
+++ b/eazy/photoz.py
@@ -6542,9 +6542,14 @@ def template_lsq(fnu_i, efnu_i, Ain, TEFz, zp, ndraws, fitter, renorm_t, hess_th
         #try:
         if ok_temp.sum() > 0:
             covar = utils.safe_invert(mat)
-            draws = np.random.multivariate_normal((coeffs_i * An)[ok_temp], 
-                                                  covar, 
-                                                  size=ndraws)
+            if np.isfinite(covar).sum() != covar.size:
+                draws = np.full((ndraws, ok_temp.sum()), np.nan)
+            else:
+                draws = np.random.multivariate_normal(
+                    (coeffs_i * An)[ok_temp],
+                    covar,
+                    size=ndraws
+                )
             coeffs_draw[:, ok_temp] = draws / An[ok_temp]
         else:
             #print('Error getting coeffs draws')
diff --git a/eazy/utils.py b/eazy/utils.py
index 2e1d36e..86802b0 100644
--- a/eazy/utils.py
+++ b/eazy/utils.py
@@ -369,14 +369,37 @@ def fill_between_steps(x, y, z, ax=None, *args, **kwargs):
     ax.fill_between(xfull[so], yfull[so], zfull[so], *args, **kwargs)
 
 
-def safe_invert(arr):
+def safe_invert(arr, tol=None):
     """
     Version-safe matrix inversion using `numpy.linalg.inv` or `numpy.matrix.I`
+
+    Parameters
+    ----------
+    arr : array-like
+        Square matrix to invert
+
+    tol : float
+        Tolerance parameter for calculating if arr is singular using the rank from
+        ``np.linalg.matrix_rank(arr, tol=tol)``.  If None, use the numpy default
+        derived from the matrix size and machine epsilon.  If ``arr`` is found to be
+        singular (rank < max(arr.shape)), don't try to invert it and just return an
+        array of NaN values.
+
+    Returns
+    -------
+    inv : array-like
+        Inverse of ``arr``
+
     """
+    # Check that matrix is not singular
+    rank = np.linalg.matrix_rank(arr, tol=tol)
+    if rank < np.max(arr.shape):
+        return arr * np.nan
+
     try:
         from numpy.linalg import inv
         _inv = inv(arr)
-    except:
+    except ImportError:
         _inv = np.matrix(arr).I.A
     
     return _inv


Test patch:
diff --git a/eazy/tests/test_utils.py b/eazy/tests/test_utils.py
index a14504a..879b340 100644
--- a/eazy/tests/test_utils.py
+++ b/eazy/tests/test_utils.py
@@ -92,6 +92,11 @@ def test_invert():
     ainv = utils.safe_invert(a)
     assert(np.allclose(np.dot(a, ainv), np.eye(2)))
 
+    # Singular matrix returns np.nan inverse without raising exception
+    a = np.ones((2,2))
+    ainv = utils.safe_invert(a)
+    assert np.all(~np.isfinite(utils.safe_invert(a)))
+
 
 def test_query_string():
     """