diff --git a/docs/src/api.md b/docs/src/api.md
index 17fe130e..f0c28e92 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -58,4 +58,4 @@ and its derivatives. Namely,
 |---------------------|---|
 | ``F(x)``            | [`residual`](@ref), [`residual!`](@ref) |
 | ``J_F(x)``          | [`jac_residual`](@ref), [`jac_coord_residual`](@ref), [`jac_coord_residual!`](@ref), [`jac_structure_residual`](@ref), [`jac_structure_residual!`](@ref), [`jprod_residual`](@ref), [`jprod_residual!`](@ref), [`jtprod_residual`](@ref), [`jtprod_residual!`](@ref), [`jac_op_residual`](@ref), [`jac_op_residual!`](@ref) |
-| ``\nabla^2 F_i(x)`` | [`hess_residual`](@ref), [`hess_coord_residual`](@ref), [`hess_coord_residual!`](@ref), [`hess_structure_residual`](@ref), [`hess_structure_residual!`](@ref), [`jth_hess_residual`](@ref), [`hprod_residual`](@ref), [`hprod_residual!`](@ref), [`hess_op_residual`](@ref), [`hess_op_residual!`](@ref) |
+| ``\nabla^2 F_i(x)`` | [`hess_residual`](@ref), [`hess_coord_residual`](@ref), [`hess_coord_residual!`](@ref), [`hess_structure_residual`](@ref), [`hess_structure_residual!`](@ref), [`jth_hess_residual`](@ref), [`jth_hess_residual_coord`](@ref), [`jth_hess_residual_coord!`](@ref), [`hprod_residual`](@ref), [`hprod_residual!`](@ref), [`hess_op_residual`](@ref), [`hess_op_residual!`](@ref) |
diff --git a/src/nls/api.jl b/src/nls/api.jl
index 0ac862b3..53292ebe 100644
--- a/src/nls/api.jl
+++ b/src/nls/api.jl
@@ -2,7 +2,7 @@ export residual, residual!, jac_residual, jac_structure_residual, jac_structure_
 export jac_coord_residual!, jac_coord_residual, jprod_residual, jprod_residual!
 export jtprod_residual, jtprod_residual!, jac_op_residual, jac_op_residual!
 export hess_residual, hess_structure_residual, hess_structure_residual!
-export hess_coord_residual!, hess_coord_residual, jth_hess_residual
+export hess_coord_residual!, hess_coord_residual, jth_hess_residual, jth_hess_residual_coord, jth_hess_residual_coord!
 export hprod_residual, hprod_residual!, hess_op_residual, hess_op_residual!
 
 """
@@ -325,10 +325,39 @@ Computes the Hessian of the j-th residual at x.
 """
 function jth_hess_residual(nls::AbstractNLSModel{T, S}, x::AbstractVector, j::Int) where {T, S}
   @lencheck nls.meta.nvar x
+  @rangecheck 1 nls.nls_meta.nequ j
+  rows, cols = hess_structure_residual(nls)
+  vals = jth_hess_residual_coord(nls, x, j)
+  return Symmetric(sparse(rows, cols, vals, nls.meta.nvar, nls.meta.nvar), :L)
+end
+
+"""
+    vals = jth_hess_residual_coord(nls, x, j)
+
+Evaluate the Hessian of j-th residual at `x` in sparse coordinate format.
+Only the lower triangle is returned.
+"""
+function jth_hess_residual_coord(nls::AbstractNLSModel{T, S}, x::AbstractVector, j::Int) where {T, S}
+  @lencheck nls.meta.nvar x
+  @rangecheck 1 nls.nls_meta.nequ j
+  vals = S(undef, nls.nls_meta.nnzh)
+  return jth_hess_residual_coord!(nls, x, j, vals)
+end
+
+"""
+    vals = jth_hess_residual_coord!(nls, x, j, vals)
+
+Evaluate the Hessian of j-th residual at `x` in sparse coordinate format, with `vals` of
+length `nls.nls_meta.nnzh`, in place. Only the lower triangle is returned.
+"""
+function jth_hess_residual_coord!(nls::AbstractNLSModel{T, S}, x::AbstractVector, j::Int, vals::AbstractVector) where {T, S}
+  @lencheck nls.meta.nvar x
+  @rangecheck 1 nls.nls_meta.nequ j
+  @lencheck nls.nls_meta.nnzh vals
   increment!(nls, :neval_jhess_residual)
   decrement!(nls, :neval_hess_residual)
   v = [i == j ? one(T) : zero(T) for i = 1:(nls.nls_meta.nequ)]
-  return hess_residual(nls, x, v)
+  return hess_coord_residual!(nls, x, v, vals)
 end
 
 """
@@ -342,8 +371,9 @@ function hprod_residual(
   i::Int,
   v::AbstractVector,
 ) where {T, S}
-  @lencheck nls.meta.nvar x
-  Hv = S(undef, nls_meta(nls).nvar)
+  @lencheck nls.meta.nvar x v
+  @rangecheck 1 nls.nls_meta.nequ i
+  Hv = S(undef, nls.meta.nvar)
   hprod_residual!(nls, x, i, v, Hv)
 end
 
@@ -361,6 +391,7 @@ Computes the Hessian of the i-th residual at x, in linear operator form.
 """
 function hess_op_residual(nls::AbstractNLSModel{T, S}, x::AbstractVector, i::Int) where {T, S}
   @lencheck nls.meta.nvar x
+  @rangecheck 1 nls.nls_meta.nequ i
   Hiv = S(undef, nls.meta.nvar)
   return hess_op_residual!(nls, x, i, Hiv)
 end
@@ -377,6 +408,7 @@ function hess_op_residual!(
   Hiv::AbstractVector,
 ) where {T, S}
   @lencheck nls.meta.nvar x Hiv
+  @rangecheck 1 nls.nls_meta.nequ i
   prod! = @closure (res, v, α, β) -> begin
     hprod_residual!(nls, x, i, v, Hiv)
     if β == 0