# Including Gradient Term in Optimization.jl

**URL:** <https://discourse.julialang.org/t/including-gradient-term-in-optimization-jl/126965>\
**Category:** General Usage\
**Tags:** optimization\
**Created:** [March 14, 2025, 7:56pm UTC](https://discourse.julialang.org/t/including-gradient-term-in-optimization-jl/126965 "2025-03-14T19:56:37Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![fdekerme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fdekerme/32/43574_2.png) [@fdekerme](https://discourse.julialang.org/u/fdekerme)\
**Post date:** [March 14, 2025, 7:56pm UTC](https://discourse.julialang.org/t/including-gradient-term-in-optimization-jl/126965/1 "2025-03-14T19:56:37Z")

</div>

Hi everyone 😀,

I’m trying to solve an optimization problem using `Optimization.jl`, but I’m struggling with incorporating a gradient term \nabla in the process. My cost function to minimize is:

\text{Cost} \left(T^{1}, T^{2}\right) = \|r^{A}\|^2 + \|r^{B}\|^2 + \|\nabla(T^{1})\|^2 + \|\nabla(T^{2})\|^2

where r^{A} and r^{B} involve log-sum-exponential terms:

r^{i}(T^{1}, T^{2}) = \log \left(\sum\_K K \exp \left(-T^{1} v\_K^{1} - T^{2} v\_K^{2} \right) \right) - \log \left(Y^{i}\right)

Y^{i}, T^{1} and T^{2} are matrix of size n x m

I first tried without the gradient term \nabla to solve the problem “pixel by pixel” (basically with 2 nested for loop), and it works. I then tried to vectorize my residual function r over the whole matrix, but then the optimization is extremely slow and doesn’t work anymore.

```julia
function f(K, v₁, v₂, T₁, T₂) # the theoretical value, compared to $Y^{i}$ the measurments
    x = sum(K[i] * exp.(-v₁[k] * T₁ - v₂[k] * T₂) for k in axes(K)[1])
    return -log.(x / length(axes(K)[1]))
end

function r(T, p)
    K, v₁, v₂, matrix_A, matrix_B, s = p
    T₁, T₂ = view(T, :, 1:s), view(T, :, (s+1):size(T, 2)) # I concatenated T₁, T₂ to be able to use Optimization.jl
    sum((f(K.A, v₁, v₂, T₁, T₂) .- matrix_A).^2 .+
    (f(K.B, v₁, v₂, T₁, T₂) .- matrix_B).^2)
end 

```

I think (?) vectorizing like that this is the only way to integrate the gradient term \nabla in the optimization, since the “pixel by pixel” approach does not give access to neighboring values in the matrix.  
I don’t know if all this is clear ( I tried haha 🙂 ), but does anyone have experience with this kind of problem or any relevant example?

Thanks in advance!  
fdekerm

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 14, 2025, 9:28pm UTC](https://discourse.julialang.org/t/including-gradient-term-in-optimization-jl/126965/2 "2025-03-14T21:28:44Z")

</div>

Hi! Can you provide a complete example of your code (with all imports, function and variable definitions), and explain what you mean by “gradient” here?

---

<div class="post-metadata">

**Author:** ![fdekerme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fdekerme/32/43574_2.png) [@fdekerme](https://discourse.julialang.org/u/fdekerme)\
**Post date:** [March 18, 2025, 9:34pm UTC](https://discourse.julialang.org/t/including-gradient-term-in-optimization-jl/126965/3 "2025-03-18T21:34:03Z")

</div>

I can’t share the inputs for the functions, unfortunately.  
For example in this [article](https://doi.org/10.1118/1.4866386), they propose the following optimization problem:

 ![Screenshot 2025-03-18 at 5.21.37 PM](https://global.discourse-cdn.com/julialang/original/3X/b/4/b4fc2e019f7f9bf3a649a23bc8a8a6eb9bb05942.png)  
with  
 ![Screenshot 2025-03-18 at 5.21.02 PM](https://global.discourse-cdn.com/julialang/original/3X/7/f/7fa6979efbba6cb492836045eb13c0d07a39206d.png)

and the “gradient”

 ![Screenshot 2025-03-18 at 5.25.06 PM](https://global.discourse-cdn.com/julialang/original/3X/8/e/8e477de5ef31b42fbcf78271243ff9d7034e7f6c.png)

They use images of N = 1024 \* 768 = 786432 pixels, so the matrix A (size 2N \*2N) is huge.  
Likewise, the operation ![Screenshot 2025-03-18 at 5.29.52 PM](https://global.discourse-cdn.com/julialang/original/3X/8/b/8b5db128b08a9504ddaa7052f93a9c94276034b8.png) is extremely heavy to realize on the whole matrix.

I’m not sure how to implement/vectorize the code to solve this type of problem with Optimiation.jl.  
Thanks in advance for your help,  
fdekerm

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 18, 2025, 11:28pm UTC](https://discourse.julialang.org/t/including-gradient-term-in-optimization-jl/126965/4 "2025-03-18T23:28:23Z")

</div>

> [@fdekerme](#):
>
> I’m not sure how to implement/vectorize the code to solve this type of problem with [Optimiation.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/Optimiation).  
> Thanks in advance for your help,

It’s a convex quadratic problem. Just take the derivative, set it equal to zero, and you have a system of linear equations that you can solve in a variety of ways, even for huge problems. In particular, it looks like the matrices here are all very sparse, so you can use sparse (or iterative) methods.

---

<div class="post-metadata">

**Author:** ![fdekerme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fdekerme/32/43574_2.png) [@fdekerme](https://discourse.julialang.org/u/fdekerme)\
**Post date:** [April 3, 2025, 9:39pm UTC](https://discourse.julialang.org/t/including-gradient-term-in-optimization-jl/126965/5 "2025-04-03T21:39:40Z")

</div>

Thank you very much for your answer! That is indeed what they are proposing in the article. I tried to implement it with [SparseMatrix](https://docs.julialang.org/en/v1/stdlib/SparseArrays/).  
The code (with an test example) is as follows:

```julia
using LinearAlgebra
using SparseArrays

function ∇R(nrows, ncols)
    N = nrows *ncols
    spdiagm(0 => 1.0*ones(N),
    -1 => -0.25*ones(N-1),
    1 => -0.25*ones(N-1),
    nrows => -0.25*ones(N-nrows),
    -nrows => -0.25*ones(N-nrows)
    )

end

function linear_solve_reg(images, μ₁, μ₂, λ)
    L = vcat(vec(images.top), vec(images.bottom))

    nrows, ncols = size(images.top)
    N = nrows * ncols
    R1 = ∇R(nrows, ncols)
    R = blockdiag(R1, R1)

    M = [μ₁.low μ₂.low; μ₁.high μ₂.high]
    A = kron(sparse(M), I(N)) + λ * R

    δ = A\L

    mat1 = reshape(δ[1:N], (nrows, ncols))
    mat2 = reshape(δ[N+1:end], (nrows, ncols))

    return mat1, mat2
end

# Test the function with fake data
nrows, ncols = 100, 100
images = (top = rand(nrows, ncols), bottom = rand(nrows, ncols))

μ₁ = (low = 0.3, high = 0.21)
μ₂ = (low = 0.19, high = 0.17)

λ = 0.1

mat1, mat2 = linear_solve_reg(images, μ₁, μ₂, λ)

```

I tried to do it with LinearSolve.jl, but the built-in Julia operator `\` seems to work even for sparse matrix. For some images, however, I get `NaN` in output, I don’t know why.

In addition to that, I tried to implement the convex quadratic problem in JuMP.jl to train myself. I manage to do it without the regularization, but it does not work anymore when I add it (divergence).

```julia
using JuMP, Ipopt, Images
function R(δ, i, j)
    R1 = δ[i, j, 1] - (δ[i-1, j, 1] + δ[i+1, j, 1] + δ[i, j-1, 1] + δ[i, j+1, 1])
    R2 = δ[i, j, 2] - (δ[i-1, j, 2] + δ[i+1, j, 2] + δ[i, j-1, 2] + δ[i, j+1, 2])
    return R1 + R2
end

model = Model(Ipopt.Optimizer)
A = [0.3 0.21; 0.19 0.17]
λ = 0.1
pad_top_image = collect(padarray(rand(100,150), Fill(0, (1, 1))))
pad_bottom_image = collect(padarray(rand(100,150), Fill(0, (1, 1))))
L = cat(pad_top_image, pad_bottom_image; dims=3)

x,y = size(pad_top_image)
@variable(model, δ[1:x, 1:y, 1:2])
@variable(model, z[1:x, 1:y, 1:2])

@constraint(model, [i=2:x-1, j=2:y-1], z[i,j,:] == A * δ[i, j, :] - L[i, j, :])
@objective(model, Min, sum(z.^2) + λ * sum(R(δ, i, j) for i=2:x-1, j=2:y-1))

optimize!(model)

δ_opt = value.(δ)
mat1 = δ_opt[:, :, 1]
mat2 = δ_opt[:, :, 2]

```

Thanks in advance for your help!  
fdekerm