5 ms·

Just now I was thinking of moving a long calculation from R to Julia (non-linear optimisation of a simple function with multiple local minima, for a lot of diff

Just now I was thinking of moving a long calculation from R to Julia (non-linear optimisation of a simple function with multiple local minima, for a lot of different datasets). No loops. Embarrassingly parallel. And to my great surprise, R and Julia took the same time.

ChrisRackauckas 4y ago
Did you use JuMP? It would be interesting to see the JuMP code.
- BrandonS113 4y ago
  No, nlopt. Why we could easy port from R to Julia as nlopt exists for both (its c)
  - ChrisRackauckas 4y ago
    Ahh then it's all based on the speed of the objective function and the gradient calculation code. Do you have an example to look at for that? Or did you try modelingtoolkitizing and running an auto-simplified form?
    - BrandonS113 4y ago
      Yes, well no gradients. It a simple parametric function (6 parameters) applied to a vector of lengths 10 to 20. All vectorizes in the language of R. All powers, logs, exponentials, ratios, and sums. Large number of local maxima, and places where function cannot be evaluated. So probabilistic optimizer is very important. Come to think of it, for this sort of calculations, R and Julia should take the same time.
      - DNF2 4y ago
        Of course, this is piques one's curiosity. It might be, if the function is simple enough, that there is little advantage to Julia here. But if you are combining multiple operations on a vector, there could be opportunities for Julia, in-place operations, fusing, simd. Maybe even StaticArrays. Any chance of sharing that little piece of code?
        BrandonS113 4y ago
        sure, for parameters par (what is optimized for), data vector x (typical length from 10 to 20), constants n and n2, a typical function is if((1 - par[3]^2)<0 ) return(100) if(par[1] + par[2] \* par[5] \* sqrt(1 - par[3]^2)<0) return(500 ) tmp <- (x - par[4])^2 + par[5]^2 tmp2 <- par[1] + par[2] \* (par[3] \* (x - par[4]) + sqrt(tmp)) result <- ((sqrt(tmp2) - Ce) / n2)^2 result <- sqrt(sum(n \* result) / n) if(is.na(result)) return(1111) return(result)
        DNF2 4y ago
        I'm not that good at reading R. But if the Julia code is similar, then this code is type unstable, sometimes returning an Int, sometimes a Float. That harms performance. Generally, it looks like a function where Julia could have a significant performance advantage.
        BrandonS113 4y ago
        R looks pretty much the same as Julia here. Don't think its unstable. Its all float. Nlopt passes parameters as float, data is float. Can't see how it could return int. But I'll take a second look at the julia version of it. We optimised R as much as possible. might have missed a julia trick. simd perhaps. Except we also run it on ARM
        DNF2 4y ago
        Maybe it's all floats in R, but in Julia, `return 500` means an Int is returned. But it's really hard to determine whether the Julia code is unstable based only on the R code.
    - celrod 4y ago
      I've noticed that R often defaults to much higher tolerances than Julia, even when it's wrappers to the same C library, like cubature. R cubature[0]: 1e-5 Cubature.jl [1]: 1e-8 The difference for NLopt in R vs Julia is smaller. `NLopt.DEFAULT_OPTIONS`[2] in Julia shows `1e-7` for `ftol_rel`, `xtol_rel`, and `constrtol_abs`, while in R `xtol_rel` is `1e-6` and the others are `0.0`[3]. So, the options at least aren't the same with nlopt. Anyway, I always recommend confirming that you're comparing the same settings. And of course, in Julia, you'll probably want to `JET.report_opt` your function and fix and glaring performance issues. NLopt seems like it may be a bit of an exception, but I noticed this is pretty common pattern elsewhere, uniroot[4] being another example, with eps()^(1/4) default tolerance, far higher than Julia root solvers will use. [0] https://cran.r-project.org/web/packages/cubature/cubature.pdf https://cran.r-project.org/web/packages/cubature/cubature.pd... [1] https://github.com/JuliaMath/Cubature.jl https://github.com/JuliaMath/Cubature.jl [2] https://github.com/JuliaOpt/NLopt.jl/blob/6ade25740362895bbfff1aee07d35911a6e2df17/src/MOI_wrapper.jl#L315 https://github.com/JuliaOpt/NLopt.jl/blob/6ade25740362895bbf... [3] https://cran.r-project.org/web/packages/nloptr/nloptr.pdf https://cran.r-project.org/web/packages/nloptr/nloptr.pdf [4] https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/uniroot https://www.rdocumentation.org/packages/stats/versions/3.6.2...
      - BrandonS113 4y ago
        We take of that by setting all the NLOPT options to be the same across calling languages. We get pretty much the same number of calls in julia, r and c.
  - harshreality 4y ago
    If you use the same native C library from both Julia and some other language, and (presumably) that library is where the algorithm spends the vast majority of its time, why would you expect Julia to be faster?
    - BrandonS113 4y ago
      no. the code spends 99% of its time in 5 lines in the objective function. And my usual experience is that Julia is much faster than R. Just not always, apparently.
Lyngbakr 4y ago
I only used Julia for a short time, but I didn't see the blazing fast speeds I was promised. I've seen the benchmarks, of course, on which the claims are founded, but the C-like speeds weren't obvious to me in everyday data science workflows. In the end, there wasn't sufficient motivation for me to switch to Julia as my weapon of choice. I do like Pluto[0], though... [0]https://plutojl.org/ https://plutojl.org/
- BrandonS113 4y ago
  If doing data science, I find Julia's tools to be inferior to Python and R. But in my work, when it comes to long computations, not only does Julia usually vastly outperform both, we write Julia code faster with fewer errors.
  - thejosh 4y ago
    What kind of questions are you trying to answer? What's the area for the long computations you're doing? Just curious :)
- rightbyte 4y ago
  You need to know how the compiler propagates types in detail to write performant code. It is quite hard.
  - CyberDildonics 4y ago
    It's not 'quite hard' you just need to put types on your variable and function declarations. If you leave them off you can make a function generic but you need to make sure you don't have multiple returns that could each return different types.
tkuraku 4y ago
I think there are a lot of performance pitfalls with Julia. In my experience you have to pay close attention https://docs.julialang.org/en/v1/manual/performance-tips/ https://docs.julialang.org/en/v1/manual/performance-tips/. Otherwise you can easily get lack luster performance. However, there are times when the speed of Julia really shines and it just feels magical compared to R, python, etc.

ChrisRackauckas 4y ago

Did you use JuMP? It would be interesting to see the JuMP code.

BrandonS113 4y ago
No, nlopt. Why we could easy port from R to Julia as nlopt exists for both (its c)
- ChrisRackauckas 4y ago
  Ahh then it's all based on the speed of the objective function and the gradient calculation code. Do you have an example to look at for that? Or did you try modelingtoolkitizing and running an auto-simplified form?
  - BrandonS113 4y ago
    Yes, well no gradients. It a simple parametric function (6 parameters) applied to a vector of lengths 10 to 20. All vectorizes in the language of R. All powers, logs, exponentials, ratios, and sums. Large number of local maxima, and places where function cannot be evaluated. So probabilistic optimizer is very important. Come to think of it, for this sort of calculations, R and Julia should take the same time.
    - DNF2 4y ago
      Of course, this is piques one's curiosity. It might be, if the function is simple enough, that there is little advantage to Julia here. But if you are combining multiple operations on a vector, there could be opportunities for Julia, in-place operations, fusing, simd. Maybe even StaticArrays. Any chance of sharing that little piece of code?
      - BrandonS113 4y ago
        sure, for parameters par (what is optimized for), data vector x (typical length from 10 to 20), constants n and n2, a typical function is if((1 - par[3]^2)<0 ) return(100) if(par[1] + par[2] \* par[5] \* sqrt(1 - par[3]^2)<0) return(500 ) tmp <- (x - par[4])^2 + par[5]^2 tmp2 <- par[1] + par[2] \* (par[3] \* (x - par[4]) + sqrt(tmp)) result <- ((sqrt(tmp2) - Ce) / n2)^2 result <- sqrt(sum(n \* result) / n) if(is.na(result)) return(1111) return(result)
        DNF2 4y ago
        I'm not that good at reading R. But if the Julia code is similar, then this code is type unstable, sometimes returning an Int, sometimes a Float. That harms performance. Generally, it looks like a function where Julia could have a significant performance advantage.
        BrandonS113 4y ago
        R looks pretty much the same as Julia here. Don't think its unstable. Its all float. Nlopt passes parameters as float, data is float. Can't see how it could return int. But I'll take a second look at the julia version of it. We optimised R as much as possible. might have missed a julia trick. simd perhaps. Except we also run it on ARM
        DNF2 4y ago
        Maybe it's all floats in R, but in Julia, `return 500` means an Int is returned. But it's really hard to determine whether the Julia code is unstable based only on the R code.
  - celrod 4y ago
    I've noticed that R often defaults to much higher tolerances than Julia, even when it's wrappers to the same C library, like cubature. R cubature[0]: 1e-5 Cubature.jl [1]: 1e-8 The difference for NLopt in R vs Julia is smaller. `NLopt.DEFAULT_OPTIONS`[2] in Julia shows `1e-7` for `ftol_rel`, `xtol_rel`, and `constrtol_abs`, while in R `xtol_rel` is `1e-6` and the others are `0.0`[3]. So, the options at least aren't the same with nlopt. Anyway, I always recommend confirming that you're comparing the same settings. And of course, in Julia, you'll probably want to `JET.report_opt` your function and fix and glaring performance issues. NLopt seems like it may be a bit of an exception, but I noticed this is pretty common pattern elsewhere, uniroot[4] being another example, with eps()^(1/4) default tolerance, far higher than Julia root solvers will use. [0] https://cran.r-project.org/web/packages/cubature/cubature.pdf https://cran.r-project.org/web/packages/cubature/cubature.pd... [1] https://github.com/JuliaMath/Cubature.jl https://github.com/JuliaMath/Cubature.jl [2] https://github.com/JuliaOpt/NLopt.jl/blob/6ade25740362895bbfff1aee07d35911a6e2df17/src/MOI_wrapper.jl#L315 https://github.com/JuliaOpt/NLopt.jl/blob/6ade25740362895bbf... [3] https://cran.r-project.org/web/packages/nloptr/nloptr.pdf https://cran.r-project.org/web/packages/nloptr/nloptr.pdf [4] https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/uniroot https://www.rdocumentation.org/packages/stats/versions/3.6.2...
    - BrandonS113 4y ago
      We take of that by setting all the NLOPT options to be the same across calling languages. We get pretty much the same number of calls in julia, r and c.
- harshreality 4y ago
  If you use the same native C library from both Julia and some other language, and (presumably) that library is where the algorithm spends the vast majority of its time, why would you expect Julia to be faster?
  - BrandonS113 4y ago
    no. the code spends 99% of its time in 5 lines in the objective function. And my usual experience is that Julia is much faster than R. Just not always, apparently.

Lyngbakr 4y ago

I only used Julia for a short time, but I didn't see the blazing fast speeds I was promised. I've seen the benchmarks, of course, on which the claims are founded, but the C-like speeds weren't obvious to me in everyday data science workflows. In the end, there wasn't sufficient motivation for me to switch to Julia as my weapon of choice. I do like Pluto[0], though... [0]https://plutojl.org/ https://plutojl.org/

BrandonS113 4y ago
If doing data science, I find Julia's tools to be inferior to Python and R. But in my work, when it comes to long computations, not only does Julia usually vastly outperform both, we write Julia code faster with fewer errors.
- thejosh 4y ago
  What kind of questions are you trying to answer? What's the area for the long computations you're doing? Just curious :)
rightbyte 4y ago
You need to know how the compiler propagates types in detail to write performant code. It is quite hard.
- CyberDildonics 4y ago
  It's not 'quite hard' you just need to put types on your variable and function declarations. If you leave them off you can make a function generic but you need to make sure you don't have multiple returns that could each return different types.

tkuraku 4y ago

I think there are a lot of performance pitfalls with Julia. In my experience you have to pay close attention https://docs.julialang.org/en/v1/manual/performance-tips/ https://docs.julialang.org/en/v1/manual/performance-tips/. Otherwise you can easily get lack luster performance. However, there are times when the speed of Julia really shines and it just feels magical compared to R, python, etc.