2

Could you please provide a solution for finding a Jacobian matrix for a vector function using the "IntervalArithmetic.jl" package in Julia? The function representation and the required Jacobian matrix is shown below.

f = [f1(x,y,z), f2(x,y,z), f3(x,y,z)]
Jacobian = [df1/dx, df1/dy, df1/dz; df2/dx, df2/dy, df2/dz; df3/dx, df3/dy, df3/dz]
James Z
  • 12,209
  • 10
  • 24
  • 44
Kanishk Sharma
  • 173
  • 7

1 Answers1

0

Here you go:

(the explanation of this code is left as a self exercise to the user of this code, which should be able to do it, otherwise using this code is a little risky - it isn't so well tested)

using IntervalArithmetic, FiniteDiff, Optim

function interval_jacobian(F, x)
    x0 = mid.(x)
    F0 = F(x0)
    dummy = zero(eltype(x0))/one(eltype(F0)) 
    restype = typeof(dummy)
    M = length(F0)
    N = length(x0)
    lows = fill(dummy, M, N)
    highs = fill(dummy, M, N)
    xlo = getproperty.(x, :lo)
    xhi = getproperty.(x, :hi)
    inner_optimizer = GradientDescent()
    for i in 1:M
        for j in 1:N
            ff = (x -> FiniteDiff.finite_difference_derivative(t -> F(Base.setindex(tuple(x...), t, j))[i], x[j]))
            gg! = ((G, x) -> ( G .= FiniteDiff.finite_difference_hessian(x -> F(x)[i], x)[:, j]))
            res = optimize(ff, gg!, xlo, xhi, x0, Fminbox(inner_optimizer))
            lows[i,j] = res.minimum
            ff2 = (x -> -ff(x))
            gg2! = ((G, x) -> (gg!(G, x); G .= -G))
            res = optimize(ff2, gg2!, xlo, xhi, x0, Fminbox(inner_optimizer))
            highs[i,j] = -res.minimum
        end
    end
    return [l..h for (l,h) in zip(lows, highs)]
end

With this definition:

# defining a function taking a single tuple parameter

g((x,y)) = [ 3x^2 + 4y^3 + x*y^2, 2x + y^4, x^2 + y^2 ]
# g (generic function with 1 method)

# and another one:

g2((x,y)) = [ 3sin(x)^4*cos(y), 2cos(x)^5+y^3 ]
# g2 (generic function with 1 method)

# some intervals as inputs:

xintervals = [1.0..2.0, 2.0..3.0]
# 2-element Vector{Interval{Float64}}:
# [1, 2]
# [2, 3]

# the interval Jacobian calculated:

interval_jacobian(g, xintervals)
# 3×2 Matrix{Interval{Float64}}:
#  [10, 21]            [52, 120]
#     [1.99999, 2]  [32, 108]
#   [2, 4]               [4, 6]

interval_jacobian(g2, xintervals)
# 2×2 Matrix{Interval{Float64}}:
# [-3.85812, 3.71688]       [-2.7279, -0.212259]
# [-0.717112, 0]       [12, 27]

This is not meant to be performant, or symbolic. Issues with infinities not dealt with, or just issues such as functions being too unstable for Optim optimizer.

Dan Getz
  • 17,002
  • 2
  • 23
  • 41
  • Thank you @Dan. The provided code is working with the constant coefficients. But it fails to provide a solution if coefficients are in interval form. For example: g((x,y)) = [ (3..4)*x^2 + (4..5)*y^3 + x*y^2, (2..3)*x + y^4, x^2 + y^2 ]. How can this issue be addressed? – Kanishk Sharma Jan 24 '23 at 03:41
  • There are enough details in the solution to address this. Apart from running the function, do you follow its operation? – Dan Getz Jan 24 '23 at 04:40
  • 1
    Okay sir, I will try. I am new to learning Julia, specifically interval-based analysis. So, I am trying to understanding the code you gave. – Kanishk Sharma Jan 25 '23 at 18:15