@@ -477,7 +477,8 @@ adjust it properly through the XML.
    gradient norm. For Newton, it is additionally applied to the Newton decrement :math:`\tfrac{1}{2} g^T H^{-1} g`, the
    predicted cost improvement of the next iteration. Before the first iteration, CG and Newton also apply it to a
    :ref:`convergence certificate<soAlgorithms>` of the warmstarted solution, possibly terminating with zero iterations.
-   Set the tolerance to 0 to disable early termination.
+   The cost improvement and the gradient are scaled by the inverse trace of the inertia matrix, of the whole model or
+   of a constraint :ref:`island<soIsland>`. Set the tolerance to 0 to disable early termination.
 
 .. _option-ls_iterations:
 
@@ -1253,9 +1254,11 @@ parameters.
 .. _statistic-meaninertia:
 
 :at:`meaninertia`: :at-val:`real, optional`
-   If this attribute is specified, it replaces the value of mjModel.stat.meaninertia computed by the compiler. The
-   computed value is the average diagonal element of the joint-space inertia matrix when the model is in qpos0. At
-   runtime this value scales the solver cost and gradient used for early termination.
+   This attribute is deprecated and will be removed in a future release; compiling a model which sets it gives a
+   warning. If it is specified, it replaces the value of mjModel.stat.meaninertia computed by the compiler. The
+   computed value is the average diagonal element of the joint-space inertia matrix when the model is in qpos0. It
+   does not affect simulation: the solvers scale their termination criteria by the trace of the current inertia
+   matrix, see :ref:`tolerance<option-tolerance>`.
 
 .. _statistic-meansize:
 
 
@@ -10,6 +10,12 @@ Engine
 - Trees marked :ref:`sleep="init"<body-sleep>` which cannot be put to sleep now cause a warning rather than an error,
   and start awake. See :ref:`Sleeping<siSleepSleeping>` for details.
 - Added single-shot :ref:`multicontact<coMultiCCD>` for geoms with positive contact :ref:`margins<body-geom-margin>`.
+- The constraint solvers scale their termination criteria by the trace of the current inertia matrix, of the model or
+  of a constraint island, instead of ``mjModel.stat.meaninertia`` times ``nv``, which was computed in ``qpos0`` and used
+  by CG and Newton only without islands. Models solved without islands, or with PGS or noslip, can stop at slightly
+  different iterations. The :ref:`statistic<statistic>` attribute :ref:`meaninertia<statistic-meaninertia>`, which no
+  longer affects simulation, is deprecated and will be removed in a future release (compiling a model that sets it now
+  gives a warning; remove the attribute to migrate).
 
 Compiler
 ^^^^^^^^
 
@@ -64,7 +64,7 @@ typedef struct mjTimerStat_ {  // timer statistics
   int     number;              // how many times was timer called
 } mjTimerStat;
 typedef struct mjSolverStat_ {  // per-iteration solver statistics
-  mjtNum  improvement;          // cost reduction, scaled by 1/trace(M(qpos0))
+  mjtNum  improvement;          // cost reduction, scaled by 1/trace(M)
   mjtNum  gradient;             // gradient norm (primal only, scaled)
   mjtNum  lineslope;            // slope in linesearch
   int     nactive;              // number of active constraints
 
@@ -84,7 +84,7 @@ typedef struct mjTimerStat_ {  // timer statistics
 
 
 typedef struct mjSolverStat_ {  // per-iteration solver statistics
-  mjtNum  improvement;          // cost reduction, scaled by 1/trace(M(qpos0))
+  mjtNum  improvement;          // cost reduction, scaled by 1/trace(M)
   mjtNum  gradient;             // gradient norm (primal only, scaled)
   mjtNum  lineslope;            // slope in linesearch
   int     nactive;              // number of active constraints
 
@@ -322,16 +322,8 @@ def _put_option(
   raise NotImplementedError(f'Unsupported implementation: {impl}')
 
 
-def _put_statistic(
-    s: mujoco.MjStatistic, impl: types.Impl
-) -> Union[types.Statistic, types.StatisticWarp]:
+def _put_statistic(s: mujoco.MjStatistic) -> types.Statistic:
   """Puts mujoco.MjStatistic onto a device, resulting in mjx.Statistic."""
-  if impl == types.Impl.WARP:
-    fields = {
-        f.name: _wp_to_np_type(getattr(s, f.name, None), f.name)
-        for f in types.StatisticWarp.fields()
-    }
-    return types.StatisticWarp(**fields)
   return types.Statistic(
       meaninertia=s.meaninertia,
       meanmass=s.meanmass,
@@ -424,7 +416,7 @@ def _put_model_jax(
   fields = {f: getattr(m, f) for f in mj_field_names}
   fields['cam_mat0'] = fields['cam_mat0'].reshape((-1, 3, 3))
   fields['opt'] = _put_option(m.opt, types.Impl.JAX)
-  fields['stat'] = _put_statistic(m.stat, types.Impl.JAX)
+  fields['stat'] = _put_statistic(m.stat)
 
   fields_jax = {}
   fields_jax['dof_hasfrictionloss'] = fields['dof_frictionloss'] > 0
@@ -495,7 +487,7 @@ def _put_model_warp(
   private_options = {k: getattr(mw.opt, k) for k in option_keys}
   private_options['graph_mode'] = graph_mode
   fields['opt'] = _put_option(m.opt, types.Impl.WARP, private_options)
-  fields['stat'] = _put_statistic(m.stat, types.Impl.WARP)
+  fields['stat'] = _put_statistic(m.stat)
 
   # Use MJW fields directly instead of MjModel, so that shape and dtype are
   # always compatible with MJXW (e.g. cam_mat0/geom_aabb).
@@ -544,7 +536,7 @@ def _put_model_cpp(
   fields = {f: getattr(m, f) for f in mj_field_names}
   fields['cam_mat0'] = fields['cam_mat0'].reshape((-1, 3, 3))
   fields['opt'] = _put_option(m.opt, impl=types.Impl.JAX)
-  fields['stat'] = _put_statistic(m.stat, impl=types.Impl.JAX)
+  fields['stat'] = _put_statistic(m.stat)
 
   # get the pointer address
   # we use a 0-d array
 
@@ -418,8 +418,21 @@ def _update_gradient(m: Model, d: Data, ctx: Context) -> Context:
   return ctx
 
 
-def _rescale(m: Model, value: jax.Array) -> jax.Array:
-  return value / (m.stat.meaninertia * max(1, m.nv))
+def _trace_m(m: Model, d: Data) -> jax.Array:
+  if not isinstance(d._impl, DataJAX):
+    raise ValueError('_trace_m requires JAX backend implementation.')
+  if not m.nv:
+    return jp.array(mujoco.mjMINVAL)
+  trace = (
+      d._impl.M[m.M_rowadr + m.M_rownnz - 1].sum()
+      if support.is_sparse(m)
+      else jp.trace(d._impl.M)
+  )
+  return jp.maximum(mujoco.mjMINVAL, trace)
+
+
+def _rescale(m: Model, d: Data, value: jax.Array) -> jax.Array:
+  return value / _trace_m(m, d)
 
 
 def _linesearch(m: Model, d: Data, ctx: Context) -> Context:
@@ -440,7 +453,7 @@ def _linesearch(m: Model, d: Data, ctx: Context) -> Context:
   ):
     raise ValueError('_lineasearch requires JAX backend implementation.')
 
-  smag = math.norm(ctx.search) * m.stat.meaninertia * max(1, m.nv)
+  smag = math.norm(ctx.search) * _trace_m(m, d)
   gtol = m.opt.tolerance * m.opt.ls_tolerance * smag
 
   # compute Mv, Jv
@@ -582,8 +595,8 @@ def solve(m: Model, d: Data) -> Data:
     raise ValueError('solve requires JAX backend implementation.')
 
   def cond(ctx: Context) -> jax.Array:
-    improvement = _rescale(m, ctx.prev_cost - ctx.cost)
-    gradient = _rescale(m, math.norm(ctx.grad))
+    improvement = _rescale(m, d, ctx.prev_cost - ctx.cost)
+    gradient = _rescale(m, d, math.norm(ctx.grad))
 
     done = ctx.solver_niter >= m.opt.iterations
     done |= improvement < m.opt.tolerance
 
@@ -475,13 +475,6 @@ class Statistic(PyTreeNode):
   center: jax.Array
 
 
-class StatisticWarp(mjxw_types.StatisticWarp, Statistic):
-  """Warp-specific model statistics."""
-
-  # NB: StatisticWarp type annotations may not match those on Statistic.
-  pass
-
-
 class OptionJAX(PyTreeNode):
   """JAX-specific option."""
 
@@ -654,7 +647,7 @@ class Model(PyTreeNode):
   npluginstate: int
   nhistory: int
   opt: Option
-  stat: Union[Statistic, StatisticWarp]
+  stat: Statistic
   qpos0: jax.Array
   qpos_spring: jax.Array
   body_parentid: np.ndarray
 
@@ -196,7 +196,6 @@ def _global_assignments():
       ('_m', 'mjwarp.Model'),
       ('_d', 'mjwarp.Data'),
       ('_o', 'mjwarp.Option'),
-      ('_s', 'mjwarp.Statistic'),
       ('_c', 'mjwarp.Contact'),
       ('_e', 'mjwarp.Constraint'),
       ('_cb', 'mjwp_types.Callback'),
@@ -325,7 +324,7 @@ def _jax_shim_fn(
     elif mjwarp_field_info[arg].param_source == 'Model':
       arg_jax = (
           arg.replace('__', '.')
-          if arg.startswith('opt') or arg.startswith('stat')
+          if arg.startswith('opt')
           else arg
       )
       public_field = arg.split('__')[0] in mjx_types.Model.__annotations__
@@ -451,7 +450,6 @@ def create_jax_warp_shim(
 def _{fn_name}_shim(
 {fn_args_raw_str}
 ):
-  _m.stat = _s
   _m.opt = _o
   _m.callback = _cb
   _d.efc = _e
 
@@ -49,7 +49,6 @@
     'Model': mjwarp.Model,
     'Data': mjwarp.Data,
     'Option': mjwarp.Option,
-    'Statistic': mjwarp.Statistic,
 }
 
 
@@ -413,7 +412,7 @@ def write_core_cls(
     extra_annotations: dict[str, type] | None = None,
     defaults: Optional[Dict[str, Any]] = None,
 ):
-  """Writes a core API class (e.g. Model/Data/Option/Statistic)."""
+  """Writes a core API class (e.g. Model/Data/Option)."""
   cls = _CLS_MAP[cls_name]
 
   annotations = dict(cls.__annotations__)
@@ -562,7 +561,6 @@ def main(argv):
   write_nested_dataclass(target_fpath, mjwarp._src.types.TileSet)
   write_nested_dataclass(target_fpath, mjwarp._src.types.BlockDim)
 
-  write_core_cls('Statistic', target_fpath, mjx_types_fpath, set_diff=False)
   write_core_cls(
       'Option', target_fpath, mjx_types_fpath,
       extra_annotations={'graph_mode': GraphMode},
 
@@ -287,7 +287,6 @@ def get_mjwarp_field_info(
   """Return field info for mujoco_warp/_src/types.py."""
   dataclass_map = {
       'opt': 'Option',
-      'stat': 'Statistic',
       'efc': 'Constraint',
       'contact': 'Contact',
   }
@@ -314,7 +313,7 @@ def get_fields(cls_name, base_order, sub_order):
   # to Model or Data based solely on the field name lookup in mjwarp_field_info
   # (via field_info.param_source) rather than inspecting base variable names.
   #
-  # Nested sub-dataclasses (Option, Statistic in Model; Constraint, Contact in
+  # Nested sub-dataclasses (Option in Model; Constraint, Contact in
   # Data) are flattened with container prefixes (e.g. 'opt__timestep',
   # 'contact__pos'). This prefixing prevents collisions between sub-dataclasses
   # that share generic field names (such as 'pos' or 'type' in Contact and