diff --git a/python/lsst/multiprofit/fitting/fit_psf.py b/python/lsst/multiprofit/fitting/fit_psf.py
index fae875d..feab73b 100644
--- a/python/lsst/multiprofit/fitting/fit_psf.py
+++ b/python/lsst/multiprofit/fitting/fit_psf.py
@@ -29,6 +29,7 @@
 from abc import abstractmethod
 from functools import cached_property
 import logging
+import math
 import time
 from typing import Any, ClassVar, Mapping, Type
 
@@ -83,6 +84,11 @@ class CatalogPsfFitterConfig(CatalogFitterConfig):
         doc="PSF model configuration",
     )
     prior_axrat_mean = pexConfig.Field[float](default=0.95, doc="Mean for axis ratio prior")
+    sigma_min = pexConfig.Field[float](
+        default=0.8,
+        doc="Minimum sigma in pixels for PSF components. Must be >=0.8 to avoid undersampling.",
+        check=lambda x: x >= 0.8,
+    )
 
     def make_psf_model(
         self,
@@ -565,8 +571,13 @@ def fit(
             for param in get_params_uniq(model_source, linear=False, channel=g2f.Channel.NONE, fixed=False)
             if isinstance(param, g2f.ProperFractionParameterD)
         )
-        # We're fitting the PSF, so there's nothing to convolve with
+        # We're fitting the PSF, so make a single Gaussian
         model_psf = make_psf_model_null()
+        # Set the size to the minimum sigma to avoid undersampling
+        ellipse = model_psf.components[0].ellipse
+        ellipse.sigma_x_param.value = config.sigma_min
+        ellipse.sigma_y_param.value = config.sigma_min
+        sigma_min_sq = config.sigma_min**2
 
         catalog = catexp.get_catalog()
         n_rows = len(catalog)
@@ -574,7 +585,12 @@ def fit(
 
         results, columns = config.make_catalog(n_rows)
         prefix = config.prefix_column
-        columns_param = {f"{prefix}{key}": param for key, param in params.items()}
+        columns_param = {}
+        for key, param in params.items():
+            is_sigma = isinstance(param, g2f.SigmaXParameterD) or isinstance(param, g2f.SigmaYParameterD)
+            columns_param[f"{prefix}{key}"] = param, is_sigma
+            if is_sigma:
+                param.value = math.sqrt(max(param.value**2 - sigma_min_sq, 0.1))
 
         # dummy size for first iteration
         size, size_new = 0, 0
@@ -631,9 +647,12 @@ def fit(
                 if config.config_fit.eval_residual:
                     results[f"{prefix}n_eval_jac"][idx] = result_full.n_eval_jac
 
-                for (key, param), value in zip(columns_param.items(), result_full.params_best):
+                for (key, (param, is_sigma)), value in zip(columns_param.items(), result_full.params_best):
                     param.value_transformed = value
-                    results[key][idx] = param.value
+                    value = param.value
+                    if is_sigma:
+                        value = math.sqrt(sigma_min_sq**2 + value**2)
+                    results[key][idx] = value
 
                 time_final = time.process_time()
                 results[f"{prefix}time_full"][idx] = time_final - time_init
diff --git a/python/lsst/multiprofit/fitting/fit_source.py b/python/lsst/multiprofit/fitting/fit_source.py
index bf6a35f..21bc531 100644
--- a/python/lsst/multiprofit/fitting/fit_source.py
+++ b/python/lsst/multiprofit/fitting/fit_source.py
@@ -908,6 +908,9 @@ def fit(
                     _ = self.modeller.fit_model_linear(model_psf)
                     model_psf.setup_evaluators(evaluatormode=g2f.EvaluatorMode.loglike)
                     loglike_psfmodel = model_psf.evaluate()
+                    # Reset fluxes for the next fit
+                    for param in fluxes_psmodel.values():
+                        param.value = 1.
                     results[f"{prefix}delta_lnL_fit_ps"][idx] = loglike_final[0] - loglike_psfmodel[0]
 
                 if compute_errors:
diff --git a/tests/test_fit_bootstrap_model.py b/tests/test_fit_bootstrap_model.py
index 7c9b7b1..fb6c9d1 100644
--- a/tests/test_fit_bootstrap_model.py
+++ b/tests/test_fit_bootstrap_model.py
@@ -19,6 +19,8 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
+import math
+
 import astropy.table
 import lsst.gauss2d.fit as g2f
 from lsst.multiprofit.componentconfig import (
@@ -223,11 +225,13 @@ def test_fit_psf(config_fitter_psfs, tables_psf_fits):
         params_init = psf_model_init.parameters()
         params_fit = psf_model_fit.parameters()
         assert len(params_init) == len(params_fit)
+        sigma_min_sq = config_data_psf.config.sigma_min**2
         for p_init, p_meas in zip(params_init, params_fit):
             assert p_meas.fixed == p_init.fixed
             if p_meas.fixed:
                 assert p_init.value == p_meas.value
             else:
+                value = p_meas.value
                 # TODO: come up with better (noise-dependent) thresholds here
                 if isinstance(p_init, g2f.IntegralParameterD):
                     atol, rtol = 0, 0.02
@@ -235,9 +239,11 @@ def test_fit_psf(config_fitter_psfs, tables_psf_fits):
                     atol, rtol = 0.1, 0.01
                 elif isinstance(p_init, g2f.RhoParameterD):
                     atol, rtol = 0.05, 0.1
+                elif isinstance(p_init, g2f.SigmaXParameterD) or isinstance(p_init, g2f.SigmaYParameterD):
+                    value = math.sqrt(value**2 + sigma_min_sq)
                 else:
                     atol, rtol = 0.01, 0.1
-                assert np.isclose(p_init.value, p_meas.value, atol=atol, rtol=rtol)
+                assert np.isclose(p_init.value, value, atol=atol, rtol=rtol)
 
 
 def test_fit_source(config_fitter_source, config_data_sources):