diff --git a/.github/workflows/demos.yml b/.github/workflows/demos.yml
index 3a19609e..67ca0853 100644
--- a/.github/workflows/demos.yml
+++ b/.github/workflows/demos.yml
@@ -69,11 +69,21 @@ jobs:
           pkg_path = dirname(Base.find_package("RegularizedOptimization"))
           include(joinpath(pkg_path, "..", "examples", "demo-bpdn.jl"))
         shell: julia --project="examples" --color=yes {0}
+      - name: Run contrained BPDN demo
+        run: |
+          pkg_path = dirname(Base.find_package("RegularizedOptimization"))
+          include(joinpath(pkg_path, "..", "examples", "demo-bpdn-constr.jl"))
+        shell: julia --project="examples" --color=yes {0}
       - name: Run FH demo
         run: |
           pkg_path = dirname(Base.find_package("RegularizedOptimization"))
           include(joinpath(pkg_path, "..", "examples", "demo-fh.jl"))
         shell: julia --project="examples" --color=yes {0}
+      - name: Run NNMF demo
+        run: |
+          pkg_path = dirname(Base.find_package("RegularizedOptimization"))
+          include(joinpath(pkg_path, "..", "examples", "demo-nnmf-constr.jl"))
+        shell: julia --project="examples" --color=yes {0}
       - name: Upload results
         uses: actions/upload-artifact@v3
         with:
diff --git a/examples/demo-bpdn-constr.jl b/examples/demo-bpdn-constr.jl
new file mode 100644
index 00000000..f5388bee
--- /dev/null
+++ b/examples/demo-bpdn-constr.jl
@@ -0,0 +1,30 @@
+using Random
+using LinearAlgebra
+using ProximalOperators
+using NLPModels, NLPModelsModifiers, RegularizedProblems, RegularizedOptimization
+using Printf
+
+include("plot-utils-bpdn.jl")
+
+Random.seed!(1234)
+
+function demo_solver(f, sol, h, χ, suffix = "l0-linf")
+  options = ROSolverOptions(ν = 1.0, β = 1e16, ϵa = 1e-6, ϵr = 1e-6, verbose = 10)
+
+  @info " using TR to solve with" h χ
+  reset!(f)
+  TR_out = TR(f, h, χ, options, x0 = f.meta.x0)
+  @info "TR relative error" norm(TR_out.solution - sol) / norm(sol)
+  plot_bpdn(TR_out, sol, "constr-tr-r2-$(suffix)")
+end
+
+function demo_bpdn_constr(compound = 1)
+  model, nls_model, sol = bpdn_model(compound, bounds = true)
+  f = LSR1Model(model)
+  λ = norm(grad(model, zeros(model.meta.nvar)), Inf) / 10
+  demo_solver(f, sol, NormL0(λ), NormLinf(1.0))
+  λ = norm(grad(model, zeros(model.meta.nvar)), Inf) / 3
+  demo_solver(f, sol, NormL1(λ), NormLinf(1.0), "l1-linf")
+end
+
+demo_bpdn_constr()
diff --git a/examples/demo-nnmf-constr.jl b/examples/demo-nnmf-constr.jl
new file mode 100644
index 00000000..0bbf3375
--- /dev/null
+++ b/examples/demo-nnmf-constr.jl
@@ -0,0 +1,29 @@
+using Random
+using LinearAlgebra
+using ProximalOperators
+using NLPModels, NLPModelsModifiers, RegularizedProblems, RegularizedOptimization
+using Printf
+
+include("plot-utils-nnmf.jl")
+
+Random.seed!(1234)
+
+function demo_solver(f, h, χ, selected, Avec, m, n, k, suffix = "l0-linf")
+  options = ROSolverOptions(ν = 1.0, β = 1e16, ϵa = 1e-6, ϵr = 1e-6, verbose = 10, maxIter = 500)
+  @info " using TR to solve with" h χ
+  reset!(f)
+  TR_out = TR(f, h, χ, options, selected = selected)
+  plot_nnmf(TR_out, Avec, m, n, k, "tr-r2-$suffix")
+end
+
+function demo_nnmf()
+  m, n, k = 100, 50, 5
+  model, A, selected = nnmf_model(m, n, k)
+  f = LSR1Model(model)
+  λ = norm(grad(model, rand(model.meta.nvar)), Inf) / 200
+  demo_solver(f, NormL0(λ), NormLinf(1.0), selected, A, m, n, k, "l0-linf")
+  λ = norm(grad(model, rand(model.meta.nvar)), Inf) / 10000
+  demo_solver(f, NormL1(λ), NormLinf(1.0), selected, A, m, n, k, "l1-linf")
+end
+
+demo_nnmf()
diff --git a/examples/plot-utils-bpdn.jl b/examples/plot-utils-bpdn.jl
index b49df969..6e57f6d0 100644
--- a/examples/plot-utils-bpdn.jl
+++ b/examples/plot-utils-bpdn.jl
@@ -30,4 +30,11 @@ function plot_bpdn(outstruct, sol, name = "tr-qr")
     ymode = "log",
   )
   save("bpdn-objdec-$(name).pdf", c)
+
+  d = Axis(
+    Plots.Linear(1:length(x), abs.(x - sol), mark = "none"),
+    xlabel = "index",
+    ylabel = "error \$|x - x^*|\$",
+  )
+  save("bpdn-error-$(name).pdf", d)
 end
diff --git a/examples/plot-utils-nnmf.jl b/examples/plot-utils-nnmf.jl
new file mode 100644
index 00000000..81137387
--- /dev/null
+++ b/examples/plot-utils-nnmf.jl
@@ -0,0 +1,34 @@
+using PGFPlots
+
+function plot_nnmf(outstruct, Avec, m, n, k, name = "tr-qr")
+  Comp_pg = outstruct.solver_specific[:SubsolverCounter]
+  objdec = outstruct.solver_specific[:Fhist] + outstruct.solver_specific[:Hhist]
+  x = outstruct.solution
+  A = reshape(Avec, m, n)
+  W = reshape(x[1:(m * k)], m, k)
+  H = reshape(x[(m * k + 1):end], k, n)
+  WH = W * H
+  
+  a = GroupPlot(2, 2, groupStyle = "horizontal sep = 2.5cm")
+  push!(a, Axis(Plots.Image(A, (1, m), (1, n), colormap=ColorMaps.Named("Jet")), xlabel = "A matrix (reference)"))
+  push!(a, Axis(Plots.Image(WH, (1, m), (1, n), colormap=ColorMaps.Named("Jet")), xlabel = "WH matrix"))
+  push!(a, Axis(Plots.Image(H, (1, k), (1, n), colormap=ColorMaps.Named("Jet")), xlabel = "H matrix"))
+  push!(a, Axis(Plots.Image(abs.(A - WH), (1, m), (1, n), colormap=ColorMaps.Named("Jet")), xlabel = "|A-WH| matrix"))
+  save("nnmf-$(name).pdf", a)
+
+  b = Axis(
+    Plots.Linear(1:length(Comp_pg), Comp_pg, mark = "none"),
+    xlabel = "outer iterations",
+    ylabel = "inner iterations",
+    ymode = "log",
+  )
+  save("nnmf-inner-outer-$(name).pdf", b)
+
+  c = Axis(
+    Plots.Linear(1:length(objdec), objdec, mark = "none"),
+    xlabel = "\$ k^{th}\$  \$ \\nabla f \$ Call",
+    ylabel = "Objective Value",
+    ymode = "log",
+  )
+  save("nnmf-objdec-$(name).pdf", c)
+end
diff --git a/src/TR_alg.jl b/src/TR_alg.jl
index ef925e11..3c0f68dc 100644
--- a/src/TR_alg.jl
+++ b/src/TR_alg.jl
@@ -34,7 +34,8 @@ The Hessian is accessed as an abstract operator and need not be the exact Hessia
 * `x0::AbstractVector`: an initial guess (default: `nlp.meta.x0`)
 * `subsolver_logger::AbstractLogger`: a logger to pass to the subproblem solver (default: the null logger)
 * `subsolver`: the procedure used to compute a step (`PG` or `R2`)
-* `subsolver_options::ROSolverOptions`: default options to pass to the subsolver (default: all defaut options).
+* `subsolver_options::ROSolverOptions`: default options to pass to the subsolver (default: all defaut options)
+* `selected::AbstractVector{<:Integer}`: (default `1:f.meta.nvar`).
 
 ### Return values
 
@@ -52,6 +53,7 @@ function TR(
   subsolver_logger::Logging.AbstractLogger = Logging.NullLogger(),
   subsolver = R2,
   subsolver_options = ROSolverOptions(),
+  selected::AbstractVector{<:Integer} = 1:f.meta.nvar,
 ) where {H, X}
   start_time = time()
   elapsed_time = 0.0
@@ -69,6 +71,12 @@ function TR(
   θ = options.θ
   β = options.β
 
+  local l_bound, u_bound
+  if has_bounds(f)
+    l_bound = f.meta.lvar
+    u_bound = f.meta.uvar
+  end
+
   if verbose == 0
     ptf = Inf
   elseif verbose == 1
@@ -81,11 +89,11 @@ function TR(
 
   # initialize parameters
   xk = copy(x0)
-  hk = h(xk)
+  hk = h(xk[selected])
   if hk == Inf
     verbose > 0 && @info "TR: finding initial guess where nonsmooth term is finite"
     prox!(xk, h, x0, one(eltype(x0)))
-    hk = h(xk)
+    hk = h(xk[selected])
     hk < Inf || error("prox computation must be erroneous")
     verbose > 0 && @debug "TR: found point where h has value" hk
   end
@@ -93,7 +101,8 @@ function TR(
 
   xkn = similar(xk)
   s = zero(xk)
-  ψ = shifted(h, xk, Δk, χ)
+  ψ = has_bounds(f) ? shifted(h, xk, max.(-Δk, l_bound - xk), min.(Δk, u_bound - xk), selected) :
+    shifted(h, xk, Δk, χ)
 
   Fobj_hist = zeros(maxIter)
   Hobj_hist = zeros(maxIter)
@@ -160,7 +169,8 @@ function TR(
 
     subsolver_options.ϵa = k == 1 ? 1.0e-5 : max(ϵ, min(1e-2, sqrt(ξ1)) * ξ1)
     ∆_effective = min(β * χ(s), Δk)
-    set_radius!(ψ, ∆_effective)
+    has_bounds(f) ? set_bounds!(ψ, max.(-∆_effective, l_bound - xk), min.(∆_effective, u_bound - xk)) :
+      set_radius!(ψ, ∆_effective)
     s, iter, _ = with_logger(subsolver_logger) do
       subsolver(φ, ∇φ!, ψ, subsolver_options, s)
     end
@@ -169,7 +179,7 @@ function TR(
     sNorm = χ(s)
     xkn .= xk .+ s
     fkn = obj(f, xkn)
-    hkn = h(xkn)
+    hkn = h(xkn[selected])
     hkn == -Inf && error("nonsmooth term is not proper")
 
     Δobj = fk + hk - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
@@ -191,11 +201,12 @@ function TR(
 
     if η2 ≤ ρk < Inf
       Δk = max(Δk, γ * sNorm)
-      set_radius!(ψ, Δk)
+      !has_bounds(f) && set_radius!(ψ, Δk)
     end
 
     if η1 ≤ ρk < Inf
       xk .= xkn
+      has_bounds(f) && set_bounds!(ψ, max.(-Δk, l_bound - xk), min.(Δk, u_bound - xk))
 
       #update functions
       fk = fkn
@@ -214,7 +225,7 @@ function TR(
 
     if ρk < η1 || ρk == Inf
       Δk = Δk / 2
-      set_radius!(ψ, Δk)
+      has_bounds(f) ? set_bounds!(ψ, max.(-Δk, l_bound - xk), min.(Δk, u_bound - xk)) : set_radius!(ψ, Δk)
     end
     tired = k ≥ maxIter || elapsed_time > maxTime
   end
diff --git a/test/runtests.jl b/test/runtests.jl
index 328e2e9f..1dc1db2f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,6 +7,7 @@ const global compound = 1
 const global nz = 10 * compound
 const global options = ROSolverOptions(ν = 1.0, β = 1e16, ϵa = 1e-6, ϵr = 1e-6, verbose = 10)
 const global bpdn, bpdn_nls, sol = bpdn_model(compound)
+const global bpdn2, bpdn_nls2, sol2 = bpdn_model(compound, bounds = true)
 const global λ = norm(grad(bpdn, zeros(bpdn.meta.nvar)), Inf) / 10
 
 for (mod, mod_name) ∈ ((x -> x, "exact"), (LSR1Model, "lsr1"), (LBFGSModel, "lbfgs"))
@@ -109,3 +110,5 @@ for (h, h_name) ∈ ((NormL1(λ), "l1"),)
     @test LMTR_out.status == :first_order
   end
 end
+
+include("test_bounds.jl")
diff --git a/test/test_bounds.jl b/test/test_bounds.jl
new file mode 100644
index 00000000..cbd97955
--- /dev/null
+++ b/test/test_bounds.jl
@@ -0,0 +1,27 @@
+for (mod, mod_name) ∈ ((x -> x, "exact"), (LSR1Model, "lsr1"), (LBFGSModel, "lbfgs"))
+  for (h, h_name) ∈ ((NormL0(λ), "l0"), (NormL1(λ), "l1"))
+    for solver_sym ∈ (:TR,)
+      solver_sym == :TR && mod_name == "exact" && continue
+      solver_name = string(solver_sym)
+      solver = eval(solver_sym)
+      @testset "bpdn-with-bounds-$(mod_name)-$(solver_name)-$(h_name)" begin
+        x0 = zeros(bpdn2.meta.nvar)
+        p = randperm(bpdn2.meta.nvar)[1:nz]
+        args = solver_sym == :R2 ? () : (NormLinf(1.0),)
+        @test has_bounds(mod(bpdn2))
+        out = solver(mod(bpdn2), h, args..., options, x0 = x0)
+        @test typeof(out.solution) == typeof(bpdn2.meta.x0)
+        @test length(out.solution) == bpdn2.meta.nvar
+        @test typeof(out.solver_specific[:Fhist]) == typeof(out.solution)
+        @test typeof(out.solver_specific[:Hhist]) == typeof(out.solution)
+        @test typeof(out.solver_specific[:SubsolverCounter]) == Array{Int, 1}
+        @test typeof(out.dual_feas) == eltype(out.solution)
+        @test length(out.solver_specific[:Fhist]) == length(out.solver_specific[:Hhist])
+        @test length(out.solver_specific[:Fhist]) == length(out.solver_specific[:SubsolverCounter])
+        @test obj(bpdn2, out.solution) == out.solver_specific[:Fhist][end]
+        @test h(out.solution) == out.solver_specific[:Hhist][end]
+        @test out.status == :first_order
+      end
+    end
+  end
+end