CamDavidsonPilon/lifelines — pull request #70 by spacecowboy (2014-06-15 14:33:10+00:00)
Title: Make CoxPHFitter linear instead of quadratic

This commit optimizes the algorithm which calculates Efron gradients and
hessian during the Newton-Raphson fitting.

The change basically comes down to noting the fact that:

R(t(i-1)) = [t(i-1)] + R(t(i))

where (the last risk set) R(t(n)) = [t(n)]

The math and original implementation translates to something like:

# Nested sums as in original equations: O(n^2) running time
for i in range(n):
phi_sum = 0
for j in riskset(i):
phi_sum += phi(j)

Whereas taking advantage of the recursive relationship is akin to
something like:

# Use recursive relationship: O(n) running time
phi_sum = 0
for i in reversed(range(n)):
phi_sum += phi(i)

The key is to iterate over the data backwards. For this to work, the
data must be sorted on time (O(nlog(n))). I have compared the old and
new implementations on a real data set (2695,12) and they give the same
results, ignoring normal floating point differences due to order of
summation.

Performance-wise, for data with shape (2695,12), the fitting speed goes
from 6min 27sec to 1.3sec on my machine (9 iterations of N-R).

I am still a bit unfamiliar with Pandas and therefor had some indexing troubles.
That's why the arrays are forced to numpy arrays in _newton_. This could
perhaps be improved.