Is it better to compute Average Precision using the trapezoidal rule or the rectangle method? Background
Average precision is a popular and important performance metric widely used for, e.g., retrieval and detection tasks. It measures the area under the precision-recall curve, which plots the precision values for all possible detection thresholds against the respective recall values.
A typical precision-recall curve would look somewhat like this.

Confusion
Since there is only a limited number of thresholds to be evaluated on any given finite dataset, the precision-recall curve is a piece-wise linear function. Thus, I always thought the correct way of computing the area under that curve would be using the trapezoidal rule:
ap = sum( (recall[k+1] - recall[k]) * (precision[k+1] - precision[k]) / 2 )

Not only I thought so, apparently, since the official code for evaluating results on the Oxford Buildings dataset, which is a widely used benchmark in content-based image retrieval, computes AP in exactly this way, too.
However, the scikit-learn package for Python computes average precision differently in the function sklearn.metrics.average_precision_score, following the definition of AP from Wikipedia. They use the rectangle method:
ap = sum( (recall[k+1] - recall[k]) * precision[k+1] )

In the example given above, that would approximate the area under the precision-recall curve with the red function in the following figure:

The documentation of scikit-learn says the following about this:

This implementation is different from computing the area under the precision-recall curve with the trapezoidal rule, which uses linear interpolation and can be too optimistic.

One of the curators of the Oxford Buildings dataset, on the other hand, explained in a related question on StackOverflow that the rectangle method would be a "commonly used worse approximation".
Question
Bad enough that different benchmarks and different packages use different variants of Average Precision to compare methods, but now I wonder:

*

*Which of the two versions is the "better" way of doing it? Trapezoidal rule or rectangle method?

*What are the pros and cons for each?

*What does the scikit-learn documentation mean with the claim the the trapezoidal rule is "too optimistic"?

 A: For a piecewise linear function, using the trapezoidal rule with endpoints on each of the ends of the "pieces" will yield the exact area under the curve --- i.e., it is equivalent to integration under the curve.  This occurs when the trapezoids correspond exactly with the lines in the piecewise linear function.  (Of course, this does not hold if there are endpoints of the pieces in the function that are not endpoints of the trapezoids.)  Contrarily, the rectangular method will not give an exact area under the curve, though it should be close if you use a large number of rectangles.
As to which method is better, the exact method (trapezoidal) is better if it is computationally feasible.  I am not aware of any particular reason why it should not be computationally more expensive than the rectangular method, since the only difference is that it uses the average height of each endpoint instead of the maximum height.  If we partition the recall values using the endpoints $r_0 < r_1 < \cdots < r_n$ then we have:
$$\begin{align}
\text{Rectangular area} 
&= \sum_{k=1}^n (r_k - r_{k-1}) \times \max (P(r_k), P(r_{k-1})), \\[10pt]
\text{Trapezoidal area} 
&= \sum_{k=1}^n (r_k - r_{k-1}) \times \frac{P(r_k) + P(r_{k-1})}{2}. \\[6pt]
\end{align}$$
Assuming that these endpoints contain the endpoints of the piecewise linear function, it is simple to show that the trapezoidal area is the exact area under the curve.
A: UPDATE
This journal article explains why linear interpolation is both "too optimist" and is also "incorrect," due to the non-linear properties of the Precision-Recall curve: https://www.biostat.wisc.edu/~page/rocpr.pdf

*

*Due to the non-linear nature of the precision-response curve, linear interpolation results in erroneous overestimations.

*With an averaging rule, missed changes in slope average out. With interpolation, they do not. Thus if points don't cover all the slope changes in the real curve, interpolation errors add up.  This is why interpolation is "too optimistic," and why the midpoint rule generally has half the error that the trapezoidal rule has.
https://math.libretexts.org/Courses/Mount_Royal_University/MATH_2200%3A_Calculus_for_Scientists_II/2%3A_Techniques_of_Integration/2.5%3A_Numerical_Integration_-_Midpoint%2C_Trapezoid%2C_Simpson%27s_rule
http://math.cmu.edu/~mittal/Recitation_notes.pdf
https://activecalculus.org/single/sec-5-6-num-int.html
This article is referenced as "[Davis2006]" in the scikit-learn documentation as the explanation as to why linear interpolation is inappropriate and "too optimistic" here.  See: https://scikit-learn.org/stable/modules/model_evaluation.html#precision-recall-f-measure-metrics
Also,
The function sklearn.metrics.average_precision_score does not use the rectangle rule, or any Riemann sum, right or otherwise.  It uses "average precision." The formulas are very different.

Note that f(x) is very, very different than Pi.
Due to the formulas for precision and recall, the Average Precision is actually computing an average, with discrete values between 0 and 1. Regarding Riemann, f(x) = y.  This gives you the height to multiply the delta with.  There is no averaging there.
Average precision is most analogous to the midpoint rule, as they are both doing averages.
Note that R uses the same formula for Average Precision: https://www.rdocumentation.org/packages/yardstick/versions/0.0.4/topics/average_precision
