finite_difference_jacobian! and forwarddiff_color_jacobian!

[NicolasLegouy]

Hi, I was wondering what is the difference between finite_difference_jacobian! and forwarddiff_color_jacobian! can I make use of both functions to solve my system f ?

using LinearAlgebra, ForwardDiff, SparseArrays, SparsityDetection, SparseDiffTools, BenchmarkTools, FiniteDiff

fcalls = 0 function f(dx,x) # in-place global fcalls += 1 for i in 2:length(x)-1 for j = i-1:i+1 dx[i] += 2x[i]*x[j] end end dx[1] = -2x[1] + x[2] dx[end] = x[end-1] - 2x[end]

nothing

dx end

input = rand(10) output = similar(input) sparsity_pattern = jacobian_sparsity(f,output,input) jac = Float64.(sparse(sparsity_pattern));

colors = matrix_colors(jac);

@time FiniteDiff.finite_difference_jacobian!(jac, f, rand(10), colorvec=colors) #Does not work

jac 10×10 SparseMatrixCSC{Float64,Int64} with 28 stored entries: [1 , 1] = -2.0 [2 , 1] = 0.813004 [1 , 2] = 1.0 [2 , 2] = 2.72799e7 [3 , 2] = 1.08469 [2 , 3] = 0.813004 [3 , 3] = 3.42811e7 [4 , 3] = 1.60417 [3 , 4] = 1.08469 [4 , 4] = 5.36339e7 [5 , 4] = 0.138877 [4 , 5] = 1.60417 ⋮ [6 , 6] = 8.37575e6 [7 , 6] = 1.14902 [6 , 7] = 0.265017 [7 , 7] = 3.84164e7 [8 , 7] = 0.0726532 [7 , 8] = 1.14902 [8 , 8] = 2.43784e6 [9 , 8] = 1.95844 [8 , 9] = 0.0726532 [9 , 9] = 6.18958e7 [10, 9] = 1.0 [9 , 10] = 1.95844 [10, 10] = -2.0

@time forwarddiff_color_jacobian!(jac, f, input, colorvec = colors) #Works

jac 10×10 SparseMatrixCSC{Float64,Int64} with 28 stored entries: [1 , 1] = -2.0 [2 , 1] = 1.65884 [1 , 2] = 1.0 [2 , 2] = 5.08291 [3 , 2] = 0.171734 [2 , 3] = 1.65884 [3 , 3] = 2.02608 [4 , 3] = 0.0237699 [3 , 4] = 0.171734 [4 , 4] = 0.895975 [5 , 4] = 0.676701 [4 , 5] = 0.0237699 ⋮ [6 , 6] = 5.01525 [7 , 6] = 1.32135 [6 , 7] = 1.5086 [7 , 7] = 4.84081 [8 , 7] = 0.689509 [7 , 8] = 1.32135 [8 , 8] = 4.03742 [9 , 8] = 1.33705 [8 , 9] = 0.689509 [9 , 9] = 3.42855 [10, 9] = 1.0 [9 , 10] = 1.33705 [10, 10] = -2.0

[ChrisRackauckas]

Try ]add FiniteDiff@2.4.1. The latest release may have had a bug?

[NicolasLegouy]

I get the same bug.

[ChrisRackauckas]

@kanav99 could you take a look at this one?

[kanav99]

Yes I will check this out at once

[kanav99]

@NicolasLegouy I changed your function and it works now I think, can you try too?

function f(dx,x) # in-place
    global fcalls += 1
    for i in 2:length(x)-1
        dx[i] = 0 # changed line
        for j = i-1:i+1
            dx[i] += 2x[i]*x[j]
        end
    end
    dx[1] = -2x[1] + x[2]
    dx[end] = x[end-1] - 2x[end]
    #nothing
    dx
end

also

@time FiniteDiff.finite_difference_jacobian!(jac, f, input, colorvec=colors) #Does not work

[ChrisRackauckas]

Was it just that explicit zeroing was needed?

[kanav99]

I think so.

With explicit zeroing, difference between the two jacobians was

Explored path: SparsityDetection.Path(Bool[], 1)
  0.906845 seconds (3.92 M allocations: 201.378 MiB, 4.81% gc time)
jac = 
  [1 ,  1]  =  -2.0
  [2 ,  1]  =  0.397597
  [1 ,  2]  =  1.0
  [2 ,  2]  =  2.49012
  [3 ,  2]  =  0.946422
  [2 ,  3]  =  0.397597
  [3 ,  3]  =  2.52617
  [4 ,  3]  =  0.235731
  [3 ,  4]  =  0.946422
  [4 ,  4]  =  1.72453
  [5 ,  4]  =  0.306642
  [4 ,  5]  =  0.235731
  [5 ,  5]  =  1.23027
  [6 ,  5]  =  0.381259
  [5 ,  6]  =  0.306642
  [6 ,  6]  =  2.12979
  [7 ,  6]  =  1.06063
  [6 ,  7]  =  0.381259
  [7 ,  7]  =  4.43299
  [8 ,  7]  =  1.93046
  [7 ,  8]  =  1.06063
  [8 ,  8]  =  5.88211
  [9 ,  8]  =  0.960544
  [8 ,  9]  =  1.93046
  [9 ,  9]  =  5.45547
  [10,  9]  =  1.0
  [9 , 10]  =  0.960544
  [10, 10]  =  -2.0
  1.421464 seconds (6.55 M allocations: 308.905 MiB, 6.20% gc time)
jac = 
  [1 ,  1]  =  -2.0
  [2 ,  1]  =  0.397597
  [1 ,  2]  =  1.0
  [2 ,  2]  =  2.49012
  [3 ,  2]  =  0.946422
  [2 ,  3]  =  0.397597
  [3 ,  3]  =  2.52617
  [4 ,  3]  =  0.235731
  [3 ,  4]  =  0.946422
  [4 ,  4]  =  1.72453
  [5 ,  4]  =  0.306642
  [4 ,  5]  =  0.235731
  [5 ,  5]  =  1.23027
  [6 ,  5]  =  0.381259
  [5 ,  6]  =  0.306642
  [6 ,  6]  =  2.12979
  [7 ,  6]  =  1.06063
  [6 ,  7]  =  0.381259
  [7 ,  7]  =  4.43299
  [8 ,  7]  =  1.93046
  [7 ,  8]  =  1.06063
  [8 ,  8]  =  5.88211
  [9 ,  8]  =  0.960544
  [8 ,  9]  =  1.93046
  [9 ,  9]  =  5.45547
  [10,  9]  =  1.0
  [9 , 10]  =  0.960544
  [10, 10]  =  -2.0
jac .- jac1 = 
  [1 ,  1]  =  -2.92319e-9
  [2 ,  1]  =  1.98521e-9
  [2 ,  2]  =  -3.17342e-8
  [3 ,  2]  =  -2.00617e-9
  [2 ,  3]  =  1.16866e-9
  [3 ,  3]  =  -2.85374e-8
  [4 ,  3]  =  3.83272e-10
  [3 ,  4]  =  4.34155e-9
  [4 ,  4]  =  -3.20952e-8
  [5 ,  4]  =  -3.09079e-10
  [4 ,  5]  =  3.83272e-10
  [5 ,  5]  =  -2.85287e-8
  [6 ,  5]  =  4.84518e-10
  [5 ,  6]  =  1.13424e-9
  [6 ,  6]  =  -3.14964e-8
  [7 ,  6]  =  1.79061e-9
  [6 ,  7]  =  -1.57888e-9
  [7 ,  7]  =  -2.66053e-8
  [8 ,  7]  =  2.62397e-8
  [7 ,  8]  =  1.79061e-9
  [8 ,  8]  =  -9.20157e-9
  [9 ,  8]  =  -1.71485e-8
  [8 ,  9]  =  1.79793e-8
  [9 ,  9]  =  -4.14528e-8
  [9 , 10]  =  1.70303e-8
  [10, 10]  =  -2.92319e-9

without -

Explored path: SparsityDetection.Path(Bool[], 1)
  0.887397 seconds (3.92 M allocations: 201.211 MiB, 4.47% gc time)
jac = 
  [1 ,  1]  =  -2.0
  [2 ,  1]  =  1.04234
  [1 ,  2]  =  1.0
  [2 ,  2]  =  3.49751e7
  [3 ,  2]  =  0.932773
  [2 ,  3]  =  1.04234
  [3 ,  3]  =  3.07857e7
  [4 ,  3]  =  1.8549
  [3 ,  4]  =  0.932773
  [4 ,  4]  =  4.61985e7
  [5 ,  4]  =  0.0755714
  [4 ,  5]  =  1.8549
  [5 ,  5]  =  2.53576e6
  [6 ,  5]  =  0.869497
  [5 ,  6]  =  0.0755714
  [6 ,  6]  =  2.86973e7
  [7 ,  6]  =  1.97536
  [6 ,  7]  =  0.869497
  [7 ,  7]  =  4.91986e7
  [8 ,  7]  =  1.11953
  [7 ,  8]  =  1.97536
  [8 ,  8]  =  3.75652e7
  [9 ,  8]  =  1.62704
  [8 ,  9]  =  1.11953
  [9 ,  9]  =  5.36995e7
  [10,  9]  =  1.0
  [9 , 10]  =  1.62704
  [10, 10]  =  -2.0
  1.467007 seconds (6.43 M allocations: 302.901 MiB, 13.38% gc time)
jac = 
  [1 ,  1]  =  -2.0
  [2 ,  1]  =  1.04234
  [1 ,  2]  =  1.0
  [2 ,  2]  =  4.38564
  [3 ,  2]  =  0.932773
  [2 ,  3]  =  1.04234
  [3 ,  3]  =  4.76279
  [4 ,  3]  =  1.8549
  [3 ,  4]  =  0.932773
  [4 ,  4]  =  4.71815
  [5 ,  4]  =  0.0755714
  [4 ,  5]  =  1.8549
  [5 ,  5]  =  2.87554
  [6 ,  5]  =  0.869497
  [5 ,  6]  =  0.0755714
  [6 ,  6]  =  3.78993
  [7 ,  6]  =  1.97536
  [6 ,  7]  =  0.869497
  [7 ,  7]  =  5.93975
  [8 ,  7]  =  1.11953
  [7 ,  8]  =  1.97536
  [8 ,  8]  =  5.84146
  [9 ,  8]  =  1.62704
  [8 ,  9]  =  1.11953
  [9 ,  9]  =  6.36431
  [10,  9]  =  1.0
  [9 , 10]  =  1.62704
  [10, 10]  =  -2.0
jac .- jac1 = 
  [1 ,  1]  =  3.23007e-9
  [2 ,  1]  =  1.10373e-8
  [2 ,  2]  =  -3.49751e7
  [3 ,  2]  =  -1.07887e-8
  [2 ,  3]  =  -4.50717e-9
  [3 ,  3]  =  -3.07857e7
  [4 ,  3]  =  8.59434e-9
  [3 ,  4]  =  -5.51638e-9
  [4 ,  4]  =  -4.61985e7
  [5 ,  4]  =  -9.88048e-11
  [4 ,  5]  =  -6.11676e-9
  [5 ,  5]  =  -2.53576e6
  [6 ,  5]  =  -1.13723e-8
  [5 ,  6]  =  3.33252e-10
  [6 ,  6]  =  -2.86973e7
  [7 ,  6]  =  -1.09474e-8
  [6 ,  7]  =  -1.065e-8
  [7 ,  7]  =  -4.91986e7
  [8 ,  7]  =  2.47997e-9
  [7 ,  8]  =  6.4964e-9
  [8 ,  8]  =  -3.75652e7
  [9 ,  8]  =  -5.8593e-9
  [8 ,  9]  =  -4.56148e-9
  [9 ,  9]  =  -5.36995e7
  [10,  9]  =  5.64277e-9
  [9 , 10]  =  -1.50162e-8
  [10, 10]  =  3.23007e-9

[ChrisRackauckas]

If we add an option to zero before each application of f, how much of a cost is that?

[kanav99]

That would be 4N additional broadcast calls, I think that's quiet a lot

[NicolasLegouy]

It works! Thank you for your prompt reply.

[ChrisRackauckas]

Let's close this, but remember that it could be an issue.