diff --git a/.github/workflows/commit_message.yml b/.github/workflows/commit_message.yml
index 5a6666b84b2b..05a3ad9dc02f 100644
--- a/.github/workflows/commit_message.yml
+++ b/.github/workflows/commit_message.yml
@@ -18,6 +18,7 @@ jobs:
       message: ${{ steps.skip_check.outputs.message }}
     steps:
       - name: Checkout scipy
+        if: ${{ !contains(github.actor, 'nektos/act') }}
         uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2
         # Gets the correct commit message for pull request
         with:
@@ -29,10 +30,14 @@ jobs:
         # and changing the output here based on the combination of tags (if desired).
         run: |
           set -xe
-          COMMIT_MSG=$(git log --no-merges -1)
           RUN="1"
-          if [[ "$COMMIT_MSG" == *"[lint only]"* || "$COMMIT_MSG" == *"[docs only]"* ]]; then
-              RUN="0"
+          # For running a job locally with Act (https://github.com/nektos/act),
+          # always run the job rather than skip based on commit message
+          if [[ $ACT != true ]]; then
+            COMMIT_MSG=$(git log --no-merges -1)
+            if [[ "$COMMIT_MSG" == *"[lint only]"* || "$COMMIT_MSG" == *"[docs only]"* ]]; then
+                RUN="0"
+            fi
           fi
           echo "message=$RUN" >> $GITHUB_OUTPUT
           echo github.ref ${{ github.ref }}
diff --git a/scipy/cluster/vq.py b/scipy/cluster/vq.py
index a791e2956070..34045c1357fe 100644
--- a/scipy/cluster/vq.py
+++ b/scipy/cluster/vq.py
@@ -310,20 +310,18 @@ def _kmeans(obs, guess, thresh=1e-5, xp=None):
     code_book = guess
     diff = xp.inf
     prev_avg_dists = deque([diff], maxlen=2)
+
+    np_obs = np.asarray(obs)
     while diff > thresh:
         # compute membership and distances between obs and code_book
         obs_code, distort = vq(obs, code_book, check_finite=False)
         prev_avg_dists.append(xp.mean(distort, axis=-1))
         # recalc code_book as centroids of associated obs
-        obs = np.asarray(obs)
         obs_code = np.asarray(obs_code)
-        code_book, has_members = _vq.update_cluster_means(obs, obs_code,
+        code_book, has_members = _vq.update_cluster_means(np_obs, obs_code,
                                                           code_book.shape[0])
-        obs = xp.asarray(obs)
-        obs_code = xp.asarray(obs_code)
-        code_book = xp.asarray(code_book)
-        has_members = xp.asarray(has_members)
         code_book = code_book[has_members]
+        code_book = xp.asarray(code_book)
         diff = xp.abs(prev_avg_dists[0] - prev_avg_dists[1])
 
     return code_book, prev_avg_dists[1]
@@ -814,7 +812,7 @@ def kmeans2(data, k, iter=10, thresh=1e-5, minit='random',
 
     data = np.asarray(data)
     code_book = np.asarray(code_book)
-    for i in range(iter):
+    for _ in range(iter):
         # Compute the nearest neighbor for each obs using the current code book
         label = vq(data, code_book, check_finite=check_finite)[0]
         # Update the code book by computing centroids
diff --git a/scipy/differentiate/_differentiate.py b/scipy/differentiate/_differentiate.py
index 2b50b1732dca..0e104a071055 100644
--- a/scipy/differentiate/_differentiate.py
+++ b/scipy/differentiate/_differentiate.py
@@ -802,27 +802,47 @@ def jacobian(f, x, *, tolerances=None, maxiter=10, order=8, initial_step=0.5,
     Notes
     -----
     Suppose we wish to evaluate the Jacobian of a function
-    :math:`f: \mathbf{R}^m \rightarrow \mathbf{R}^n`, and assign to variables
+    :math:`f: \mathbf{R}^m \rightarrow \mathbf{R}^n`. Assign to variables
     ``m`` and ``n`` the positive integer values of :math:`m` and :math:`n`,
-    respectively. If we wish to evaluate the Jacobian at a single point,
-    then:
+    respectively, and let ``...`` represent an arbitrary tuple of integers.
+    If we wish to evaluate the Jacobian at a single point, then:
 
     - argument `x` must be an array of shape ``(m,)``
-    - argument `f` must be vectorized to accept an array of shape ``(m, p)``.
-      The first axis represents the :math:`m` inputs of :math:`f`; the second
-      is for evaluating the function at multiple points in a single call.
-    - argument `f` must return an array of shape ``(n, p)``. The first
-      axis represents the :math:`n` outputs of :math:`f`; the second
-      is for the result of evaluating the function at multiple points.
+    - argument `f` must be vectorized to accept an array of shape ``(m, ...)``.
+      The first axis represents the :math:`m` inputs of :math:`f`; the remainder
+      are for evaluating the function at multiple points in a single call.
+    - argument `f` must return an array of shape ``(n, ...)``. The first
+      axis represents the :math:`n` outputs of :math:`f`; the remainder
+      are for the result of evaluating the function at multiple points.
     - attribute ``df`` of the result object will be an array of shape ``(n, m)``,
       the Jacobian.
 
     This function is also vectorized in the sense that the Jacobian can be
     evaluated at ``k`` points in a single call. In this case, `x` would be an
     array of shape ``(m, k)``, `f` would accept an array of shape
-    ``(m, k, p)`` and return an array of shape ``(n, k, p)``, and the ``df``
+    ``(m, k, ...)`` and return an array of shape ``(n, k, ...)``, and the ``df``
     attribute of the result would have shape ``(n, m, k)``.
 
+    Suppose the desired callable ``f_not_vectorized`` is not vectorized; it can
+    only accept an array of shape ``(m,)``. A simple solution to satisfy the required
+    interface is to wrap ``f_not_vectorized`` as follows::
+
+        def f(x):
+            return np.apply_along_axis(f_not_vectorized, axis=0, arr=x)
+
+    Alternatively, suppose the desired callable ``f_vec_q`` is vectorized, but
+    only for 2-D arrays of shape ``(m, q)``. To satisfy the required interface,
+    consider::
+
+        def f(x):
+            m, batch = x.shape[0], x.shape[1:]  # x.shape is (m, ...)
+            x = np.reshape(x, (m, -1))  # `-1` is short for q = prod(batch)
+            res = f_vec_q(x)  # pass shape (m, q) to function
+            n = res.shape[0]
+            return np.reshape(res, (n,) + batch)  # return shape (n, ...)
+
+    Then pass the wrapped callable ``f`` as the first argument of `jacobian`.
+
     References
     ----------
     .. [1] Jacobian matrix and determinant, *Wikipedia*,
diff --git a/scipy/linalg/tests/test_lapack.py b/scipy/linalg/tests/test_lapack.py
index f32ce213a5aa..3e4294602f03 100644
--- a/scipy/linalg/tests/test_lapack.py
+++ b/scipy/linalg/tests/test_lapack.py
@@ -3479,9 +3479,9 @@ def test_sy_hetrs(mtype, dtype, lower):
     names = f'{mtype}trf', f'{mtype}trf_lwork', f'{mtype}trs'
     trf, trf_lwork, trs = get_lapack_funcs(names, dtype=dtype)
     lwork = trf_lwork(n, lower=lower)
-    ldu, ipiv, info = trf(A, lwork=lwork)
+    ldu, ipiv, info = trf(A, lwork=lwork, lower=lower)
     assert info == 0
-    x, info = trs(a=ldu, ipiv=ipiv, b=b)
+    x, info = trs(a=ldu, ipiv=ipiv, b=b, lower=lower)
     assert info == 0
     eps = np.finfo(dtype).eps
     assert_allclose(A@x, b, atol=100*n*eps)
diff --git a/scipy/optimize/_lsq/dogbox.py b/scipy/optimize/_lsq/dogbox.py
index 6bb5abbe7902..b986929626f2 100644
--- a/scipy/optimize/_lsq/dogbox.py
+++ b/scipy/optimize/_lsq/dogbox.py
@@ -45,6 +45,8 @@
 
 from scipy.sparse.linalg import LinearOperator, aslinearoperator, lsmr
 from scipy.optimize import OptimizeResult
+from scipy._lib._util import _call_callback_maybe_halt
+
 
 from .common import (
     step_size_to_bound, in_bounds, update_tr_radius, evaluate_quadratic,
@@ -147,7 +149,7 @@ def dogleg_step(x, newton_step, g, a, b, tr_bounds, lb, ub):
 
 
 def dogbox(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale,
-           loss_function, tr_solver, tr_options, verbose):
+           loss_function, tr_solver, tr_options, verbose, callback=None):
     f = f0
     f_true = f.copy()
     nfev = 1
@@ -322,6 +324,18 @@ def dogbox(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale,
             actual_reduction = 0
 
         iteration += 1
+        
+        # Call callback function and possibly stop optimization
+        if callback is not None:
+            intermediate_result = OptimizeResult(
+                x=x, fun=f, nit=iteration, nfev=nfev)
+            intermediate_result["cost"] = cost_new
+
+            if _call_callback_maybe_halt(
+                callback, intermediate_result
+            ):
+                termination_status = -2
+                break
 
     if termination_status is None:
         termination_status = 0
diff --git a/scipy/optimize/_lsq/least_squares.py b/scipy/optimize/_lsq/least_squares.py
index e01ce04a8992..71fbba41762e 100644
--- a/scipy/optimize/_lsq/least_squares.py
+++ b/scipy/optimize/_lsq/least_squares.py
@@ -14,8 +14,11 @@
 from .dogbox import dogbox
 from .common import EPS, in_bounds, make_strictly_feasible
 
+    
+from scipy.optimize._optimize import _wrap_callback
 
 TERMINATION_MESSAGES = {
+    -2: "Stopped because `callback` function raised `StopIteration` or returned `True`",
     -1: "Improper input parameters status returned from `leastsq`",
     0: "The maximum number of function evaluations is exceeded.",
     1: "`gtol` termination condition is satisfied.",
@@ -242,7 +245,9 @@ def least_squares(
         fun, x0, jac='2-point', bounds=(-np.inf, np.inf), method='trf',
         ftol=1e-8, xtol=1e-8, gtol=1e-8, x_scale=1.0, loss='linear',
         f_scale=1.0, diff_step=None, tr_solver=None, tr_options=None,
-        jac_sparsity=None, max_nfev=None, verbose=0, args=(), kwargs=None):
+        jac_sparsity=None, max_nfev=None, verbose=0, args=(), kwargs=None,
+        callback=None
+):
     """Solve a nonlinear least-squares problem with bounds on the variables.
 
     Given the residuals f(x) (an m-D real function of n real
@@ -444,6 +449,26 @@ def least_squares(
         Additional arguments passed to `fun` and `jac`. Both empty by default.
         The calling signature is ``fun(x, *args, **kwargs)`` and the same for
         `jac`.
+    callback : None or callable, optional
+        Callback function that is called by the algorithm on each iteration.
+        This can be used to print or plot the optimization results at each
+        step, and to stop the optimization algorithm based on some user-defined
+        condition.  Only implemented for the `trf` and `dogbox` methods.
+
+        The signature is ``callback(intermediate_result: OptimizeResult)``
+
+        `intermediate_result is a `scipy.optimize.OptimizeResult`
+        which contains the intermediate results of the optimization at the
+        current iteration.
+
+        The callback also supports a signature like: ``callback(x)``
+
+        Introspection is used to determine which of the signatures is invoked.
+
+        If the `callback` function raises `StopIteration` the optimization algorithm
+        will stop and return with status code -2.
+
+        .. versionadded:: 1.16.0
 
     Returns
     -------
@@ -487,6 +512,7 @@ def least_squares(
         status : int
             The reason for algorithm termination:
 
+            * -2 : terminated because callback raised StopIteration.
             * -1 : improper input parameters status returned from MINPACK.
             *  0 : the maximum number of function evaluations is exceeded.
             *  1 : `gtol` termination condition is satisfied.
@@ -939,14 +965,20 @@ def jac_wrapped(x, f):
             else:
                 tr_solver = 'lsmr'
 
+    # Wrap callback function.  If callback is None, callback_wrapped also is None
+    callback_wrapped = _wrap_callback(callback)
+
     if method == 'lm':
+        if callback is not None:
+            warn("Callback function specified, but not supported with `lm` method.",
+                 stacklevel=2)
         result = call_minpack(fun_wrapped, x0, jac_wrapped, ftol, xtol, gtol,
                               max_nfev, x_scale, diff_step)
 
     elif method == 'trf':
         result = trf(fun_wrapped, jac_wrapped, x0, f0, J0, lb, ub, ftol, xtol,
                      gtol, max_nfev, x_scale, loss_function, tr_solver,
-                     tr_options.copy(), verbose)
+                     tr_options.copy(), verbose, callback=callback_wrapped)
 
     elif method == 'dogbox':
         if tr_solver == 'lsmr' and 'regularize' in tr_options:
@@ -958,7 +990,7 @@ def jac_wrapped(x, f):
 
         result = dogbox(fun_wrapped, jac_wrapped, x0, f0, J0, lb, ub, ftol,
                         xtol, gtol, max_nfev, x_scale, loss_function,
-                        tr_solver, tr_options, verbose)
+                        tr_solver, tr_options, verbose, callback=callback_wrapped)
 
     result.message = TERMINATION_MESSAGES[result.status]
     result.success = result.status > 0
diff --git a/scipy/optimize/_lsq/trf.py b/scipy/optimize/_lsq/trf.py
index 9154bdba5b2c..c72fbfae00f0 100644
--- a/scipy/optimize/_lsq/trf.py
+++ b/scipy/optimize/_lsq/trf.py
@@ -107,10 +107,11 @@
     CL_scaling_vector, compute_grad, compute_jac_scale, check_termination,
     update_tr_radius, scale_for_robust_loss_function, print_header_nonlinear,
     print_iteration_nonlinear)
+from scipy._lib._util import _call_callback_maybe_halt
 
 
 def trf(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale,
-        loss_function, tr_solver, tr_options, verbose):
+        loss_function, tr_solver, tr_options, verbose, callback=None):
     # For efficiency, it makes sense to run the simplified version of the
     # algorithm when no bounds are imposed. We decided to write the two
     # separate functions. It violates the DRY principle, but the individual
@@ -118,11 +119,11 @@ def trf(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale,
     if np.all(lb == -np.inf) and np.all(ub == np.inf):
         return trf_no_bounds(
             fun, jac, x0, f0, J0, ftol, xtol, gtol, max_nfev, x_scale,
-            loss_function, tr_solver, tr_options, verbose)
+            loss_function, tr_solver, tr_options, verbose, callback=callback)
     else:
         return trf_bounds(
             fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale,
-            loss_function, tr_solver, tr_options, verbose)
+            loss_function, tr_solver, tr_options, verbose, callback=callback)
 
 
 def select_step(x, J_h, diag_h, g_h, p, p_h, d, Delta, lb, ub, theta):
@@ -203,7 +204,8 @@ def select_step(x, J_h, diag_h, g_h, p, p_h, d, Delta, lb, ub, theta):
 
 
 def trf_bounds(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev,
-               x_scale, loss_function, tr_solver, tr_options, verbose):
+               x_scale, loss_function, tr_solver, tr_options, verbose, 
+               callback=None):
     x = x0.copy()
 
     f = f0
@@ -385,8 +387,20 @@ def trf_bounds(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev,
         else:
             step_norm = 0
             actual_reduction = 0
-
+            
         iteration += 1
+            
+        # Call callback function and possibly stop optimization
+        if callback is not None:
+            intermediate_result = OptimizeResult(
+                x=x, fun=f_true, nit=iteration, nfev=nfev)
+            intermediate_result["cost"] = cost
+
+            if _call_callback_maybe_halt(
+                callback, intermediate_result
+            ):
+                termination_status = -2
+                break
 
     if termination_status is None:
         termination_status = 0
@@ -399,7 +413,8 @@ def trf_bounds(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev,
 
 
 def trf_no_bounds(fun, jac, x0, f0, J0, ftol, xtol, gtol, max_nfev,
-                  x_scale, loss_function, tr_solver, tr_options, verbose):
+                  x_scale, loss_function, tr_solver, tr_options, verbose, 
+                  callback=None):
     x = x0.copy()
 
     f = f0
@@ -549,6 +564,18 @@ def trf_no_bounds(fun, jac, x0, f0, J0, ftol, xtol, gtol, max_nfev,
             actual_reduction = 0
 
         iteration += 1
+        
+        # Call callback function and possibly stop optimization
+        if callback is not None:
+            intermediate_result = OptimizeResult(
+                x=x, fun=f_true, nit=iteration, nfev=nfev)
+            intermediate_result["cost"] = cost
+
+            if _call_callback_maybe_halt(
+                callback, intermediate_result
+            ):
+                termination_status = -2
+                break
 
     if termination_status is None:
         termination_status = 0
diff --git a/scipy/optimize/tests/test_least_squares.py b/scipy/optimize/tests/test_least_squares.py
index 30401cd9dce2..14731ca13812 100644
--- a/scipy/optimize/tests/test_least_squares.py
+++ b/scipy/optimize/tests/test_least_squares.py
@@ -13,6 +13,7 @@
 from scipy.optimize._lsq.least_squares import IMPLEMENTED_LOSSES
 from scipy.optimize._lsq.common import EPS, make_strictly_feasible, CL_scaling_vector
 
+from scipy.optimize import OptimizeResult
 
 def fun_trivial(x, a=0):
     return (x - a)**2 + 5.0
@@ -787,6 +788,73 @@ def test_basic():
     assert_allclose(res.x, 0, atol=1e-10)
 
 
+def test_callback():
+    # test that callback function works as expected
+
+    results = []
+    
+    def my_callback_optimresult(intermediate_result: OptimizeResult):
+        results.append(intermediate_result)
+        
+    def my_callback_x(x):
+        r = OptimizeResult()
+        r.nit = 1
+        r.x = x
+        results.append(r)
+        return False
+        
+    def my_callback_optimresult_stop_exception(
+        intermediate_result: OptimizeResult):
+        results.append(intermediate_result)
+        raise StopIteration
+
+    def my_callback_x_stop_exception(x):
+        r = OptimizeResult()
+        r.nit = 1
+        r.x = x
+        results.append(r)
+        raise StopIteration
+
+    # Try for different function signatures and stop methods
+    callbacks_nostop = [my_callback_optimresult, my_callback_x]
+    callbacks_stop = [my_callback_optimresult_stop_exception,
+                      my_callback_x_stop_exception]
+    
+    # Try for all the implemented methods: trf, trf_bounds and dogbox
+    calls = [
+        lambda callback: least_squares(fun_trivial, 5.0, method='trf', 
+                                       callback=callback),
+        lambda callback: least_squares(fun_trivial, 5.0, method='trf', 
+                                       bounds=(-8.0, 8.0), callback=callback),
+        lambda callback: least_squares(fun_trivial, 5.0, method='dogbox', 
+                                       callback=callback)
+    ]
+
+    for mycallback, call in product(callbacks_nostop, calls):
+        results.clear()
+        # Call the different implemented methods
+        res = call(mycallback)
+        # Check that callback was called
+        assert len(results) > 0
+        # Check that results data makes sense
+        assert results[-1].nit > 0
+        # Check that it didn't stop because of the callback
+        assert res.status != -2
+        # final callback x should be same as final result
+        assert_allclose(results[-1].x, res.x)
+
+    for mycallback, call in product(callbacks_stop, calls):
+        results.clear()
+        # Call the different implemented methods
+        res = call(mycallback)
+        # Check that callback was called
+        assert len(results) > 0
+        # Check that only one iteration was run
+        assert results[-1].nit == 1
+        # Check that it stopped because of the callback
+        assert res.status == -2
+
+
 def test_small_tolerances_for_lm():
     for ftol, xtol, gtol in [(None, 1e-13, 1e-13),
                              (1e-13, None, 1e-13),