diff --git a/deerlab/snlls.py b/deerlab/snlls.py
index 8e018f212..f15246f43 100644
--- a/deerlab/snlls.py
+++ b/deerlab/snlls.py
@@ -293,6 +293,9 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx
         elif nnlsSolver == 'cvx':
             linSolver = lambda AtA, Aty: cvxnnls(AtA, Aty, tol=lin_tol, maxiter=lin_maxiter)
         parseResult = lambda result: result
+    
+    # Ensure correct formatting and shield against float-point errors
+    validateResult = lambda result: np.maximum(lbl,np.minimum(ubl,np.atleast_1d(result)))
     # ----------------------------------------------------------
 
     # Containers for alpha-update checks
@@ -319,9 +322,9 @@ def linear_problem(A,optimize_alpha,alpha):
          
         # Solve the linear least-squares problem
         result = linSolver(AtA, Aty)
-        linfit = parseResult(result)
-        linfit = np.atleast_1d(linfit)
-        
+        result = parseResult(result)
+        linfit = validateResult(result)
+
         return linfit, alpha
     #===========================================================================
 
@@ -566,4 +569,3 @@ def _plot(subsets,y,yfit,show):
         plt.close()
     return fig
 # ===========================================================================================
-
diff --git a/deerlab/utils/utils.py b/deerlab/utils/utils.py
index b9272b6e8..d61476bc3 100644
--- a/deerlab/utils/utils.py
+++ b/deerlab/utils/utils.py
@@ -409,7 +409,7 @@ def Jacobian(fcn, x0, lb, ub):
 
     """
     J = opt._numdiff.approx_derivative(fcn,x0,method='2-point',bounds=(lb,ub))
-    J = np.atleast_2d(J)
+    J = np.atleast_2d(J) 
     return J
 #===============================================================================