Suggestions on performance

Question

I have following piece of code that is called by a program many thousands of times as part of Monte-Carlo simulation. using gprof, I see that 38% of the time was spent on this function. Are there any obvious areas of improvement in this code?

double StudentTCopula::uniform_to_default_time_student(double u, const 
std::vector<double>& times, const std::vector<double>& values)
{
  if (u == 0.0)
    return 99999.0;
  if (u == 1.0)
    return 0.0;
  size_t num_points = times.size();
  size_t index = 0;
  for (size_t i{1};i<num_points;++i){
    if (u <= values[i - 1] && u > values[i]){
      index = i;
      break;
    }
  }   
  double tau = 0.0;
  if (index == num_points + 1) {
    auto t1 = times[num_points - 1];
    auto q1 = values[num_points - 1];
    auto t2 = times[num_points];
    auto q2 = values[num_points];
    auto lam = log(q1 / q2) / (t2 - t1);
    tau = t2 - log(u / q2) / lam;
  } else if (index == 0){
    auto t1 = times.back();
    auto q1 = values.back();
    auto t2 = times[index];
    auto q2 = values[index];
    tau = (t1 * log(q2 / u) + t2 * log(u / q1)) / log(q2 / q1);
  } else {
    auto t1 = times[index - 1];
    auto q1 = values[index - 1];
    auto t2 = times[index];
    auto q2 = values[index];
    tau = (t1 * log(q2 / u) + t2 * log(u / q1)) / log(q2 / q1);
 }
 return tau;
}

I have tried removing the intermediate variables and directly genearting a long return expression with the calculations, which did not affect gprof results.

Thank you for your time.

Answer 1

• Comparing floating point representations for equality has its own risks - consider

      if (u <= 0 + EPS)
        return ‹BIG_VALUE›;
      if (1 <= u + EPS)
        return 0;
  

  For suitable values of ‹BIG_VALUE›(there should be a better name) and EPS (possibly 0).

• Most of the handling for index != num_points + 1 is common. Don't repeat yourself:

    } else {
      auto t1 = (index == 0) ? times.back() : times[index - 1];
      auto q1 = (index == 0) ? values.back() : values[index - 1];
      auto t2 = times[index];
      auto q2 = values[index];
      tau = (t1 * log(q2 / u) + t2 * log(u / q1)) / log(q2 / q1);
    }
  

Just before the end of uniform_to_default_time_student(), code indentation is slightly inconsistent.

---

G. Sliepen: use log(a/b)=log(a)−log(b) to avoid the division:

auto log_u = log(u),
    log_q1 = log(q1),
    log_q2 = log(q2);
tau = (t1 * (log_q2 - log_u) + t2 * (log_u - log_q1)) / (log_q2 - log_q1);

(implementation imperfections by greybeard)

---

I have no idea whether substituting two more divisions is useful:

  if (index == num_points + 1) {
    auto t1 = times[num_points - 1];
    auto q1 = values[num_points - 1];
    auto t2 = times[num_points];
    auto _q2 = 1 / values[num_points];
    auto _lam = (t2 - t1) / log(q1 * _q2);
    tau = t2 - log(u * _q2) * _lam;
  }

Answer 2

Disclaimer: Review based on assumptions`

From the code snippet, I am guessing that values are in descending order. In this case you may use std::lower_bound to do a binary search instead of a for loop. Apart from this you can do a micro optimization by replacing log(q2/q1) with log(q2/u) + log(u/q1).

Answer 3

Division is usually a relatively slow operation. You're using the same divisor for a couple of division operations in a few places, so you might gain a little by converting something like this:

tau = (t1 * log(q2 / u) + t2 * log(u / q1)) / log(q2 / q1);

...to something more like this:

double q1i = 1.0 / q1;
tau = t1 * (log(q2 / u) + t2 * log(u * q1i)) / log(q2 * q1i);

May help a little (but no guarantee--your compiler may recognize the same possibility on its own).

At least in my testing, a * log(b) + c * log(d) is slower than just using pow to compute the power, so changing:

double tau1 = (t1 * log(q2 / u) + t2 * log(u / q1)) / log(q2 / q1);

...to:

double tau2 = std::log(std::pow(q2 / u, t1) * std::pow(u / q1, t2)) / log(q2 / q1);

...seems to improve speed by a factor of about 2. I'd guess there's less likelihood of the compiler recognizing the identity involved here, but I suppose it's possible (or it could turn out that with the CPU/flags you're using, the original equation is just as fast; impossible to be sure).

Answer 4

I don't understand this stanza:

  if (index == num_points + 1) {

    auto t1 = times[num_points - 1];
    auto q1 = values[num_points - 1];
    auto t2 = times[num_points];
    auto q2 = values[num_points];
    auto lam = log(q1 / q2) / (t2 - t1);
    tau = t2 - log(u / q2) / lam;
  } else ...

After initializing at zero, we have the invariant that i, and therefore index, shall always be strictly less than the number of points. As @Pkkm observed, a [num_points] de-reference Would Be Bad, if it ran.

Surely the stanza is equivalent to:

  if (0) {
    /* NOTREACHED */
  } else ...

Why did your toolchain not issue "C4702 unreachable code" or a similar warning? Does it lack a -Wunreachable-code switch? (Gcc, alas, removed that switch in 2010, though it does offer -Wunreachable-code-ctrl.)

Recommend deleting it, or changing that conditional.

---

Please assert times.size() == values.size(). Or convert to passing around an array of structs, each holding a time and a value.

---

It wasn't obvious to me how to go from this definition to the OP's calculations. Within the source code, please cite a reference.