By: The Julia Language

Re-posted from: https://www.youtube.com/watch?v=H1NRkUlCtJo

By: Tamás K. Papp

Re-posted from: https://tamaspapp.eu/post/log1p/

edit: fixed bogus interaction of MathJax and code highlighting.

edit2: added benchmarks.

This is a follow-up on a question I asked on the Julia forums about calculating

\[
\text{log1p}(x) = \log(1+x)
\]

This calculation is tricky because if \(x \approx 0\),

\[
\log(1+x) \approx x
\]

while as \(x \to \infty\), \(\log(1+x)\) approaches \(\log(x)\), so simply using log(1+x) will not be as accurate as it could be. Numerical analysts have developed specialized methods for these calculations that try to get an accurate answer, and all programming languages serious about numerical calculations have an implementation either in the core language or a library.

Julia’s Base.log1p currently suggests that Base.Math.JuliaLibm.log1p would be preferable, but then I was wondering why isn’t that the default? So I decided to perform a trivial numerical experiment, calculating the error for both, and also benchmark the two methods.

Accuracy

The key question is what to compare the results with. One could compare to an existing "gold standard" implementation, or simply calculate the results using a higher precision representation. Fortunately, Julia has BigFloats available right out of the box.

The graph below shows the base-2 logarithm of the relative error for Base.log1p vs \(\log\_2(1+x)\); horizontal lines are log2(eps()) and log2(eps())+1. This suggests that Base.log1p is very accurate, but not as good as it could be when \(x \approx 0\).

The next plot shows the relative accuracy of the relative error above, comparing Base.Math.JuliaLibm.log1p to Base.log1p (lower values better). In these simulations, Base.Math.JuliaLibm.log1p is never worse, but sometimes significantly better (resulting in an extra binary digit of accuracy). This matters especially when \(x \approx 0\).

The next plot confirms this.

Speed

I also evaluated relative speed. The graph below shows the relative runtimes, obtained using BenchmarkTools.@belapsed. Values below \(1\) mean that Base.Math.JuliaLibm.log1p is faster: indeed, this seems to be the case except for values very close to \(0\), where there is a 10–20% performance penalty. At other values, Base.Math.JuliaLibm.log1p is 30–40% faster.

Conclusion

1. For values near \(0\), Base.Math.JuliaLibm.log1p is more accurate, at a slight performance cost.

2. For values away from \(0\), it is at least as accurate as Base.log1p, and 30—40% faster.

To me, this suggests that Base.Math.JuliaLibm.log1p should be the default method — mostly because the extra accuracy is more important to me than the slight performance cost.

Code is available below.

# consistent random numbers
srand(UInt32[0xfd909253, 0x7859c364, 0x7cd42419, 0x4c06a3b6])

"""
    err(x, [prec])

Return two values, which are the log2 relative errors for calculating
`log1p(x)`, using `Base.log1p` and `Base.Math.JuliaLibm.log1p`.

The errors are calculated by compating to `BigFloat` calculations with the given
precision `prec`.
"""
function err(x, prec = 1024)
    yb = log(1+BigFloat(x, prec))
    e2(y) = Float64(log2(abs(y-yb)/abs(yb)))
    e2(log1p(x)), e2(Base.Math.JuliaLibm.log1p(x))
end

z = exp.(randn(20000)*10)       # z > 0, Lognormal(0, 10)
x = z .- 1                      # x > -1
es = map(err, x)                # errors

using Plots; gr()               # plots

scatter(log2.(z), first.(es), xlab = "log2(x+1)", ylab = "log2 error of Base.log1p",
        legend = false)
hline!(log2(eps())-[0,1])
Plots.svg("Base_log1p_error.svg")
scatter(log2.(z), last.(es) .- first.(es), xlab = "log2(x+1)",
        ylab = "improvement (Base.Math.JuliaLibm.log1p)", legend = false)
Plots.svg("JuliaLibm_improvement.svg")
scatter(log2.(z), last.(es), xlab = "log2(x+1)",
        ylab = "log2 error of Base.Math.JuliaLibm.log1p", legend = false)
hline!(log2(eps())-[0,1])
Plots.svg("JuliaLibm_log1p_error.svg")

######################################################################
# WARNING: these run for a very long time
######################################################################
using BenchmarkTools

z = exp.(vcat(randn(200)*10, rand(200)*0.1)) # z > 0, more values around 
x = z .- 1                                   # x > -1
b1 = [@belapsed log1p($x) for x in x]        # WARNING: takes forever
b2 = [@belapsed Base.Math.JuliaLibm.log1p($x) for x in x] # ditto

scatter(log2.(z), b2 ./ b1, xlab = "log2(x+1)",
        ylab = "time Math.JuliaLibm.log1p / log1p", legend = false, yticks = 0:0.2:1.2)
hline!([1])
Plots.svg("relative_time.svg")

Berkeley, CA – Julia has joined the rarefied ranks of computing languages that have achieved peak performance exceeding one petaflop per second – the so-called ‘Petaflop Club.’

The Julia application that achieved this milestone is called Celeste. It was developed by a team of astronomers, physicists, computer engineers and statisticians from UC Berkeley, Lawrence Berkeley National Laboratory, National Energy Research Scientific Computing Center (NERSC), Intel, Julia Computing and the Julia Lab at MIT.

Celeste uses the Sloan Digital Sky Survey (SDSS), a dataset of astronomical images from the Apache Point Observatory in New Mexico that includes every visible object from over 35% of the sky – hundreds of millions of stars and galaxies. Light from the most distant of these galaxies has been traveling for billions of years and lets us see how the universe appeared in the distant past.

Since SDSS data collection began in 1998, the process of cataloging these stars and galaxies was painstaking and laborious.

So the Celeste team developed a new parallel computing method to process the entire SDSS dataset. Celeste is written entirely in Julia, and the Celeste team loaded an aggregate of 178 terabytes of image data to produce the most accurate catalog of 188 million astronomical objects in just 14.6 minutes with state-of-the-art point and uncertainty estimates.

Celeste achieved peak performance of 1.54 petaflops using 1.3 million threads on 9,300 Knights Landing (KNL) nodes of the Cori supercomputer at NERSC – a performance improvement of 1,000x in single-threaded execution.

The Celeste research team is already looking to new challenges. For example, the Large Synoptic Survey Telescope (LSST), scheduled to begin operation in 2019, is 14 times larger than the Apache Point telescope and will produce 15 terabytes of images every night. This means that every few days, the LSST will produce more visual data than the Apache Point telescope has produced in 20 years. With Julia and the Cori supercomputer, the Celeste team can analyze and catalog every object in those nightly images in as little as 5 minutes.

The Celeste team is also working to:

• Further increase the precision of point and uncertainty estimates
• Identify ever-fainter points of light near the detection limit
• Improve the quality of native code for high performance computing

The Celeste project is a shining example of:

• High performance computing applied to real-world problems
• Cross-institutional collaboration including researchers from UC Berkeley, Lawrence Berkeley National Laboratory, National Energy Research Scientific Computing Center (NERSC), Intel, Julia Computing and the Julia Lab at MIT
• Cross-departmental collaboration including astronomy, physics, computer science, engineering and mathematics
• Julia, the fastest modern open source high performance programming language for scientific computing
• Parallel and multithreading supercomputing capabilities
• Public support for basic and applied scientific research

About Julia and Julia Computing

Julia is the fastest modern high performance open source computing language for data, analytics, algorithmic trading, machine learning and artificial intelligence. Julia combines the functionality and ease of use of Python, R, Matlab, SAS and Stata with the speed of C++ and Java. Julia delivers dramatic improvements in simplicity, speed, capacity and productivity. Julia provides parallel computing capabilities out of the box and unlimited scalability with minimal effort. With more than 1 million downloads and +161% annual growth, Julia is one of the top 10 programming languages developed on GitHub and adoption is growing rapidly in finance, insurance, energy, robotics, genomics, aerospace and many other fields.

Julia users, partners and employers hiring Julia programmers in 2017 include Amazon, Apple, BlackRock, Capital One, Comcast, Disney, Facebook, Ford, Google, Grindr, IBM, Intel, KPMG, Microsoft, NASA, Oracle, PwC, Raytheon and Uber.

1. Julia is lightning fast. Julia provides speed improvements up to
   1,000x for insurance model estimation, 225x for parallel
   supercomputing image analysis and 11x for macroeconomic modeling.

2. Julia is easy to learn. Julia’s flexible syntax is familiar and
   comfortable for users of Python, R and Matlab.

3. Julia provides unlimited scalability. Julia applications can be deployed on large clusters with a click of a button and can run parallel and distributed computing quickly and easily on tens of thousands of nodes.

4. Julia integrates well with existing code and platforms. Users of
   Python, R, Matlab and other languages can easily integrate their
   existing code into Julia.

5. Elegant code. Julia was built from the ground up for
   mathematical, scientific and statistical computing, and has advanced
   libraries that make coding simple and fast, and dramatically reduce
   the number of lines of code required – in some cases, by 90%
   or more.

6. Julia solves the two language problem. Because Julia combines
   the ease of use and familiar syntax of Python, R and Matlab with the
   speed of C, C++ or Java, programmers no longer need to estimate
   models in one language and reproduce them in a faster
   production language. This saves time and reduces error and cost.

Julia Computing was founded in 2015 by the creators of the open source Julia language to develop products and provide support for businesses and researchers who use Julia.