Used gprof with grpof2dot.py to get the following diagram: gprof.pdf.

iso_mvn_lpdf is getting hit hard:

    double iso_mvn_lpdf(const double mu, const double y, double epsilon, size_t D)
    {
        double accum = 0;
        double d;
        for(size_t i = 0; i < D; ++i)
        {
            d = (mu++ - y++) / epsilon;
            accum += dd;
        }
        return -0.5  (accum + Dlog(2M_PIepsilonepsilon));
    }

It's already pretty lean (no alloc, all c-style). But we can move the divide-by-epsilon out of the loop for an easy 1.8x speedup.

    double iso_mvn_lpdf(const double mu, const double y, double epsilon, size_t D)
    {
        double accum = 0;
        double d;
        for(size_t i = 0; i < D; ++i)
        {
            d = mu++ - y++;
            accum += dd;
        }
        accum /= epsilonepsilon;
        return -0.5  (accum + Dlog(2M_PIepsilon*epsilon));
    }

Now the bottleneck is all the allocation, copying and freeing of kjb::Vector temporaries.

I tweaked the code for evaluating a piecewise linear function to avoid creating kjb::Vector temporaries, and running time dropped dramatically from 11.4s to 1.26s. This is in production mode, so its surprising that more temporaries aren't optimized out.

GProf after eliminating temporaries: gprof.pdf.

Since 1 iterations takes about 1.2s, bumping up to 10 iterations.

Remaining speed-up opportuntiies: exploit gradient independence, parallel gradient

Parallel gradient

Enabled 8-way parallel gradient evaluation, and got worse performance!

single threaded:  0:09.63
multi threaded: 0:15.48

This despote top displaying 550% CPU utilization.

Maybe gnuprof is affecting performance

Tuning

step size gradient size