All packages for numerical math

[Hans W Borchers]

I would like to use Julia mostly on applied numerical math problems. I think what I need is an overview of what is available in Base and contributed packages. Therefore I decided to write a short "vignette" listing all relevant functions and packages, possibly with a short description, usage, and one or two examples. 

I have identified the items in the following list. I would welcome any suggestions of functions I have overlooked in this area.

(Optimization at the moment is not covered by intention.)

### Julia Base  0.2.0

    -------------------   --------------------------------------------------

    sin, cos, ...         Trigonometric, hyperbolic functions, square root,

                              exponential and logarith functions, etc.

    factor, gcd, primes   Number-theoretic and combinatorial functions

    factorial, binomial

    ...

    gamma, airy, bessel   Special functions

    beta, eta, zeta

    ...

    quadgk                Gauss-Kronrod adaptive integration

### Calculus  0.1.3

    -------------------   --------------------------------------------------

    derivative, ...       finite-differences numerical derivatives

    gradient, hessian         second-order, gradients and hessians

    [unfortunately:]      finite_difference and complex_step not exported!

    integrate             adaptive Simpson, Monte-Carlo integration

    differentiate         symbolic differentiation

### DualNumbers  0.0.0 (??)

    -------------------   --------------------------------------------------

    Dual, epsilon         Automated Differentiation

                              [if the function accepts dual numbers]

### Cubature  1.0.1

    -------------------   --------------------------------------------------

    h/pquadrature         one- and multidimensional adaptive integration

    h/pcubature               (Gauss-Kronrod, Genz-Malik, Clenshaw-Curtis)

### Grid  0.2.8

    -------------------   --------------------------------------------------

    InterpGrid            function interpolation on (regular) grids

                          [a simpler interp1d() function for irregular

                           grids might still be helpful]

### BSplines  0.0.0 (??)

    -------------------   --------------------------------------------------

    linear/quadratic/     

        cubicSplineBFE

### ApproxFun  0.0.0 (??)

    -------------------   --------------------------------------------------

    Fun                   generates an approximating function

    diff, cumsum          differentiate or integrate the approximation

### Polynomial  0.0.0 (??)

    -------------------   --------------------------------------------------

    Poly, poly            construct polynomial from coefficients resp. roots

    polyval               evaluate the polynomial (Horner scheme)

    polyint, polyder      derivative, anti-derivative of a polynomial

    roots                 determine the roots/zeros of the polynomial 

                              (as eigenvalues of the companion matrix)

                          [missing is a polyfit() function]

### PowerSeries

    -------------------   --------------------------------------------------

    series, restrict      generates or truncates a (finite) power series

                              (used to determine the Taylor series)

### Roots  0.1.0

    -------------------   --------------------------------------------------

    find_zero, fzero      bracketing/derivative-free root finding methods

                              [Ridders' method is not implemented !]

    newton, halley        derivative-based root finding methods

    multroot              multiple roots of (inexact) polynomials

### NLsolve  0.1.1

    -------------------   --------------------------------------------------

    nlsolve               solves systems of nonlinear equations, based on

                              Newton and trust-region approaches

### ODE  0.0.0 (??)

    -------------------   --------------------------------------------------

    ode23                 Runge-Kutta (2, 3)-method with variable step size

    ode4                  Runge-Kutta of order 4 (with fixed step size?)

### Elliptic  0.2.0

    -------------------   --------------------------------------------------

    K, F, E, Pi           (in)complete elliptic integrals of 1. and 2. kind

### GSL  0.1.1

    -------------------   --------------------------------------------------

                          interface to the GNU Scientific Library(GSL)

    e.g., hypergeom       Gauss' hypergeometric function 2F1

[Tim Holy]

On Monday, April 21, 2014 01:20:39 AM Hans W Borchers wrote:
> [a simpler interp1d() function for irregular
> grids might still be helpful]

It's actually there, I just need to documented it.

--Tim

[Tomas Lycken]

ODE.jl has released a v0.1.0, that contained also

* ode45, a (4,5)-order adaptive method (there's also a few unexported methods with a few choices for the RK coefficients, e.g. Dormand-Prince, Fehlberg or Cash-Karp coeffs)
* ode4s, a 4th order solver for stiff problems (based on the Rosenbrock method, also with a few unexported choices for the coeffients)
* ode4ms, a 4th order fixed-step multistep method with Adams-Bashforth-Moulton coefficients

However, the API for ODE.jl is under discussion, and although development happens in irregular bursts of activity, the last of which is a while ago, it is worthwhile to keep that in mind when using the library. There *will* be breaking changes in the next release, whenever that happens.

// Tomas

[Robert J Goedman]

Hi Hans,

Two additional packages are TaylorSeries (not in METADATA right now) and HyperDualNumbers, which I 'derived' from the C code by the authors (see below links).

The entry in your list could look like:

### HyperDualNumbers  0.1.1

    -------------------   --------------------------------------------------

    Hyper, eps1,eps1eps2  Automated 1st & 2nd order Differentiation

                              [if the function accepts dual numbers]

The c++ code can be found at http://adl.stanford.edu/hyperdual/hyperdual.h . The paper itself at http://adl.stanford.edu/hyperdual/Fike_AIAA-2011-886.pdf .

The Julia package can be found at: https://github.com/goedman/HyperDualNumbers.jl.git .

Regards,

Rob J. Goedman

goe...@icloud.com

[Evgeny Shevchenko]

Hi 

Is there a package for numeric integration of obtained arrays, say x and y? Couldn't find one, the led me to use @pyimport and numpy.trapz.

--

pupoque@IRC

[John Myles White]

Have you tried the quadgk function?

 -- John

[Евгений Шевченко]

Hi, John. 

No, I didn't. I didn't find it and it seems to be not what i need:

"no method quadgk(Array{Float64,1}, Array{Float64,1})"

quadgk(f,a,b,...) expects a function as its first argument but I mean the case when y = f(x), but i don't have f, e.g. obtained experimental data, so x and y are 1-D arrays of floats.

[Tomas Lycken]

The trapezoidal rule (http://en.wikipedia.org/wiki/Trapezoidal_rule) would probably be almost trivial to implement.

function trapz{T<:Real}(x::Vector{T}, y::Vector{T})

   if (length(y) != length(x))

       error("Vectors must be of same length")

   end

   sum( (x[2:end] .- x[1:end-1]).*(y[2:end].+y[1:end-1]) ) / 2

end

    
x = [0:0.01:pi]
y = sin(x)

trapz(x,y) # 1.9999820650436642

This, of course, only currently works on vectors of real numbers, but it's easy to extend it if you want.

And there might be more accurate methods as well, of course (see e.g. http://en.wikipedia.org/wiki/Simpson%27s_rule) but this one's often "good enough".

// T

[Cameron McBride]

Or you can use the non-vectorized version and save the overhead of the temporary arrays being created by the addition and multiplication steps.

function trapz{T<:Real}(x::Vector{T}, y::Vector{T})

    local len = length(y)

    if (len != length(x))

        error("Vectors must be of same length")

    end

    r = 0.0

    for i in 2:len

        r += (x[i] - x[i-1]) * (y[i] + y[i-1])

    end

    r/2.0

end

BTW, another possibility is to use a spline interpolation on the original data and integrate the spline evaluation  with quadgk().  (Ideally, integrations can be incorporated into a spline interface.) This could be useful depending on the functional shape between the grid points. 

Cameron

[Tomas Lycken]

On Wednesday, April 23, 2014 11:10:15 PM UTC+2, Cameron McBride wrote:

Or you can use the non-vectorized version and save the overhead of the temporary arrays being created by the addition and multiplication steps.

There's really no way I can hide that I learnt scientific computing in Matlab, is there? :P

[Hans W Borchers]

Just being curious, I compared these two version on accuracy and execution time, having grown up in Matlab as well. Let's call the first, vectorized version by Tomas as 'trapz1', the second version by Cameron as 'trapz2', see below. Then I get the following results:

    julia> x = linspace(0, pi, 1000000);

    julia> y = sin(x);

    julia> trapz1(x, y)
    1.999999999998355

    julia> trapz2(x, y)
    1.9999999999984162

    julia> @time for i in 1:100; trapz1(x, y); end
    elapsed time: 4.858258575 seconds (5600105600 bytes allocated)

    julia> @time for i in 1:100; trapz2(x, y); end
    elapsed time: 0.41866061 seconds (1600 bytes allocated)

and in general the unvectorized version is more than 10 times faster. But why are the two results different (though only by 6.1e-14)? It looks like they are performing exactly the same computations. Does the sequence of elementary operations play a role?

Is it to be expected that a version with a localized length is faster? I didn't see significantly different running times.

PS: The function definitions I used here are:

    function trapz1{T<:Number}(x::Vector{T}, y::Vector{T})
        local n = length(x)
        if length(y) != n
            error("Vectors 'x', 'y' must be of same length.")
        end
        if n == 1; return 0.0; end
        sum((x[2:end]-x[1:end-1]) .* (y[2:end]+y[1:end-1])) / 2.0
    end


    function trapz2{T<:Number}(x::Vector{T}, y::Vector{T})
        local n = length(x)
        if (length(y) != n)
            error("Vectors 'x', 'y' must be of same length")
        end
        if n == 1; return 0.0; end
        r = 0.0
        for i in 2:n
            r += (x[i] - x[i-1]) * (y[i] + y[i-1])

        end
        r / 2.0
    end

[Tobias Knopp]

In the vectorized version the data is written back to memory while the unvectorized version can hold everything in the CPU registers. This can explain slightly different results.

[Alex]

I couldn't find the definition of sum right now, but I guess it is calling BLAS or something like that, which might explain the difference?

Consider for example

function trapz3{T<:Number}(x::Vector{T}, y::Vector{T})

    local n = length(x)
    if length(y) != n
        error("Vectors 'x', 'y' must be of same length.")
    end
    if n == 1; return 0.0; end

    int = (x[2:end]-x[1:end-1]) .* (y[2:end]+y[1:end-1])

    r = 0.0
    for i in 1:length(int)
      r += int[i]
    end
    r / 2.0, BLAS.asum(int)/2.0, sum(int)/2.0
end

which gives me

julia> trapz3(x, y)
(1.9999999999984162,1.9999999999983549,1.999999999998355)

Best,

Alex.

PS: Just a minor remark: one should probably use zero(T) instead of 0.0 (or even zero(zero(T)*zero(T))).

[Tobias Knopp]

The usage of construct like zero(zero(T)*zero(T)) is not really easy to understand. I assume you propose this due to the upcasting arithmetic integer behavior of non-native integer types (int16, ...), right?

[Alex]

The usage of construct like zero(zero(T)*zero(T)) is not really easy to understand. I assume you propose this due to the upcasting arithmetic integer behavior of non-native integer types (int16, ...), right?

Basically, I wanted to avoid that r changes its type in the loop (this might be more evident/relevant if the elements of x and y have different types, say Real and Complex). I don't know if it matters much in this case, though.

Best,

Alex.

[Gunnar Farnebäck]

Sum uses pairwise summation, see
https://github.com/JuliaLang/julia/pull/4039

To find the definition:

julia> @which sum(int)
sum{T}(a::AbstractArray{T,N}) at reduce.jl:276

[Tobias Knopp]

Ah yes indeed that is a very good reason when the function is parametric in two different types. Maybe it would be more elgegant to have a zero version that takes multiple arguments and calls promote internally.

[Tomas Lycken]

I wanted to try to compare these versions with a version based on trapz1 but using the @devec macro from Devectorize.jl. However, it seems that development of Devectorize.jl has sort of stopped, and it's REQUIRE specifies Julia 0.2-, so I guess that package isn't really the way to go anymore (I'm on the current master...)

Is there a replacement for it somewhere? If so, where?

// T

[Ivar Nesje]

Have you tried to Pkg.add("Devectorize"). It seems to work for me, and Julia assumes forward compatibility for packages (even though the rate deprecation makes that assumption kind of wrong), so the "0.2-" version simply says it should not be installed on `0.1`. I think the development rate reflects that the low hanging fruits have been picked and nobody seems to have ideas for how to solve the more difficult problems.

It seems like nobody has fixed the deprecation for the `array[2:]` yet, so it does not look like it is used a lot on the master branch.

WARNING: deprecated syntax "x[i:]" at /Users/ivarne/.julia/Devectorize/src/texpr.jl:493.
Use "x[i:end]" instead.


Ivar
 

[Cameron McBride]

On Thu, Apr 24, 2014 at 4:28 AM, Hans W Borchers <hwbor...@gmail.com> wrote:

    function trapz2{T<:Number}(x::Vector{T}, y::Vector{T})
        local n = length(x)
        if (length(y) != n)
            error("Vectors 'x', 'y' must be of same length")
        end
        if n == 1; return 0.0; end
        r = 0.0
        for i in 2:n
            r += (x[i] - x[i-1]) * (y[i] + y[i-1])

        end
        r / 2.0
    end

I'm not sure it's behavior we should "rely on", but the if branch for n == 1 isn't necessary in this for loop, although perhaps it was included to aid the comparison. A range of 2:1 is a zero-length iteration, so the loop will not run even once.  Try this: 

julia> [2:1]

0-element Array{Int64,1}

If  we are being pedantic on types, then the result of the integral should at least be a floating point.  Granted, I am not sure of what the use case for integer vector inputs would be, but I doubt *if* someone used them that they'd want an integer result in return.  (This would be similar to sqrt(), for example.)

Cameron

[Hans W Borchers]

Understood about the "if"; I often used this extra statement because of differences between Matlab and R.

The function as is returns a floating point number also for integer input 'x', 'y'.

But I am more worried about the case that one of the vectors is real, the other complex.

Do I need two types T1 and T2 for defining different types for the input vectors?

[Tomas Lycken]

Yes; if the vectors aren't of the same type, this method isn't applicable. (In other words, you need to have two type arguments in order to support one vector of ints and one vector of floating points). In this case, though, it's probably reasonable to expect only floating point arguments (of the same size), at least as long as they're real.

And as soon as you start working with complex analysis, I'm not entirely sure the trapezoidal rule is valid at all. It might just be because the article author was lazy, but the Wikipedia article only talks about integrals of real-valued functions of one (real, scalar) variable. If you need complex numbers for something more than curiosity about Julia's type system, you probably want another approach altogether...

// Tomas

[Tobias Knopp]

One could also remove all types and only rely on duck typing. Then one will want to initialize r with something like zero(zero(x[1])+zero(y[1])). And the same in the if statement. This function would run at the same. I still like to specify types as this allows me to restrict input parameters.

[Jason Merrill]

On Thursday, April 24, 2014 11:57:33 PM UTC-7, Tomas Lycken wrote:

And as soon as you start working with complex analysis, I'm not entirely sure the trapezoidal rule is valid at all. It might just be because the article author was lazy, but the Wikipedia article only talks about integrals of real-valued functions of one (real, scalar) variable. If you need complex numbers for something more than curiosity about Julia's type system, you probably want another approach altogether...

The trapezoid rule is valid in complex analysis, and in fact, it converges exponentially for integrals around a circular contour for functions that are analytic in an annulus containing the contour. Trefethen has a beautiful paper on this subject:

https://people.maths.ox.ac.uk/trefethen/trefethen_weideman.pdf

This exponential convergence also applies for periodic functions integrated over a full period on the real line, for the same reasons.

[Hans W Borchers]

Exactly. That's what I do quite often with complex line integrals, where one vector is real: a [0,1]-parametrization of the curve, and the other complex: the value of a complex function on the curve. And indeed, this works very well for closed circle integrals in the complex plane, for example returning the residuum of a meromorphic function.

[Luis Benet]

Hi Hans,

Two additional packages are TaylorSeries (not in METADATA right now) and HyperDualNumbers, which I 'derived' from the C code by the authors (see below links).

The link where TaylorSeries.jl can be found is  https://github.com/lbenet/TaylorSeries.jl .

It implements handling polynomial expansions (type `Taylor`) of basic functions in one variable, including differentiation and integration. We are extending it to more variables (type `TaylorN`), which is not yet in github. We use extensively automatic differentiation techniques to manage everything. We plan to make a pull-request to METADATA, but want to have certain functionality first.

Best,

Luis

[Hans W Borchers]

Sorry for reviving an old thread.

I got interested in timings for vectorized vs. unvectorized versions of
functions in MATLAB. So I wrote MATLAB versions of  the `trapz` functions
mentioned in this thread, and measured their running times by applying the
`timeit` macro several hundred times. The results did surprise me:

    trapz1  (vectorized)  :  ~ 35 µs
    trapz2  (unvectorized):  ~ 12 µs  -- per loop

for evaluating the trapezoidal rule on `x=linspace(0,pi); y=sin(x)`.

I know one should not draw conclusions from such simple benchmarks. Still, does
this comparison imply the old myth to "vectorize MATLAB functions to make them
faster" is not true any more?

The JIT compiler that MATLAB is equipped with is not as powerful as Julia's
LLVM-based compiler, but it seems to be at work. Have others had similar
experiences, or am I on a wrong path here?

On Wednesday, April 23, 2014 11:20:50 PM UTC+2, Tomas Lycken wrote:

[Tim Holy]

Quite interesting indeed. Thanks for sharing.
--Tim

On Tuesday, September 30, 2014 04:52:22 AM Hans W Borchers wrote:
> Sorry for reviving an old thread.
>
> I got interested in timings for vectorized vs. unvectorized versions of
> functions in MATLAB. So I wrote MATLAB versions of the `trapz` functions
> mentioned in this thread, and measured their running times by applying the
> `timeit` macro several hundred times. The results did surprise me:
>
> trapz1 (vectorized) : ~ 35 µs
> trapz2 (unvectorized): ~ 12 µs -- per loop
>
> for evaluating the trapezoidal rule on `x=linspace(0,pi); y=sin(x)`.
>
> I know one should not draw conclusions from such simple benchmarks. Still,
> does
> this comparison imply the old myth to "

> *vectorize MATLAB functions to make them faster*" is not true any more?