JuliaMath / JuliaMath/Calculus.jl

Calculus.hessian is not symmetric

Open
#91 6 comments 0 reactions 0 assignees View on GitHub
Dominant language
Julia
Stars
285
Forks
78
PR merge metrics
No merged PRs in 30d

Description

```
julia> @time h10 = Calculus.hessian(RB.mylike, RB.troubleParamsSim1.raw)
209.557745 seconds (2.43 G allocations: 83.867 GB, 10.31% gc time)
10x10 Array{Float64,2}:
566.715 -206.282 -50.7658 … 3135.03 -89.8427
-206.282 566.715 17.7032 -1993.75 47.1135
-50.7658 17.7032 143.437 -1163.21 33.7016
17.7032 -50.7658 -47.7584 282.533 1.3897
182.291 -122.875 -12.5876 -2943.26 27.7885
-2.95777 -5.09084 2.97638 … -452.285 15.1723
-45.7494 12.5752 39.2448 0.613878 24.8811
-1942.38 804.496 450.315 -35319.9 1935.28
3135.03 -1993.75 -1163.21 56057.6 -2498.87
-89.8427 47.1135 33.7016 -2498.87 1041.54

julia> issym(h10)
false
```

Closer inspection of h10 shows all the errors are extremely small; the largest is 4e-8.

My original assumption was that this was the result of numerical noise from calculating the derivatives in different orders. But the code (`finite_difference_hessian!` in `finite_difference.jl`) appears to compute only the upper triangle,

```
for i = 1:n
#code omitted
# and i+1:n appears to be parses (i+1):n, which is right
for j = i+1:n
```

and then ends with

```
Base.LinAlg.copytri!(H,'U')
```

So maybe the problem is in `copytri!`.

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.